{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Fitting RooFit-generated likelihoods\n", "\n", "In this tutorial, we show how to create a negative log-likelihood function with the [RooFit framework](https://root.cern/manual/roofit/) and minimize it with iminuit. \n", "\n", "RooFit is a powerful fitting framework developed by CERN's ROOT team.\n", "RooFit is very powerful and sophisticated, but there are a few reasons to use iminuit instead:\n", "\n", "- RooFit documention is extensive, but lacking in important places\n", "- RooFit interfaces are not Pythonic, they mimic the C++ interfaces (which are also dated)\n", "- Errors are difficult to understand and debug since RooFit is internally executing C++ code\n", "- You may experience segmentation faults when using RooFit from Python due to bugs in the ROOT Python layer (problems with handling life-times of dynamic C++ objects in Python correctly)\n", "\n", "For these reasons, you may consider to transition to iminuit and its cost functions for your project. As a first step, you want to convince yourself that iminuit gives you the same fitting result as you get from RooFit. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Welcome to JupyROOT 6.26/10\n" ] } ], "source": [ "# ROOT is best installed via a conda virtual environment from conda-forge\n", "import ROOT" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# fix PRNG seed for RooFit random number generation\n", "ROOT.RooRandom.randomGenerator().SetSeed(1)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We generate a Gaussian with mean 1 and width 3 and draw 10000 data points from it. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "x = ROOT.RooRealVar(\"x\", \"x\", -10, 10)\n", "mean = ROOT.RooRealVar(\"mean\", \"mean of gaussian\", 1, -10, 10)\n", "sigma = ROOT.RooRealVar(\"sigma\", \"width of gaussian\", 3, 0.1, 10)\n", "\n", "gauss = ROOT.RooGaussian(\"gauss\", \"gaussian PDF\", x, mean, sigma)\n", "\n", "data = gauss.generate({x}, 10000)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We now fit this Gaussian. We use the `createNLL` method and a simple wrapper function `evaluate`. Note that this simple wrapping function does not propagate the parameter names of the Gaussian to iminuit. A future version of iminuit will come with a builtin wrapper that will also propagate the names and limits." ] }, { "cell_type": "code", "execution_count": 4, "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 = 2.514e+04 Nfcn = 31
EDM = 2.72e-08 (Goal: 0.0001)
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", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0 1.003 0.030
1 x1 3.017 0.022
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x0 x1
x0 0.000926 0 (0.030)
x1 0 (0.030) 0.000497
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 2.514e+04 │ Nfcn = 31 │\n", "│ EDM = 2.72e-08 (Goal: 0.0001) │ │\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 │ 1.003 │ 0.030 │ │ │ │ │ │\n", "│ 1 │ x1 │ 3.017 │ 0.022 │ │ │ │ │ │\n", "└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌────┬───────────────────┐\n", "│ │ x0 x1 │\n", "├────┼───────────────────┤\n", "│ x0 │ 0.000926 0 │\n", "│ x1 │ 0 0.000497 │\n", "└────┴───────────────────┘" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from iminuit import Minuit\n", "\n", "nll = gauss.createNLL(data)\n", "\n", "def evaluate(*args):\n", " for par, arg in zip(nll.getVariables(), args):\n", " par.setVal(arg)\n", " # following RooMinimizerFcn.cxx \n", " nll.setHideOffset(False)\n", " r = nll.getVal()\n", " nll.setHideOffset(True)\n", " return r\n", "\n", "evaluate.errordef = Minuit.LIKELIHOOD\n", "\n", "m = Minuit(evaluate, *[p.getVal() for p in nll.getVariables()])\n", "m.migrad()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Let's compare this to fitting directly with the `fitTo` method." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization\n", " **********\n", " ** 1 **SET PRINT 1\n", " **********\n", " **********\n", " ** 2 **SET NOGRAD\n", " **********\n", " PARAMETER DEFINITIONS:\n", " NO. NAME VALUE STEP SIZE LIMITS\n", " 1 mean 1.00653e+00 2.00000e+00 -1.00000e+01 1.00000e+01\n", " 2 sigma 3.01930e+00 9.90000e-01 1.00000e-01 1.00000e+01\n", " **********\n", " ** 3 **SET ERR 0.5\n", " **********\n", " **********\n", " ** 4 **SET PRINT 1\n", " **********\n", " **********\n", " ** 5 **SET STR 1\n", " **********\n", " NOW USING STRATEGY 1: TRY TO BALANCE SPEED AGAINST RELIABILITY\n", " **********\n", " ** 6 **MIGRAD 1000 1\n", " **********\n", " FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.\n", " START MIGRAD MINIMIZATION. STRATEGY 1. CONVERGENCE WHEN EDM .LT. 1.00e-03\n", " FCN=25136.7 FROM MIGRAD STATUS=INITIATE 10 CALLS 11 TOTAL\n", " EDM= unknown STRATEGY= 1 NO ERROR MATRIX \n", " EXT PARAMETER CURRENT GUESS STEP FIRST \n", " NO. NAME VALUE ERROR SIZE DERIVATIVE \n", " 1 mean 1.00653e+00 2.00000e+00 2.02444e-01 3.46428e+01\n", " 2 sigma 3.01930e+00 9.90000e-01 2.22272e-01 2.14205e+01\n", " ERR DEF= 0.5\n", " MIGRAD MINIMIZATION HAS CONVERGED.\n", " MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.\n", " COVARIANCE MATRIX CALCULATED SUCCESSFULLY\n", " FCN=25136.7 FROM MIGRAD STATUS=CONVERGED 25 CALLS 26 TOTAL\n", " EDM=1.96964e-05 STRATEGY= 1 ERROR MATRIX ACCURATE \n", " EXT PARAMETER STEP FIRST \n", " NO. NAME VALUE ERROR SIZE DERIVATIVE \n", " 1 mean 1.00330e+00 3.04253e-02 3.34944e-04 1.02604e+00\n", " 2 sigma 3.01694e+00 2.22842e-02 5.41347e-04 6.16930e-01\n", " ERR DEF= 0.5\n", " EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 2 ERR DEF=0.5\n", " 9.257e-04 2.032e-05 \n", " 2.032e-05 4.966e-04 \n", " PARAMETER CORRELATION COEFFICIENTS \n", " NO. GLOBAL 1 2\n", " 1 0.02997 1.000 0.030\n", " 2 0.02997 0.030 1.000\n", " **********\n", " ** 7 **SET ERR 0.5\n", " **********\n", " **********\n", " ** 8 **SET PRINT 1\n", " **********\n", " **********\n", " ** 9 **HESSE 1000\n", " **********\n", " COVARIANCE MATRIX CALCULATED SUCCESSFULLY\n", " FCN=25136.7 FROM HESSE STATUS=OK 10 CALLS 36 TOTAL\n", " EDM=1.96993e-05 STRATEGY= 1 ERROR MATRIX ACCURATE \n", " EXT PARAMETER INTERNAL INTERNAL \n", " NO. NAME VALUE ERROR STEP SIZE VALUE \n", " 1 mean 1.00330e+00 3.04246e-02 6.69888e-05 1.00499e-01\n", " 2 sigma 3.01694e+00 2.22838e-02 1.08269e-04 -4.23243e-01\n", " ERR DEF= 0.5\n", " EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 2 ERR DEF=0.5\n", " 9.257e-04 1.984e-05 \n", " 1.984e-05 4.966e-04 \n", " PARAMETER CORRELATION COEFFICIENTS \n", " NO. GLOBAL 1 2\n", " 1 0.02926 1.000 0.029\n", " 2 0.02926 0.029 1.000\n", "[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization\n" ] } ], "source": [ "gauss.fitTo(data);" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The results are in agreement, because the results of a fit cannot depend on the minimizer. Technically, RooFit uses a different minimizer than iminuit by default. Unless you change some options, RooFit uses the original MINUIT Fortran implementation translated to C, while iminuit uses the rewritten Minuit2 C++ library." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Just doing the fitting with iminuit does not offer you a lot of advantages. Eventually, you want to switch completely. The equivalent of this exercise in pure Python is the following." ] }, { "cell_type": "code", "execution_count": 6, "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 = 4.998e+04 Nfcn = 29
EDM = 3e-06 (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", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 mu 0.96 0.03
1 sigma 2.985 0.022
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
mu sigma
mu 0.000906 0 (0.027)
sigma 0 (0.027) 0.000483
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 4.998e+04 │ Nfcn = 29 │\n", "│ EDM = 3e-06 (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 │ mu │ 0.96 │ 0.03 │ │ │ │ │ │\n", "│ 1 │ sigma │ 2.985 │ 0.022 │ │ │ │ │ │\n", "└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌───────┬───────────────────┐\n", "│ │ mu sigma │\n", "├───────┼───────────────────┤\n", "│ mu │ 0.000906 0 │\n", "│ sigma │ 0 0.000483 │\n", "└───────┴───────────────────┘" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABP0UlEQVR4nO3de3xT9f0/8NdpoIXSC5SWXkih6LRYBHRs1narE+ELAm6yUqfIlG1MhqtKwflF/DpB/A6c+hWqosjGwN+8UyMIQxERaLUFAbkWqICFNumd0vRGm/bk/P7AhIYm6Umb5Jykr+fjkceDnPNJ8j7G0hfncxMkSZJAREREpCIBShdAREREdDUGFCIiIlIdBhQiIiJSHQYUIiIiUh0GFCIiIlIdBhQiIiJSHQYUIiIiUh0GFCIiIlKdPkoX0B1msxllZWUIDQ2FIAhKl0NEREQySJKEhoYGxMXFISDA+T0SnwwoZWVliI+PV7oMIiIi6obS0lJotVqnbXwyoISGhgK4fIFhYWEKV0NERERy1NfXIz4+3vp73BmfDCiWbp2wsDAGFCIiIh8jZ3gGB8kSERGR6jCgEBERkeowoBAREZHqMKAQERGR6jCgEBERkeowoBAREZHqMKAQERGR6jCgEBERkeowoBAREZHqMKAQERGR6jCgEBERkeowoBAREZHqMKAQERGR6jCgEFGv0NTUBEEQIAgCmpqalC6HiLrAgEJERESqw4BCREREqsOAQkRERKrDgEJERESqw4BCREREqsOAQkRERKrDgEJERESqw4BCREREqsOAQkS9giiK1j/n5ubaPCci9WFAISK/p9PpkJSUZH0+depUJCQkQKfTKVgVETnDgEJEfk2n0yEjIwMGg8HmuMFgQEZGBkMKkUoxoBCR3xJFEfPnz4ckSZ3OWY5lZWWxu4dIhRhQiMhv5eXlQa/XOzwvSRJKS0uRl5fnxaqISA4GFCLyW+Xl5S61447HROrBgEJEfis2Ntat7YjIexhQiMhvpaWlQavVQhAEu+cFQUB8fDzS0tK8XBkRdYUBhYhUrSfdLhqNBtnZ2QDQKaRYnq9atQoajcY9xRKR2zCgEJFfS09PR05ODuLi4myOa7Va5OTkID09XaHKiMiZPkoXQETkbpIkYe/3tThUehFnq5pwtjoaw+a+AcOSXwEAfvP0avx2xq/wX6NcH3vS1NSEkJAQAEBjYyMGDBjg1tqJ6DIGFCLyGxcaW5FzUI/395eiuMa2O8hsMln/XNAUjX3vH8GAwOO444Zo3H/LMKRcO9jb5RKREwwoROTz2kQz3th9Fq/tOgNTu1n265pMIrYcKcOWI2WYlBSNhbcP92CVROQKBhQi8mmHSi7iyY+OoaiyoUfv8/mJSuw8es76PDc3F5MmTeIAWiKFcJAsEfkkSZLw8udFmPFGfo/DCQA0F+Xj/JvzrM+5oSCRshhQiMjntIlmLPzwCF758gzMnbfZcVlzUT6qNy2H2HjB5jg3FCRSDgMKEalax438cnNzcbGpBbP/9Q0+PmRw8ir5JLOI2p1r7Z+TJEjghoJESmBAISLV0ul0SEpKsj6fOnUqYocOwxfbtrjtM1r1hRAbahw34IaCRIpgQCEiVdLpdMjIyIDBYHunpNVYjepNy9FclG89Zja14Pzf78L5v98Fs6nFpc8RGy/Kapd3+DuX3peIeoYBhYhURxRFzJ8/H5LkeIBJ7c61kMzyu10CAvth+KKtGL5oKwIC+1mPa0IGyXr9vw7V4XBpXacuJ3b9EHkGAwoRqU5eXh70er3TNmJDDVr1hT3+rCDtKGhCI5220YRGwhw9EtOfeBmJI2+wHudMHyLPYUAhItUpLy+X1U5u94wzQoAGERPmOm0TMWEuLp3eh/MfPIeK8jKbc5zpQ+QZDChEpDqxsfL2yJHbPdOV4MRURE1/CpoQ2+XuNaGRiJr+FPpfl+x0pg/AmT5E7saAQkSqk5aWhti4oU7baEIjEaQd5bbPDE5MReyc163PozKWYui8dQhOTO1ypo/EmT5EbseAQkTqIwTgml9mOm0SMWEuhAD3LkPf8f36xd9ofS63K0lu1xQRdY0BhYgU09TUBEEQIAgCmpqu7D6c/cV30A8c47TbJTgx1Wt1yu1Kkts1RURdcymgLF261PqXieUxcuRI6/mWlhZkZmZi8ODBCAkJwYwZM1BZWWnzHiUlJZg2bRqCg4MxZMgQPPHEE2hvb3fP1RCRzzuqr8Nru84AcN7t0lHH6cYtpcddmn4sR1czfQRBQHx8PNLS0tz6uUS9mct3UEaNGoXy8nLr46uvvrKeW7BgAbZs2YKNGzdiz549KCsrQ3p6uvW8KIqYNm0aTCYT8vPz8dZbb2HDhg145pln3HM1ROTTzGYJT286brO/jqNuF4vmonyUr/uz9Xl1zlIY1syxWcitp+TM9Fm1ahV3PiZyI5cDSp8+fRATE2N9REZe/leF0WjEunXr8PLLL+OOO+7AuHHjsH79euTn52Pv3r0AgM8//xwnTpzA22+/jZtuuglTpkzBc889h9WrV8NkMrn3yojI57yz7zyO6o2y2zva5E9sqOm02mxPOZvpc+tDf8Ovf/1rt30WEXUjoJw+fRpxcXG45pprMGvWLJSUlAAADh48iLa2NkycONHaduTIkRg2bBgKCgoAAAUFBRg9ejSio6OtbSZPnoz6+noUFjpecKm1tRX19fU2DyLyL9UNrXhxe5Hs9s42+bNwdbXZrjjqciobNAbv7y912+cQkYsBJTk5GRs2bMBnn32GN954A8XFxUhLS0NDQwMqKioQGBiIgQMH2rwmOjoaFRUVAICKigqbcGI5bznnyIoVKxAeHm59xMfHu1I2EfmAl7YXob5F/ni0Ljf5g/tWm+3IUZfTim0nUdXg2j5AROSYSwFlypQpuOeeezBmzBhMnjwZ27ZtQ11dHT788ENP1QcAWLx4MYxGo/VRWsp/qRD5m0+OlHXdqAO5U3/dsdqsHPUt7Vi25YRXPouoN+jRNOOBAwfi+uuvx5kzZxATEwOTyYS6ujqbNpWVlYiJiQEAxMTEdJrVY3luaWNPUFAQwsLCbB5E1LvJnfrrymqzjjYUlGvr0XLsOlXl8uuIqLMeBZTGxkacPXsWsbGxGDduHPr27YudO3dazxcVFaGkpAQpKSkAgJSUFBw7dgxVVVd+gHfs2IGwsDAkJSX1pBQi6mXkbvLnztVm5Xh603E0m7h0AlFPuRRQ/vKXv2DPnj04d+4c8vPz8etf/xoajQYzZ85EeHg45syZg4ULF2LXrl04ePAgfv/73yMlJQW33norAGDSpElISkrCAw88gCNHjmD79u14+umnkZmZiaCgII9cIBGpV+mFRuufXV2/RO4mf+5ebbYrhrpLWLnjO69+JpE/6uNKY71ej5kzZ+LChQuIiorCz3/+c+zduxdRUVEAgJUrVyIgIAAzZsxAa2srJk+ejNdfvzLiXaPRYOvWrXj44YeRkpKCAQMGYPbs2Vi2bJl7r4qIVE+n0+GBP86zPq/OWQpNaCQiJsy1WYjN0u1ij2Xqb+0Xb9pMNbb3Pt70r6/P4e6bhuLGoeGKfD6RPxAky1acPqS+vh7h4eEwGo0cj0Lkg3Q6HWZkZAAO/vpxdSl7saUJ+ux7L782Yyn6j7jZY3dOzKYWlK7MAADEL8hxOFblZz8ajHf+eKtHaiDyVa78/nbpDgoRUU+JoojH5s93GE6Ay+uX9L8uWXbI6Gq1WXdydkeno6/PXMC+7y8g+ZrBXbYlos64WSAReYSjjQDz8vJg0OudvtYT65d4k9nUgvN/vwu3Xhtpc+1EJB8DChF5VaneIKudt9Yv8bT8M84XkyMi+xhQiMirTtbL635xZf0SNXv1y9NKl0DkkxhQiMhrTO1m7LoYocr1SzzlcKkRu4o6L97mqAuMiC5jQCEir/nwQCnKG9pUuX6JJ63iuihELmNAISKvMLWb8cbuswCurF+iCbGd4aIJjXR5irEvOKI3YseJyq4bEpEVpxkTkVdsPFgKQ90l6/PgxFQEDR/rlvVL5E79VdLKHd9h4g1DIAiC0qUQ+QTeQSEijxDFK8vW79q9B6t3du7m8Ob6Jd7Uccl+yxL+J8rr8dnxCgWrIvItDChE5LKuBnjqdDqbDUB/edc0fPP8/WguyvdmmYpoLspH+bo/W59X5yyFYc0cNBflY+UX38Fs9rnFu4kUwYBCRG6l0+mQkZEBg8F2vROxoQbVm5b7dUhpLspH9ablNvsCAVeu/XDu59h6rFyh6oh8CwMKEbmNKIqYP38+nG3xVbtzrUu7FvsKySyidudap21qd67Fys9PQpR5F4VTkak3Y0AhIrfJy8uD3s+XsXekVV8IscH5qrFiQw1OHvoGnxwx2IzRyc3NtXlORAwoRORG5eXyui/8ZRn7juRek9h4Ecte3WAzRmfq1KlISEiATqfzVHlEPocBhYjcJjY2VlY7f1nGviO519R20YDD6//aaYyOwWBARkYGQwrRDxhQiMht0tLSoNVqna710XEZe8v6JcMXbUVAYD9vlekRQdpRXS7hHxAyGI2Ht9s9Zxm3k5WVxe4eIjCgEFE3OBo/odFokJ2d7fS1/raMvYUQoOlyCf/QsXd2muHTkSRJKC0tRV5enrvLI/I5DChE5JKr1zi5evxEeno6/u/Nt3rNMvYddbWEf9+IOFnvI3csD5E/Y0AhItkcrXFy9fiJ0rDRiJ3zuvV8VMZSDJ23zq/DiUVwYqrDa5c7TkXuWB4if8aAQkSyOFvjpOP4iYq6Zmw9Wu63y9jL4ejauxqnIggC4uPjkZaWBsBxVxpRb8CAQkSydLXGiWX8xP+u+wgm0ezFynyHs3EqloHFq1atgkaj6bIrjcjfMaAQkSxyx0Vs23fSw5X4NkfjVLRaLXJycpCeni67K43InzGgEJEscsdFNPcJ9XAlvu/qcSqpmS+huLgY6enpsrvS2N1D/o4BhYhk6WqNE0EQ0G/gEOsaJ+Rcx3EqpYEJOFbWAEB+VxqnIpO/Y0AhIlk6rnFydUgRBAESgNDb/9irBsO607++KgYgvyuNU5HJ3zGgEJFs6enpyMnJQVyc7XoeWq0Wd85/oVdMI/aUT4+Xo8LYIrsrjVORyd8JkrN90VWqvr4e4eHhMBqNCAsLU7ocol7H8jMIANu2bcPo5DTc/lIu2s0+99eJYsymFpSuzAAAxC/IQUBgPzx8+7X4y39dh4SEBBgMBrvjUARBgFarRXFxMTQa3q0i3+LK72/eQSEil3X8xXjbbbfh3W/0DCdu8N43JTCJcNqVBlyZikzkzxhQiKhHWtpEvL+/VOky/EJdcxt0h/ROu9IsU5GJ/B0DChH1yNajZahtMildht94Z28JgMvjfU6cOGE9vm3bNutUZKLeoI/SBRCRb/t3wXmlS/BJAYH9MHzR1k7HT5TX40hpHcbGD+zUlcZuHepNeAeFiHrku8pGpUvwO+/vL1G6BCLFMaAQEanMJ4fL0NTarnQZRIpiQCEiUpkmk4hPjpQpXQaRohhQiMhlAwYMwL8LzmH4oq0ICOyndDl+6f1v2M1DvRsDChF1ywecWuxRR/RGnCyvV7oMIsVwFg8RuaywzIhjBqPSZfi9Twov2F1Nlqg34B0UInLZ+9/w7ok3fHzIgJY2UekyiBTBgEJELmlpE7H5sEHpMnqFhpZ2bD3KXYupd2JAISKXbDtWjvoWToH1Fg6Wpd6KAYWIXMJ9d7zrwPmLOF3ZoHQZRF7HgEJENpqamiAIAgRBQFNTk82576sb8U1xrUKV9V7vccwP9UIMKEQkG6cWK+PjQ3q0tnOwLPUuDChEJEubaMZH33JwrBIuNrfhs+MVSpdB5FUMKEQky86TlahpbFW6jF7rPQ6WpV6GAYWIZOHgWGXt/b4WxTVNXTck8hMMKETUpbK6S8j9rlrpMnq99/d37y6Ks4HPRGrFgEJEXfrwQCnMXHFdcZsOGSDyi6BeggGFiGyI4pXZIrm5uWhra8fGA3oFKyKLyvpWfHWmRukyiLyiRwHl+eefhyAIyMrKsh5raWlBZmYmBg8ejJCQEMyYMQOVlZU2ryspKcG0adMQHByMIUOG4IknnkB7O1emJFKaTqdDUlKS9fnUqVMxdNhwnN63U8GqqCPdtwyL1Dt0O6Ds378fb775JsaMGWNzfMGCBdiyZQs2btyIPXv2oKysDOnp6dbzoihi2rRpMJlMyM/Px1tvvYUNGzbgmWee6f5VEFGP6XQ6ZGRkwGCwnUpcXVGO6k3L0VyUr1Bl1NH2wgo0tLQpXQaRx3UroDQ2NmLWrFn4xz/+gUGDBlmPG41GrFu3Di+//DLuuOMOjBs3DuvXr0d+fj727t0LAPj8889x4sQJvP3227jpppswZcoUPPfcc1i9ejVMJpN7roqIXCKKIubPnw9Jsje+4fKx2p1rIZm5WJjSWtrM2HaMGwiS/+tWQMnMzMS0adMwceJEm+MHDx5EW1ubzfGRI0di2LBhKCgoAAAUFBRg9OjRiI6OtraZPHky6uvrUVhYaPfzWltbUV9fb/MgIvfJy8uDXu+860BsqEGr3v7PKHmXqwvmXT2uqONzIrVyOaC8//77+Pbbb7FixYpO5yoqKhAYGIiBAwfaHI+OjkZFRYW1TcdwYjlvOWfPihUrEB4ebn3Ex8e7WjYROVFeLu9f5GLjRQ9XQnLsP1eL0tpmWW3tjStKSEiATqfzVHlEbuFSQCktLcX8+fPxzjvvoF+/fp6qqZPFixfDaDRaH6WlXDCKyJ1iY2NltdOEDOq6EXmcJAEfyRgs62hckcFgQEZGBkMKqZpLAeXgwYOoqqrCj3/8Y/Tp0wd9+vTBnj178Morr6BPnz6Ijo6GyWRCXV2dzesqKysRExMDAIiJiek0q8fy3NLmakFBQQgLC7N5EJH7pKWlQavVQhAEh200oZEI0o7yYlXkzMeHnHfzOBtXZDmWlZXF7h5SLZcCyoQJE3Ds2DEcPnzY+vjJT36CWbNmWf/ct29f7Nx5ZUpiUVERSkpKkJKSAgBISUnBsWPHUFVVZW2zY8cOhIWF2dyGJCLv0Wg0yM7OBgCHISViwlwIARpvlkVOnL/QjP3nah2e72pckSRJKC0tRV5enifKI+qxPq40Dg0NxY033mhzbMCAARg8eLD1+Jw5c7Bw4UJEREQgLCwMjz76KFJSUnDrrbcCACZNmoSkpCQ88MADeOGFF1BRUYGnn34amZmZCAoKctNlEZGr0tPTkZOTg8cee8ymS0ATGomICXMRnJiqYHVkj+5bPX6aEGH3nNxxRXLbEXmb21eSXblyJe666y7MmDEDt912G2JiYmz6OTUaDbZu3QqNRoOUlBT89re/xYMPPohly5a5uxQiclF6ejpOnDhhfR6VsRRD561jOFGprUfL0dJmv4tG7rgiue2IvE2Q7C98oGr19fUIDw+H0WjkeBQiN2tqakJISAgAIH5BDgICvTcgnlz3ysyb8auxcZ2Oi6KIhIQEGAwGu+NQBEGAVqtFcXExNBp23ZF3uPL7m3vxEJGNQyWOxzWQ+jha+t7ZuCLL81WrVjGckGoxoBCRjU+OlCldArkg73QNqhpa7J6zjCuKi7O9w6LVapGTk2OzDQmR2jCgEJFVm2jGp8cru25IqiGaJWw+5DhUXj2uaNu2bSguLmY4IdVjQCEiq12nqmBs5kZ0vqarRds6duPcdttt7NYhn8CAQkRWmw67tscLqcOpigYcNxiVLoPIrVxaB4WI/Fd9Sxu+OFmFgMB+GL5oq9LlkIt03xpw49BwpcsgchveQSEiAMC2o+UwtZuVLoO6acvRMohmn1s1gsgh3kEhIgBd7+1C6mQ2taB0ZQbOA9h5VzEmjU1QuiQit+AdFCKCoe4SvnGyrwv5hq1HKpQugchteAeFiLDpkAG+t6Y0Xe2Lk5VoaRPRr6/tLJ0BAwbYXU2WSM14B4WIsIndO36hsbUdX5zkOjbkHxhQiHq54wYjTlc1Kl0GuckmJ4u2EfkSBhSiXo53T/zLnu+qUNdsUroMoh5jQCHqxUSzxL13fJxkFq1/bik9DlNbO/5zrFzBiojcg4NkiXqRpqYmhISEAAAaGxtx0NCMqoZWhaui7mouykftF29an1fnLIUmNBKrjQsxK3mxgpUR9RzvoBD1Yuze8V3NRfmo3rQcYuMFm+NiQw0K1j6Ff/z7PYUqI3IPBhSiXuqSScT2Qq6b4Ysks4janWudtln0+OMQRdFpGyI1Y0Ah6qW+OFmJJhN/gfmiVn0hxIYap20uVpcjLy/PSxURuR8DClEvxcGxvktsvCir3YET33u4EiLPYUAh6qUKzl7ouhGpkiZkkKx2RQ2arhsRqRQDClEv0nFMQtP5YzZTVMl3BGlHQRMa6bSNJjQSx9tjucQ9+SwGFKJeQqfTISkpyfq8OmcpDGvmoLkoX8GqqDuEAA0iJsx12iZiwlyUN7RhXzE3gSTfxIBC1AvodDpkZGTAYLCdViw21KB603KGFB8UnJiKqOlPQRMy2Oa4JjQSUdOfQnBiKgBg82FOJSffxIBC5OdEUcT8+fOd3uqv3bmW3T0+KDgxFbFzXrc+j8pYiqHz1lnDCQBsO1YBU7tZifKIeoQBhcjP5eXlQa/XO20jNtSgVV/opYrInYSAKwNh+8XfaPMcAIyX2rCrqMrbZRH1GAMKkZ8rL5e3L4vcqavke9jNQ76IAYXIz8XGxspqJ3fqKvmenSer0NjarnQZRC5hQCHyc2lpadBqtRAEwWEbTWgkgrSjvFgVeVNruxmfHee2BuRbGFCI/JxGo0F2djacrYYRMWFup7EL5F/YzUO+hgGFqBdIT0/H5Mf+3uWUVPJf+WcvoKaxVekyiGTro3QBROR5ja3tKA65EbFzXoc++14Al6ek9h9xM++c+LiAwH4Yvmhrl+1Es4StR8rwu5+N8EJVRD3HOyhEvcBnxyvQ2m7uckoq+bfN3CCSfAgDClEvwPEHBACHSupQWtusdBlEsjCgEPm5msZW5HPnYvoBwyr5CgYUIj+39UgZRDN3tKXLNh9mNw/5BgYUIj/HcQfU0emqRpwoq1e6DKIucRYPkR8rudCMQyV11udyZ3yQf9t8xICkuDClyyByindQiPwYxxuQPVuPlDvd3ZpIDRhQiPwYu3fIHkPdJew/x80hSd0YUIj8VGGZEWeqGpUug1SKd9dI7RhQiPxAU1MTBEGAIAhoamoCAHzC2RrkxLZj5WgTzUqXQeQQAwqRH5IkCVvYvUNOXGxuQ97paqXLIHKIAYXID31TXIsyY4vSZZDKcU0UUjMGFCI/xMGxJMeOE5W4ZBKVLoPILgYUIj/TJprx6bFypcsgH9BsEvH5iQqlyyCyiwGFyA+I4pV/Bb/x3hbUNrJ7h+SxDKa2N9CaSEkMKEQ+TqfTISkpyfr8qXmzYFgzB81F+QpWRb4i93Q16ppNSpdB1AkDCpEP0+l0yMjIgMFgu6aF2FCD6k3LGVKoS22ihP+wS5BUyKWA8sYbb2DMmDEICwtDWFgYUlJS8Omnn1rPt7S0IDMzE4MHD0ZISAhmzJiByspKm/coKSnBtGnTEBwcjCFDhuCJJ55Ae3u7e66GqBcRRRHz5893umR57c61kMwcBEnOcTYPqZFLAUWr1eL555/HwYMHceDAAdxxxx24++67UVhYCABYsGABtmzZgo0bN2LPnj0oKytDenq69fWiKGLatGkwmUzIz8/HW2+9hQ0bNuCZZ55x71UR9QJ5eXnQ6/VO24gNNWjVF3qpIvJV+8/VosJ4SekyiGy4FFB++ctfYurUqbjuuutw/fXX429/+xtCQkKwd+9eGI1GrFu3Di+//DLuuOMOjBs3DuvXr0d+fj727t0LAPj8889x4sQJvP3227jpppswZcoUPPfcc1i9ejVMJvaBErmivFzebXmxkXuukHOSBGw5dCXs5ubm2gy8JlJCt8egiKKI999/H01NTUhJScHBgwfR1taGiRMnWtuMHDkSw4YNQ0FBAQCgoKAAo0ePRnR0tLXN5MmTUV9fb70LY09rayvq6+ttHkS9XWxsrKx2mpBBHq6EfF1zUT6enPVf1udTp05FQkICdDqdglVRb+dyQDl27BhCQkIQFBSEefPm4eOPP0ZSUhIqKioQGBiIgQMH2rSPjo5GRcXlefYVFRU24cRy3nLOkRUrViA8PNz6iI+Pd7VsIr+TlpYGrVYLQRActtGERiJIO8qLVZGvaS7KR/Wm5Wg11tgcNxgMyMjIYEghxbgcUBITE3H48GHs27cPDz/8MGbPno0TJ054ojarxYsXw2g0Wh+lpaUe/TwiX6DRaJCdnQ0ADkNKxIS5EAI03iyLfIhkFlG7c639cz8Mvs7KymJ3DynC5YASGBiIH/3oRxg3bhxWrFiBsWPHIjs7GzExMTCZTKirq7NpX1lZiZiYGABATExMp1k9lueWNvYEBQVZZw5ZHkQEpKenIycnB2GDh9gc14RGImr6UwhOTFWoMvIFrfpCiA01Ds9LkoTS0lLk5eV5sSqiy3q8DorZbEZrayvGjRuHvn37YufOndZzRUVFKCkpQUpKCgAgJSUFx44dQ1VVlbXNjh07EBYWZrPQFBHJl56ejnEL/2V9HpWxFEPnrWM4oS7JHUAtd0A2kTv1caXx4sWLMWXKFAwbNgwNDQ149913sXv3bmzfvh3h4eGYM2cOFi5ciIiICISFheHRRx9FSkoKbr31VgDApEmTkJSUhAceeAAvvPACKioq8PTTTyMzMxNBQUEeuUAif1dU0YDT1VemiPaLv5HdOiSL3AHUcgdkE7mTSwGlqqoKDz74IMrLyxEeHo4xY8Zg+/bt+K//ujz6e+XKlQgICMCMGTPQ2tqKyZMn4/XXX7e+XqPRYOvWrXj44YeRkpKCAQMGYPbs2Vi2bJl7r4qoF9l82NB1IyI7grSjoAmNdNjNIwgCtFot0tLSvFwZESBIzpahVKn6+nqEh4fDaDRyPAr1ej//+5coqbyI0pUZAID4BTkICOyncFXkKyyzeK5mGXidk5Njs+AmUU+48vube/EQ+bCD52uhv8gVQKn7ghNTETX9KWhCBtsc12q1DCekKAYUIpVramqCIAgQBAFNTU0257iHCrlDcGIqYudc6Y4fPvM5nPzuDMMJKcqlMShEpB7tohnbftiFNiCwH4Yv2qpwReTLOg6sNsfcgNzTFzBlNAfHknJ4B4XIR311pgY1jdzDijyDd+dIaQwoRD7qE/4CIQ/aVVSF+pY2pcugXowBhcgHtbSJ+PxEZdcNibqptd2Mz4473iONyNMYUIhUruM+KLm5uRBFETtOVKKxtV3BqsjfWMYxDV+01TpN/eq7dM4GbBO5GwMKkYrpdDqbbSCmTp2KhIQEZP/zbQWrot4i/2wNyo2cxk7KYEAhUimdToeMjAwYDLYrxRoMBux8bRGai/IVqox6C7MEfHyIKxWTMhhQiFRIFEXMnz8f9hZ6thyr3bkWklnsdJ7InT46qFe6BOqlGFCIVCgvLw96vfNfDGJDDVr1hV6qiHqrs9VNOFRyeddje+OhiDyFAYVIheRuby82XvRwJUTAR9/qHY6H0ul0ClZG/owBhUiF5G5vrwkZ5OFKiIC339vocDxURkYGQwp5BAMKkQqlpaVBq9Vad5S1RxMaiSDtKC9WRb2RZBah//QNp+OhsrKy2N1DbseAQqRCGo0G2dnZAOAwpERMmGuzfwqRJ7TqCyE21Dg8L0kSSktLkZeX58WqqDdgQCFSqfT0dOTk5CAuLs7muCY0ElHTn0JwYqpClVFvIneck9xxU0RycTdjIhVLT0/HxIkTER4eDgCIyliK/iNu5p0T8hq545zkjpsikot3UIhUTqO5Ekb6xd/IcEJeFaQdBU1opMPzgiAgPj4eaWlpXqyKegMGFCIickgI0CBiwlz7534YH7Vq1SqbIE3kDgwoRCpX12xSugTq5YITUxE1/SloQgbbHNdqtcjJyUF6erpClZE/4xgUIpX79BgHH5LyghNTETR8LPTZ9wIAtm3bhkmTJvHOCXkMAwqRyv3n1EUMX7RV6TKIbMY/3XbbbQwn5FHs4iFSseKaJhwqqVO6DCIir2NAIVKxj7/lTrKkThwbRZ7GgEKkUpIk4ePDhq4bEimAY6PI0zgGhUil9p+7iNLaS0qXQWQVENjPOh5qy8mL+OMdChdEfo13UIhU6uND7N4h9TpSWoczVY1Kl0F+jAGFSIVa2kT85yhvoZO6fcQxUuRBDChEKrTzZBXqW9qVLoPIqY+/NcBslpQug/wUAwqRCn1woFTpEoi6VFHfgq/P1ihdBvkpBhQilSk3XsJXp6uVLoNIlo8OspuHPIMBhUhlcg7owbvm5Cu2F1aisZXdkeR+DChEKiJJEnI48JB8yKU2Ef85WqZ0GeSHGFCIFNTU1ARBECAIApqamrD3+1qcv9CsdFlELvnoIBcUJPdjQCFSkY0cHEs+aP/5WpQwWJObMaAQqURDSxs+PV6hdBlELpMkrolC7seAQqQgURStf37535+guZUbsJFv0h3SQ5I4upvchwGFSCE6nQ5JSUnW5y8umA3DmjloLspXsCqi7imtvYRvimuVLoP8CAMKkQJ0Oh0yMjJgMNgOLhQbalC9aTlDCvkkdvOQOzGgEHmZKIqYP3++09vhtTvXQjKLDs8TqdG2YxW4ZOL/t+QeDChEXpaXlwe93vm/NMWGGrTqC71UEZF7NLa2Y3shB3qTezCgEHlZebm8XYrFxoseroTI/T4oOG2ztg9RdzGgEHlZbGysrHaakEEeroTI/QrOXlC6BPITDChEXpaWlgatVgtBEBy20YRGIkg7yotVEbkH95Eid2FAIfIyjUaD7OzsH57ZDykRE+ZCCNB4rygiN+k4uDs3N9dmrR8iVzCgECkgPT0dOTk5CAwbbHNcExqJqOlPITgxVaHKiLqvuSgf5ev+bH0+depUJCQkQKfTKVgV+SoGFCKFDB93O4b8frX1eVTGUgydt47hhHxSc1E+qjcth9hoOwbFYDAgIyODIYVcxoBC5CFX71R8tXf2ldh04/SLv5HdOuSTJLOI2p1r7Z/7Yb2frKwsdveQS1wKKCtWrMBPf/pThIaGYsiQIZg+fTqKiops2rS0tCAzMxODBw9GSEgIZsyYgcrKSps2JSUlmDZtGoKDgzFkyBA88cQTaG9v7/nVEPmIumYT/nNU3nRjIrVr1RdCbKhxeF6SJJSWliIvL8+LVZGvcymg7NmzB5mZmdi7dy927NiBtrY2TJo0yeZfhwsWLMCWLVuwceNG7NmzB2VlZUhPT7eeF0UR06ZNg8lkQn5+Pt566y1s2LABzzzzjPuuikjlNh7Qo7XdrHQZRG4hd80euWsAEQGAIPVg+8nq6moMGTIEe/bswW233Qaj0YioqCi8++67yMjIAACcOnUKN9xwAwoKCnDrrbfi008/xV133YWysjJER0cDANasWYNFixahuroagYGBXX5ufX09wsPDYTQaERYW1t3yiTyqqakJISEhAIDGxkYMGDAAwOV/Td7xf3tQXMNFrMg/tJQcReV7T3XZbteuXbj99ts9XxCpliu/v3s0BsVoNAIAIiIiAAAHDx5EW1sbJk6caG0zcuRIDBs2DAUFBQCAgoICjB492hpOAGDy5Mmor69HYSGX9ib/0bG/veN0y6/PXGA4Ib8SpB0FTWikw/OCICA+Ph5paWlerIp8XbcDitlsRlZWFn72s5/hxhtvBABUVFQgMDAQAwcOtGkbHR2NiooKa5uO4cRy3nLOntbWVtTX19s8iNRMp9MhKSnJ+rzjdMu3955XsDIi9xMCNIiYMNdpm1WrVkGj4SBwkq/bASUzMxPHjx/H+++/78567FqxYgXCw8Otj/j4eI9/JlF36XQ6ZGRkwGAw2By3TLfctInTLcn/BCemImr6U9CEdF7b53d/fcVmLCKRHN0KKI888gi2bt2KXbt2QavVWo/HxMTAZDKhrq7Opn1lZSViYmKsba6e1WN5bmlztcWLF8NoNFofpaWl3SmbyONEUcT8+fNhb2iXJEmQJKBmx1qb1TaJ/EVwYipi57xufW5Z2+dkvxsgcg18cpFLAUWSJDzyyCP4+OOP8eWXX2LEiBE258eNG4e+ffti586d1mNFRUUoKSlBSkoKACAlJQXHjh1DVVWVtc2OHTsQFhZmc0u8o6CgIISFhdk8iNQoLy8Per3eSQsJYkMNWvUcb0X+yd7aPuXGFuw8WenkVUSd9XGlcWZmJt59911s3rwZoaGh1jEj4eHh6N+/P8LDwzFnzhwsXLgQERERCAsLw6OPPoqUlBTceuutAIBJkyYhKSkJDzzwAF544QVUVFTg6aefRmZmJoKCgtx/hUReJHcapdxpmUT+4t1vSjBplP275ET2uBRQ3njjDQDoNE1s/fr1+N3vfgcAWLlyJQICAjBjxgy0trZi8uTJeP31K7f8NBoNtm7diocffhgpKSkYMGAAZs+ejWXLlvXsSohUIDY2VlY7TcggD1dCpC6531WjtLYZ8RHBSpdCPqJH66AoheugkFqJooiEhAQYDAa741CAy4MGh85bx2Xtqdf58+3X4r/vHKl0GaQgr62DQkS2NBoNsrOzAVxe+8GeiAlzGU6oV/rwgB5tIldQJnkYUIjcLD09HTk5OYiLi7M5rgmNRNT0p7hbMfVaNY2t2F5of70roqsxoBB5QHp6Ok6cOGF9bpluyXBCvd07e0uULoF8BAMKkYd0XDXTMt2SqLcr+P4CDn1fAUEQIAiCzWazRB0xoBARkVf9v/xzSpdAPsClacZEJN+AAQOwKOcI3t/PlY+JOtp6VN56QdS78Q4KkYfUNLbi40OGrhsS9TItpjbrnzvu9E3UEQMKkYf8u+A8Wts5pZKoo+aifJSv+7P1ecedvok6YkAh8oCWNhHv7DuvdBlEqtJclI/qTcshNl6wOW7Z6ZshhTpiQCHygE2HDKhpNCldBpFqSGYRtTvX2j/3w6rLWVlZ7O4hKwYUIg9Y91Wx0iUQqUqrvhBiQ43D85IkobS0FHl5eV6sitSMAYXIzfZ8V43TVY1Kl0GkKnJ38Ja7Izj5PwYUIjf7Z973SpdApDpyd/CWuyM4+T8GFCI3+q6yAXmnHd/GJuqtgrSjoAmNdHheEATEx8cjLS3Ni1WRmjGgELkR754Q2ScEaBAxYa7TNqtWrbLZIoJ6NwYUom5oamrqtJdITWMrNh0uU7gyIvUKTkxF1PSnoAkZbHNcExqJectWIz09XaHKSI241D2Rm/y/gvMwcWE2IqeCE1MRNHws9Nn3Ari803f/ETfjaN8QiGYJmgBB4QpJLXgHhcgNWtpEvLOXC7MRydFxZ2/LTt/6i5ew40SF9bi9u5TUuzCgEHVDx8WkcnNz8dGBElxo4sJsRD3B9YOoIwYUIhfpdDokJSVZn0+dOhVzpiSjuShfwaqIfEdAYD8MX7QVwxdtRUBgP+vx/ecu4sC5WgUrIzVhQCFygU6nQ0ZGBgwG212KW43VqN60nCGFqIde/fIMgM53KbkEfu/DgEIkkyiKmD9/vnXfEHtqd66FZOZfpETdtee7aqz8x7873aXkjse9DwMKkUx5eXnQ6/VO24gNNWjVF3qpIiL/01yUj4VzH+x0l5I7Hvc+DChEMsndI0TuniNEZIs7HlNHDChEMsndI0TuniNEZIs7HlNHDChEMqWlpUGr1UIQHC8kpQmNRJB2lBerIvIf3PGYOmJAIZJJo9EgOzv78hMHISViwlybRaiISD7ueEwdMaAQuSA9PR05OTkIDrfdlVUTGomo6U8hODFVocqIfB93PKaOGFCIXHR98gREzH7N+jwqYymGzlvHcELUQ852PLZ0rXLH496DAYXIRa/sPG13LxEi6jlHOx5rtVrk5ORwx+NehLsZE7ngVEU9dpystC7VTUTud/WOx79dugYbnv4j75z0MryDQuSCV3eegZOFZInITTrelfymORrGFq590tvwDgqRTGeqGvDpcU5vJPKGjncpWwH8M+97/PedI5UtiryKd1CIZHr1yzMw8+4JkSL+XXAeF5tMSpdBXsSAQiTD99WN2HqUd0+IlNLQ2o7Xdp1RugzyIgYUIhle2XkaIm+fECnq33vPQ3+xWekyyEsYUIi6cExvxOYjZUqXQdTrmdrN+L/Pv1O6DPISBhSiLizfdpIzd4hUYvNhA06U1StdBnkBAwqRE7tOVaHg+wtKl0FEPzBLwPOfnVK6DPICBhQiB0SzhBWfnlS6DCK6Su531cg/U6N0GeRhDChEDmw8UIrvKhuVLoOI7Hj+s1NobGyEIAgQBAFNTU1Kl0RuxoBCZEezqR0v7+BgPCK1Oqo34j9HDNbnubm5EEWuNutPGFCI7PhHbjGqGlqVLoOIHGguyseDU39ufT516lQkJCRAp9MpWBW5EwMK0VXOVdQi67+ux/m/3wWzqUXpcojoKs1F+ajetBymettxKAaDARkZGQwpfoIBhegqr+68MkOgpfQ4JDNvGxOphWQWUbtzrf1zP6wHkJWVxe4eP8CAQtTB6vXvIHver6zPq3OWwrBmDpqL8hWsiogsWvWFEBscz+CRJAmlpaXIy8vzYlXkCQwoRD/Q6XR45A+/hdhou+6J2FCD6k3LGVKIVEBsvCirXXk5987ydQwo1Ks0NTXZnZYoiiIeznzU6Wtrd65ldw+RwjQhg2S1i42N9XAl5GkMKNSrdOyX7jgtcffuPaiqcL7fjthQg1Z9oUfrIyLngrSjoAmNdHheEATEx8cjLS3Ni1WRJ7gcUHJzc/HLX/4ScXFxEAQBmzZtsjkvSRKeeeYZxMbGon///pg4cSJOnz5t06a2thazZs1CWFgYBg4ciDlz5qCxkQtikWfpdDokJSVZn3eclvjenqOy3kPu7WUi8gwhQIOICXOdtlm1ahU0Go2XKiJPcTmgNDU1YezYsVi9erXd8y+88AJeeeUVrFmzBvv27cOAAQMwefJktLRcma45a9YsFBYWYseOHdi6dStyc3Mxd67z/+GIekKn0yEjIwMGg8HmuGVa4uY9B2W9j9zby0TkOcGJqYia/hQ0IYNtjmtCI/Hg09lIT09XqDJyJ0GSur9PqyAI+PjjjzF9+nQAl++exMXF4fHHH8df/vIXAIDRaER0dDQ2bNiA++67DydPnkRSUhL279+Pn/zkJwCAzz77DFOnToVer0dcXFyXn1tfX4/w8HAYjUaEhYV1t3zqJURRREJCAvR6vcM2ASGDIQCdBsh2pAmNxNB56yAE8F9mRGogtjRBn30vACAqYyn6j7gZQX374rOsNFwTFaJwdWSPK7+/3ToGpbi4GBUVFZg4caL1WHh4OJKTk1FQUAAAKCgowMCBA63hBAAmTpyIgIAA7Nu3z+77tra2or6+3uZBJFdeXp7TcAIA5sYLCLlpstM2ERPmMpwQqUjHn8d+8TdCCNDAJJrx183HFayK3MWtAaWiogIAEB0dbXM8Ojraeq6iogJDhgyxOd+nTx9ERERY21xtxYoVCA8Ptz7i4+PdWTb5ObnTDfsOGurwtnHU9KcQnJjqifKIqJsCAvth+KKtGL5oKwIC+1mPf33mAjYfvtyd62jmHqmfT8ziWbx4MYxGo/VRWlqqdEnkQ+RON9SEDEJwYipi57xuPRaVsRRD561jOCHyMf/7n5Oob2lzOHOP1M+tASUmJgYAUFlZaXO8srLSei4mJgZVVVU259vb21FbW2ttc7WgoCCEhYXZPIjkSktLg1arhSAIDttoQiMRpB0FwP5tYyLyLdUNrfjj0tUOZ+6R+rk1oIwYMQIxMTHYuXOn9Vh9fT327duHlJQUAEBKSgrq6upw8OCVWRNffvklzGYzkpOT3VkOEQBAo9EgOzsbAByGlI7jSxzdNiYi39FclI+Nf1/gcOYeQ4r6uRxQGhsbcfjwYRw+fBjA5YGxhw8fRklJCQRBQFZWFv73f/8Xn3zyCY4dO4YHH3wQcXFx1pk+N9xwA+6880489NBD+Oabb/D111/jkUcewX333SdrBg9Rd6SnpyMnJ6fT/2McX0Lkf7ihoH9weZrx7t27MX78+E7HZ8+ejQ0bNkCSJCxZsgRr165FXV0dfv7zn+P111/H9ddfb21bW1uLRx55BFu2bEFAQABmzJiBV155BSEh8qaFcZoxddfJ8xVISrg8JsUyLZFdOET+paXkKCrfe6rLdrt27cLtt9/u+YLIypXf3z1aB0UpDCjUHZIkYfb6/cj9rlrpUojIg5pO7EHNlhe7bPfuu+9i5syZXqiILBRbB4VIzdZ9VcxwQtQLcENB/8CAQr1CYZkRL3xWpHQZROQF3FDQPzCgkN+7ZBLx2HuHYBLNSpdCRF7gdEPBH2bycUNB9WNAIb+3bGshzlZzBUmi3sTRhoJ9QyOx7t/vcUNBH8CAQn7t02PleO8brjxM1BvZWxk69k//xB7TCAWrIrn6KF0AkaeUGy/hSd0xpcsgIgVp+g3A8EVbbY7tLqrG/ys4hwdTEpQpimThHRTyS2azhKz3D8N4qU3pUohIhZZvO4kzVQ0ALk99tWwo+Omnn3IBN5VgQCG/9PruM9hXXKt0GUSkUi1tZjz23mF8uDGH+/WoFAMK+Z1vSy5i1RenlS6DiFTuwO7PcO9vfsP9elSKAYX8SnVDKx599xDazT63QDIRedGV/Xo6/13B/XrUgQGF/EbVhYsYEtYP+YsnoPnsAUhm/sVCRPa16gshNtQ4PC9JEkpLS5GXl+fFqqgjBhTyCzqdDtdcl2h9Xp2zFIY1c9BclK9gVUSkVmLjRVntysvLPVwJOcKAQj5Pp9NhxowMNF203WdHbKhB9ablDClE1An361E/BhTyaaIoYu6fH4G9fmSL2p1r2d1DRDa4X4/6MaCQT/tXzjZcqHR+C1ZsqEGrvtBLFRGRL3C6Xw8679cjiiJ2796N9957D7t37+bgWS9gQCGfZai7hBUfFchqK7e/mYh6D0f79WhCB2PB82us+/XodDoMHz4c48ePx/3334/x48dzrRQv4FL35JOaWtsxZ8N+NPcJldVebn8zEfUuwYmp6H9d8uVZPY0XoQkZhCDtKGwx9sWMMzWoOJqLjIwM69RjC8taKTk5Odx40EMYUMjntIlmPPLutzhV0WDtR3Y2XVATGokg7SgvVkhEvkQI0KDfsDE2x0yiGXM27EPtukc7hRPg8jRkQRCQlZWFu+++29oVRO7DLh7yKW2iGQ+//S12FV2eseO8H/myiAlzIQTwLw8ics3F74+iqqLM4XmuleJZDCjkMyzh5IuTlTbHHfcjRyJq+lMITkz1ZplE5Ce4Voqy2MVDPqHF1IYZS/6FbwrPWvuIO94VcdSPzDsnRNRdXCtFWQwopHof5uTg93P/jOYOC7FpQiMRMWGuzd0Re/3IRETd1dUYN0EQoNVquVaKh7CLh1Ttw5wc3HvPPTbhBOAqsUTkec7GuAlC57VSyL0YUEi1Wkxt+MPcTKdtuEosEXmSozFuMbFxnGLsYeziIVVqaRMxY8m/0HSxymk7yyqx7NohIk+xN8Yt8rqbEDv2FqVL82u8g0KqU1nfgt+8WYBvCs/Kas9VYonI0yxj3AYk/QL9ho1BfasZD/5rH3IO6pUuzW8xoJCqHNMb8avXvsJRvVH2CHquEktESmgTJfxl4xG8tL3I7mJu1DPs4iGPEkUReXl5KC8vR2xsLNLS0joNKLO02bbvBN4/3gAh9gYIARquEktEPuG1XWdwvrYZL2aMQb++HDDrLgwo5DE6nQ6PPfYYDAaD9ZhWq0V2drbNJlzz58+HXn/lNmnHKcQRE+aietNyh5/BVWKJSA22HClDWd0lrH1gHAaHBCldjl8QJB+8L1VfX4/w8HAYjUaEhYUpXQ7ZodPp7G6wZZmal5OTAwB221hYVoFtLspH7c61NndS7K2DQkSktGERwfjn7J/g+mh5G5n2Nq78/mZAIbvkdM04e21CQoLNXZGOBEHA0KFD0S6aUVHueJ8LTWgkhs5bByFAA8kscpVYIlI1y99TwqU6PHjHTfjbw/egTx92VHTEgEI9Yq/b5equGWd2796N8ePHu6WW6JnLOYWYiFTP3p3e4EFD8Nqr2fj9rPtceq+e/ANR7Vz5/c1ZPGTD0jVz9d0Pg8GAjIwM6HS6Lt/DnRtncQoxEaldc1E+qjct7zSgv/liFf7w25n466p1st9Lp9MhISEB48ePx/3334/x48cjISFB1t+9/oYBhaxEUcT8+fPtjgmxHMvKyoIoOl+51Z0bZ3EKMRGpmWQWUbtzrdM2K555Ev/94SE0m9ohiiJ2796N9957D7t377b5+9Qd/0D0J+ziISu5XTO7du3C7bff7vA2pGUMisFgcDgANiBkMAQAYuMFh5/TcQwKEZEatZQcReV7T3XZLnrmcgzuY0LNjrWoqrgy9s7SfX733Xd3OXZPq9WiuLjYp7t72MVD3SK3a6a8vNzpbUiNRoPs7OwfWgt232PwxD8hYuKfnH4OpxATkdrJ7YZuPr0PJ/691CacAFfujvztb39zGE6Ay3exS0tLkZeX16N6fQkDClnJ7Zo5ffq009uQH2zMQUPMj/Gjmc9AE2q7wZYmNNI6fdi6CVdopMM2RERqJrcbuqlwl93jlrvMV/5R55w7x/ipHbt4yKqrrhnL9GAATpK+gMDwSMTM/afs6cGcQkxEvkoyizCsmeN0xeuA/mEwX6p3y+dZuth9Fbt4qFs6ds1YFlSzsDx/6KGHnN6GBCSYjNVo1Rdeft1VG2zZCx5y2hARqZEQoEHEhLlO2wwYJW/ZhbDwQZ3+7rV+jiAgPj4eaWlpLtfoqxhQyEZ6ejpycnKsd0ostFotcnJycN1118l6H04PJqLeosvu6uuSZb2PMHrKD3ev7f8DcdWqVT49QNZVXOKOOklPT8fdd99td4bO7t27Zb0HpwcTUW8SnJiK/tcl2+2ulsyirI1Pw1PuRWBkQqcF37RaLVatWmWzUKY/L+ZmwYBCdmk0mk79nKZ2M873HYbAsEiY6rnDMBFRR5buanvH5W58enXQiYgagj/fexfu+Pk11rY9Xe3bV3CQrJ+Rk6pdTd5H9XX4vLASHx8ywFB3ybpqoiOcgUNE1FlPNj7tqxGQem0kBlcfwqrFDzvdiFXNIYV78fRSclK1nDbtohn7imvxeWEFdpyoRJmxpdNncYdhIiLX9WTWYlczhq5ezE2N3UAMKL2QZYlkZ6kagNM2//PyP9AUOw5fFlWhrrmty8/k9GAiIu+Ru2rtrl27UFtbq8puIAYUH9STpGtZv8TZEsldr1/CpeWJiNSs6cQe1Gx5sct24zN+h90fvSWrG8jbd1lc+f3NQbJOeGI8hz1yBzw5+qy8vLwul0h2vnbJD+/fUINWfaHdQV5ERKQsubMj9/xH53DTV0EQkJWVhbvvvhubN29W5V0WCwYUB9w1ngNwHmIcdc1Ylo23JF1nn9Xa2uq26+b6JURE6hSkHdXldOWuVq217Onz4KNP4r01/9fl7x4lKdrFs3r1arz44ouoqKjA2LFj8eqrr+KWW27p8nWe7uJxx3gOOcFC7u6VL7/8Mn7zm990TsSCAEjAbfc9jNz3X+/RNVtEz1zOOyhERCrV1SzK0J/cjYYDm7t8n4B+ITC3NNo958mdk31iDMoHH3yABx98EGvWrEFycjJWrVqFjRs3oqioCEOGDHH6Wk8GFHeM5+gqWFhCzNKlS7FkyZIuawodGIGGulqH5wNCBkMAIDZe6FEbjkEhIlI/Z7MoA/qHyBpIK4cn9v3xiYCSnJyMn/70p3jttdcAAGazGfHx8Xj00Ufx5JNPOn2tJwPK7t27MX68vH0TuhIVFYXq6mq75wRBwMBBg3Cx1nHwcEX4z++H8at3Hdcy/fL/sFy/hIjI9zmaRSln80KhXyikloYuP+Pdd9/FzJkz3Vm2+jcLNJlMOHjwICZOnHilkIAATJw4EQUFBZ3at7a2or6+3ubhKe7cytpROAEu9wO6K5wAQN9BQ53vBZGY2vV+EQwnREQ+wdEmq3I2Lwwb9ytZnxEbG9vjOntCkUGyNTU1EEUR0dHRNsejo6Nx6tSpTu1XrFiBZ5991iu1CcEDvfI5FsGh4WhuqAdg70aWgAHhg9Bk7DrI3D9+LK6/KRnmzAdx5tgB1F+oQtjgIfjR6J8goGMf4uTErtsQEZHvmpyIwzfFIWf1/6KuusJ6eGBULDIy/wdjfjYRz9y/E3XVlXD0u2dIbJziOycr0sVTVlaGoUOHIj8/HykpKdbj//3f/409e/Zg3759Nu1bW1ttZqrU19cjPj7eo2NQDAaD3WlaHcegOGsTGRnp9A6KxbPPPoulS5cCgM17WcapfPDBB1i4cKHTz/LUYCYiIvJdcmaQAvZ/93hqFo/qu3giIyOh0WhQWVlpc7yyshIxMTGd2gcFBSEsLMzm4SkajQbZ2dkArnxRFpbn2dnZXbZZvXo1tFptp/Md28XHx+N//ud/kJOTYw09FlqtFjk5Objnnnu6/KzetgU3ERF1zbLp68yZM3H77bfb/J5IT093+rtH6SnGAABJIbfccov0yCOPWJ+LoigNHTpUWrFiRZevNRqNEgDJaDR6rL6PPvpI0mq1Ei7f/5IASPHx8dJHH30ku81HH30kCYIgCYJg08ZyrON7tbe3S7t27ZLeffddadeuXVJ7e7vL9RAREbmiq9897ubK729FpxnPnj0bb775Jm655RasWrUKH374IU6dOtVpbMrVvLXUvTtWkrW3Dkp8fDxWrVrlckJV48ZPREREcvnENGMAeO2116wLtd1000145ZVXkJyc3OXrfG0vHgYLIiIiHwoo3eVrAYWIiIh8YJAsERERkTMMKERERKQ6DChERESkOgwoREREpDoMKERERKQ6DChERESkOgwoREREpDoMKERERKQ6DChERESkOn2ULqA7LIvf1tfXK1wJERERyWX5vS1nEXufDCgNDQ0ALm+6R0RERL6loaEB4eHhTtv45F48ZrMZZWVlCA0NhSAIbnvf+vp6xMfHo7S01G/3+PH3a/T36wP8/xp5fb7P36/R368P8Nw1SpKEhoYGxMXFISDA+SgTn7yDEhAQAK1W67H3DwsL89v/6Sz8/Rr9/foA/79GXp/v8/dr9PfrAzxzjV3dObHgIFkiIiJSHQYUIiIiUh0GlA6CgoKwZMkSBAUFKV2Kx/j7Nfr79QH+f428Pt/n79fo79cHqOMafXKQLBEREfk33kEhIiIi1WFAISIiItVhQCEiIiLVYUAhIiIi1el1AeVvf/sbUlNTERwcjIEDB9ptU1JSgmnTpiE4OBhDhgzBE088gfb2dqfvW1tbi1mzZiEsLAwDBw7EnDlz0NjY6IErcM3u3bshCILdx/79+x2+7vbbb+/Uft68eV6sXL6EhIROtT7//PNOX9PS0oLMzEwMHjwYISEhmDFjBiorK71UsXznzp3DnDlzMGLECPTv3x/XXnstlixZApPJ5PR1av/+Vq9ejYSEBPTr1w/Jycn45ptvnLbfuHEjRo4ciX79+mH06NHYtm2blyp13YoVK/DTn/4UoaGhGDJkCKZPn46ioiKnr9mwYUOn76tfv35eqtg1S5cu7VTryJEjnb7Gl74/e3+fCIKAzMxMu+194bvLzc3FL3/5S8TFxUEQBGzatMnmvCRJeOaZZxAbG4v+/ftj4sSJOH36dJfv6+rPsat6XUAxmUy455578PDDD9s9L4oipk2bBpPJhPz8fLz11lvYsGEDnnnmGafvO2vWLBQWFmLHjh3YunUrcnNzMXfuXE9cgktSU1NRXl5u8/jjH/+IESNG4Cc/+YnT1z700EM2r3vhhRe8VLXrli1bZlPro48+6rT9ggULsGXLFmzcuBF79uxBWVkZ0tPTvVStfKdOnYLZbMabb76JwsJCrFy5EmvWrMFTTz3V5WvV+v198MEHWLhwIZYsWYJvv/0WY8eOxeTJk1FVVWW3fX5+PmbOnIk5c+bg0KFDmD59OqZPn47jx497uXJ59uzZg8zMTOzduxc7duxAW1sbJk2ahKamJqevCwsLs/m+zp8/76WKXTdq1CibWr/66iuHbX3t+9u/f7/Nte3YsQMAcM899zh8jdq/u6amJowdOxarV6+2e/6FF17AK6+8gjVr1mDfvn0YMGAAJk+ejJaWFofv6erPcbdIvdT69eul8PDwTse3bdsmBQQESBUVFdZjb7zxhhQWFia1trbafa8TJ05IAKT9+/dbj3366aeSIAiSwWBwe+09YTKZpKioKGnZsmVO2/3iF7+Q5s+f752iemj48OHSypUrZbevq6uT+vbtK23cuNF67OTJkxIAqaCgwAMVutcLL7wgjRgxwmkbNX9/t9xyi5SZmWl9LoqiFBcXJ61YscJu+9/85jfStGnTbI4lJydLf/rTnzxap7tUVVVJAKQ9e/Y4bOPo7yM1WrJkiTR27FjZ7X39+5s/f7507bXXSmaz2e55X/ruJEmSAEgff/yx9bnZbJZiYmKkF1980Xqsrq5OCgoKkt577z2H7+Pqz3F39Lo7KF0pKCjA6NGjER0dbT02efJk1NfXo7Cw0OFrBg4caHNHYuLEiQgICMC+ffs8XrMrPvnkE1y4cAG///3vu2z7zjvvIDIyEjfeeCMWL16M5uZmL1TYPc8//zwGDx6Mm2++GS+++KLTLrmDBw+ira0NEydOtB4bOXIkhg0bhoKCAm+U2yNGoxERERFdtlPj92cymXDw4EGb//YBAQGYOHGiw//2BQUFNu2Byz+TvvBdAZe/LwBdfmeNjY0YPnw44uPjcffddzv8+0YNTp8+jbi4OFxzzTWYNWsWSkpKHLb15e/PZDLh7bffxh/+8AenG9P60nd3teLiYlRUVNh8R+Hh4UhOTnb4HXXn57g7fHKzQE+qqKiwCScArM8rKiocvmbIkCE2x/r06YOIiAiHr1HKunXrMHny5C43W7z//vsxfPhwxMXF4ejRo1i0aBGKioqg0+m8VKl8jz32GH784x8jIiIC+fn5WLx4McrLy/Hyyy/bbV9RUYHAwMBOY5Cio6NV931d7cyZM3j11Vfx0ksvOW2n1u+vpqYGoija/Rk7deqU3dc4+plU+3cFXN55PSsrCz/72c9w4403OmyXmJiIf/3rXxgzZgyMRiNeeuklpKamorCw0KMbo3ZHcnIyNmzYgMTERJSXl+PZZ59FWloajh8/jtDQ0E7tffn727RpE+rq6vC73/3OYRtf+u7ssXwPrnxH3fk57g6/CChPPvkk/v73vzttc/LkyS4HcvmS7lyzXq/H9u3b8eGHH3b5/h3Hz4wePRqxsbGYMGECzp49i2uvvbb7hcvkyvUtXLjQemzMmDEIDAzEn/70J6xYsUK1S1F35/szGAy48847cc899+Chhx5y+lqlvz+6LDMzE8ePH3c6RgMAUlJSkJKSYn2empqKG264AW+++Saee+45T5fpkilTplj/PGbMGCQnJ2P48OH48MMPMWfOHAUrc79169ZhypQpiIuLc9jGl747X+MXAeXxxx93mnAB4JprrpH1XjExMZ1GIltmd8TExDh8zdUDg9rb21FbW+vwNT3VnWtev349Bg8ejF/96lcuf15ycjKAy/+C98YvuJ58p8nJyWhvb8e5c+eQmJjY6XxMTAxMJhPq6ups7qJUVlZ67Pu6mqvXV1ZWhvHjxyM1NRVr1651+fO8/f05EhkZCY1G02nGlLP/9jExMS61V4tHHnnEOmDe1X9J9+3bFzfffDPOnDnjoercZ+DAgbj++usd1uqr39/58+fxxRdfuHzX0Ze+O+DK77XKykrExsZaj1dWVuKmm26y+5ru/Bx3i9tGs/iYrgbJVlZWWo+9+eabUlhYmNTS0mL3vSyDZA8cOGA9tn37dlUNkjWbzdKIESOkxx9/vFuv/+qrryQA0pEjR9xcmfu9/fbbUkBAgFRbW2v3vGWQbE5OjvXYqVOnVDtIVq/XS9ddd5103333Se3t7d16DzV9f7fccov0yCOPWJ+LoigNHTrU6SDZu+66y+ZYSkqKagdZms1mKTMzU4qLi5O+++67br1He3u7lJiYKC1YsMDN1blfQ0ODNGjQICk7O9vueV/7/iyWLFkixcTESG1tbS69Tu3fHRwMkn3ppZesx4xGo6xBsq78HHerVre9k484f/68dOjQIenZZ5+VQkJCpEOHDkmHDh2SGhoaJEm6/D/XjTfeKE2aNEk6fPiw9Nlnn0lRUVHS4sWLre+xb98+KTExUdLr9dZjd955p3TzzTdL+/btk7766ivpuuuuk2bOnOn163Pkiy++kABIJ0+e7HROr9dLiYmJ0r59+yRJkqQzZ85Iy5Ytkw4cOCAVFxdLmzdvlq655hrptttu83bZXcrPz5dWrlwpHT58WDp79qz09ttvS1FRUdKDDz5obXP19UmSJM2bN08aNmyY9OWXX0oHDhyQUlJSpJSUFCUuwSm9Xi/96Ec/kiZMmCDp9XqpvLzc+ujYxpe+v/fff18KCgqSNmzYIJ04cUKaO3euNHDgQOvMuQceeEB68sknre2//vprqU+fPtJLL70knTx5UlqyZInUt29f6dixY0pdglMPP/ywFB4eLu3evdvm+2pubra2ufoan332WWn79u3S2bNnpYMHD0r33Xef1K9fP6mwsFCJS3Dq8ccfl3bv3i0VFxdLX3/9tTRx4kQpMjJSqqqqkiTJ978/Sbr8y3bYsGHSokWLOp3zxe+uoaHB+rsOgPTyyy9Lhw4dks6fPy9JkiQ9//zz0sCBA6XNmzdLR48ele6++25pxIgR0qVLl6zvcccdd0ivvvqq9XlXP8fu0OsCyuzZsyUAnR67du2ytjl37pw0ZcoUqX///lJkZKT0+OOP26ToXbt2SQCk4uJi67ELFy5IM2fOlEJCQqSwsDDp97//vTX0qMHMmTOl1NRUu+eKi4tt/huUlJRIt912mxQRESEFBQVJP/rRj6QnnnhCMhqNXqxYnoMHD0rJyclSeHi41K9fP+mGG26Qli9fbnO36+rrkyRJunTpkvTnP/9ZGjRokBQcHCz9+te/tvmlrxbr16+3+/9rx5ufvvj9vfrqq9KwYcOkwMBA6ZZbbpH27t1rPfeLX/xCmj17tk37Dz/8ULr++uulwMBAadSoUdJ//vMfL1csn6Pva/369dY2V19jVlaW9b9HdHS0NHXqVOnbb7/1fvEy3HvvvVJsbKwUGBgoDR06VLr33nulM2fOWM/7+vcnSZfvgAOQioqKOp3zxe/O8jvr6oflOsxms/TXv/5Vio6OloKCgqQJEyZ0uvbhw4dLS5YssTnm7OfYHQRJkiT3dRgRERER9RzXQSEiIiLVYUAhIiIi1WFAISIiItVhQCEiIiLVYUAhIiIi1WFAISIiItVhQCEiIiLVYUAhIiIi1WFAISIiItVhQCEiIiLVYUAhIiIi1WFAISIiItX5/7YMGJZdDrL9AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from scipy.stats import truncnorm\n", "from iminuit.cost import UnbinnedNLL\n", "import numpy as np\n", "\n", "xrange = (-10., 10.)\n", "\n", "rng = np.random.default_rng(1)\n", "x = rng.normal(1, 3, size=10000)\n", "x = x[(xrange[0] < x) & (x < xrange[1])]\n", "\n", "def model(x, mu, sigma):\n", " zrange = np.subtract(xrange, mu) / sigma\n", " return truncnorm.pdf(x, *zrange, mu, sigma)\n", "\n", "# better use numba_stats.truncnorm, which is simpler to use and faster\n", "#\n", "# from numba_stats import truncnorm\n", "#\n", "# def model(x, mu, sigma):\n", "# return truncnorm.pdf(x, *xrange, mu, sigma)\n", "\n", "nll = UnbinnedNLL(x, model)\n", "m = Minuit(nll, 1, 3)\n", "m.migrad()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We do not get the exact same fitted values as before, since the data sample is different from the one generated by RooFit.\n", "\n", "To get the exact same result, we need to convert the variable `data` which has the type `RooDataSet` into a numpy array. The ROOT Python layer offers the method `to_numpy` for this purpose." ] }, { "cell_type": "code", "execution_count": 7, "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 = 5.027e+04 Nfcn = 31
EDM = 5.43e-08 (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", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 mu 1.003 0.030
1 sigma 3.017 0.022
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
mu sigma
mu 0.000926 0 (0.030)
sigma 0 (0.030) 0.000497
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 5.027e+04 │ Nfcn = 31 │\n", "│ EDM = 5.43e-08 (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 │ mu │ 1.003 │ 0.030 │ │ │ │ │ │\n", "│ 1 │ sigma │ 3.017 │ 0.022 │ │ │ │ │ │\n", "└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌───────┬───────────────────┐\n", "│ │ mu sigma │\n", "├───────┼───────────────────┤\n", "│ mu │ 0.000926 0 │\n", "│ sigma │ 0 0.000497 │\n", "└───────┴───────────────────┘" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNiUlEQVR4nO3de1xUdf4/8NcwCojcRJCLoNhlMa+1tRHsF7voV1PbXy5O28Utvy2b1WJBtq2r26bZd7W1viWVpe66WdtdGrXIS2oqFGiKmopKaigwMFxEboMwcOb8/mAZGWFuMDPnzMzr+XjM49E55zPD+3QE3nwu749CFEURRERERDLiI3UARERERFdjgkJERESywwSFiIiIZIcJChEREckOExQiIiKSHSYoREREJDtMUIiIiEh2mKAQERGR7AyQOoC+MBgMqKioQFBQEBQKhdThEBERkQ1EUURTUxNiYmLg42O5j8QtE5SKigrExcVJHQYRERH1QVlZGWJjYy22ccsEJSgoCEDnDQYHB0scDREREdmisbERcXFxxt/jlrhlgtI1rBMcHMwEhYiIyM3YMj2Dk2SJiIhIdpigEBERkewwQSEiIiLZYYJCREREssMEhYiIiGSHCQoRERHJDhMUIiIikh0mKERERCQ7TFCIiIhIdpigEBERkewwQSEiIiLZYYJCREREssMEhYiIiGSHCQoRuT2dTgeFQgGFQgGdTid1OETkAExQiIiISHaYoBAREZHsMEEhIiIi2WGCQkRERLLDBIWIiIhkhwkKEXkFW1b6cDUQkXwwQSEiIiLZYYJCREREssMEhYiIiGSHCQoRuT1BEIz/nZuba3JMRO6JCQoRuTW1Wo0xY8YYj2fMmIH4+Hio1WoJoyKi/mKCQkSyZmlljVqthkqlgkajMTmv0WigUqmYpBC5MSYoROSWBEFARkYGRFHsca3rXGZmJod7iNwUExQickt5eXkoLy83e10URZSVlSEvLw+AbfNUOJeFSD6YoBCRW6qsrLS53eeff44bbrA8T4VzWYjkZYDUARAR9UV0dLRN7Z7953ZUfvN+j/Nd81Sys7MBACqVqsdwUfc2qamp/Q+aiGymEHsbwJW5xsZGhISEoKGhAcHBwVKHQ0ROpNPpEBgYCABobm7G4MGDAXQOx8THx0Oj0fQ6DwUAfAKHQgFAaL5o9vOHDIvCoIFKVFw10baLQqFAbGwsSkpKoFQq+3czRF7Ont/fHOIhIrekVCqRlZVlsU3QxLstJicAcKlaazY5AXrOZSEi12CCQkRuK/amOzDm4aVQBg41Oa8MCkfErMUYGBbjsK9l65wXgJsOEjkC56AQkdu5rBfw9+2n8V7BeYjRNyM67W2UZ90PAIhQLcWgUTdB4aNEa+kxh31NW+e8EJFjMEEhIlm7eulv0LU3Y/HmIpTWtRjPK3yuzA3xjxtnPPaLHQtlUDiEplqzn29tnkrXHJSUlJR+3gkR2YNDPEQkW70t/b3jlrE4XbDLpvcrfJQImzzPYpuhUx5H2JTHzX0CAGDVqlWcIEvkYkxQiEiWzJWxF5pqUbN5OVqK8236nICEZETMWmx2nkpAQrKFNkMx85lXuMSYSAJMUIjIKfozUdRSGfsudbvXQTTYVuk1ICEZ0WlvG48jVEsx/In1CEhIttrm+MDReHP3GbviJ6L+sytBWbp0qfEHTtdr9OjRxuutra1IT0/H0KFDERgYiNmzZ6OqqsrkM0pLSzFz5kwEBARg2LBheO6559DR0eGYuyEit2IuibFWxh7o7ElpKy+y+WuZm6diS5vXdv2InGMVNn8tIuo/uyfJjh07Frt2XRn/HTDgykc888wz+Oqrr7Bx40aEhIRg/vz5SE1NxXfffQeg86+imTNnIioqCvn5+aisrMQjjzyCgQMHYvny5Q64HSLyBLYu6RWaLzk5kk6iCDz72Q8YHjoIN40Y4pKvSeTt7B7iGTBgAKKiooyv8PBwAEBDQwPWr1+P1157DXfddRduvvlmvPvuu8jPz8f+/fsBAF9//TVOnjyJDz74ADfeeCOmT5+Ol156CatXr4Zer3fsnRGR27J1Sa8y0HXJQluHAY+9XwhN/WWrbbnpIFH/2Z2gnDlzBjExMbjmmmswZ84clJaWAgAKCwvR3t6OKVOmGNuOHj0aI0aMQEFBAQCgoKAA48ePR2RkpLHNtGnT0NjYiKIi8121bW1taGxsNHkRkee6LfmXCBgyzGIbZVA4/GLHuiiiTrXNbUjbcBDNbeaHpbnpIJFj2JWgJCYmYsOGDdi+fTveeecdlJSUICUlBU1NTdBqtfD19UVoaKjJeyIjI6HVagEAWq3WJDnput51zZwVK1YgJCTE+IqLi7MnbCJyM4s3FWHw7b+32CZs8jzjHBEfX3+MXJiDkQtz4OPr79TYTmub8NRHhyEYek7gNbfyqGvTQSYpRLazK0GZPn067rvvPkyYMAHTpk3D1q1bUV9fj88++8xZ8QEAFi1ahIaGBuOrrKzMqV+PiKSTXVgO9RGNTcuDpWDQt2LD7xIxQOljMrHX0sqjrnOZmZkc7iGyUb+WGYeGhuJnP/sZzp49i6ioKOj1etTX15u0qaqqQlRUFAAgKiqqx6qeruOuNr3x8/NDcHCwyYuIPE91YyteyjlpPLZlebCUPi+8stLI2sojbjpIZJ9+JSjNzc04d+4coqOjcfPNN2PgwIHYvXu38XpxcTFKS0uRlJQEAEhKSsLx48dRXV1tbLNz504EBwebjNkSkfvry0TRv2w+gYbL7SbnbFkeLJWV20+juqkVgO0rj+zZdJDIm9mVoPzxj3/Evn37cP78eeTn5+PXv/41lEolHnzwQYSEhCAtLQ0LFizAnj17UFhYiEcffRRJSUm47bbbAABTp07FmDFj8PDDD+OHH37Ajh078PzzzyM9PR1+fn5OuUEicr2+TBTNOVaBnSerzF7vL1vmqdg7l6WxtQMvftHZ42PryiNuOkhkG7sSlPLycjz44INISEjAb37zGwwdOhT79+9HREQEAOD111/HPffcg9mzZ2PSpEmIiooy+YGkVCqRk5MDpVKJpKQk/Pa3v8UjjzyCZcuWOfauiEgy9kwU7d6rsvjtz2yuDCsnXx2vxDenq5CSkoLY2FgoFIpe2ykUCsTFxXHTQSIbKURLtaRlqrGxESEhIWhoaOB8FCIZEQQB8fHxZudidO0MXFJSgi1btuDpp582SWSUQeEImzzPZI6JQd+KstdVAIC4Z7KdvkrHmt7iGR46CDsXTML2nC+gUnVe6/6jtStpyc7O5r4+5NXs+f3NvXiIyGFsnSj6t7/9zSEbAcqFpv4y/u/rH5Gamors7GzExMSYXI+Nje2RnPRnryIib8AEhYgcxtYJoFlZWQ7bCFAuNuSfx/HyBqSmpuLkySsrkbZu3YqSkhL2nBDZiQkKETmMrRNA6+rqLF63dyNAV+qeOLWWnTAeCwYRizYdg2AQoVReWWk0adIkk2Misg0TFCJyGFsmioaFhdn0WV0bAbqySqw1LcX5qFz/B+NxTfZSaNakGYekTmga8e53JVKFR+RRmKAQkcMolUpkZWUBQI8kpes4IyPDts9y4UaAtmgpzkfN5uUQmi+anL963sxrO3+E5lKLFCESeRQmKETkUNYmiv7lL39BaLj5ytGANBsBWiIaBNTtXmexTde8mRa9gJe+OuWiyIg8FxMUInI4SxNFG1sFBN1p+0aActBWXgShqdZim+7zZvYV17giLCKPxgSFiJzC3ETRt/eehc81t8lyI0BzuubDOKodEVnHBIWIXKay4TLeL7gAQP4bAXZn63wYe+bN9GWvIiJvwgSFiFwma9cZtHUYjMdy3giwO7/YsVAGhVtsY27eTLtg6HGuL3sVEXkbJihE5BLnapqxsdB8lVk5U/goETZ5nsU23efNdF8anXPStOaLPXsVEXkzJihE5BKvff0jBIPbbf1lFJCQ3Kd5M299cxat7f8p5iYIyMjI6LWKbte5zMxMDvcQgQkKEbnACU0Dtp6wrQy+nPVl3oy2sRUfHigFYPteRXl5eY4LmshNMUEhIqdbtetHuN++6b3ry7yZd/aew2W9YPNeRba2I/JkTFCIyOm+O3vReiMPVtvchg35523eq8jWdkSejAkKETnF4MGDIYoi7n3rW8n30JGDtbnncOMvbrO6V1FcXBxSUlJcHB2R/DBBISKn2X5Ci6Nl9Wavy2kjQGerb2nHhoJSq3sVrVq1irsfE4EJChE5icEg4v++LpY6DFlZ/20JJt/9K4t7FaWmpkoUHZG8MEEhIqf4+qQWZ6qbpQ5DVppaO7A295zFvYqIqBMTFCJyirW5P0kdgixtyD+Pi81tZvcqIqJOTFCIyOEOna/DkdJ6qcOQpRa9gHV5TN6IrBkgdQBE5Hk8ufeka2Jvf3x8oBRpiTHWGxJ5MfagEJFD/VTTjN2nqqQOQ9YaWzugPuye+xIRuQoTFCJyqH/klcCNt9xxmff3n5c6BCJZY4JCRA5T29zGngEblde1Sh0CkawxQSEih3k//zzaOgxSh0FEHoAJChE5xGW9gH/vvyB1GG6ja7Lt4Qt1GDx4sNThEMkOExQicoiNhWW41NIudRhu55/flkgdApEsMUEhon4zGESs5y/aPtl+QovySy1Sh0EkO0xQiKjfthdpceEif8n2hWAQseG781KHQSQ7TFCIqN/WeXBhNlf49GAZmts6pA6DSFaYoBBRv3xfUoejZfVSh+HWmto68Mn3pVKHQSQrTFCIqF/W5Z6TOgSPsCH/PARWuCMyYoJCRH12vlaH3aerpQ7DI5RfuoztJ7RSh0EkG0xQiMhuOp0OCoUCoyICIbSxIqqj/PNbzuUh6sIEhYhIJo6U1qPwwiWpwyCSBSYoREQy8l7+ealDIJIFJihERDKy/YQWdTo9gCtDaQqFAjqdTuLIiFyLCQoR2e20tlHqEDyWXjAgu7BM6jCIJMcEhYjsxpodzvXx92UQRS45Ju/GBIWI7KJr68BXx7gc1plKanUoOHdR6jCIJMUEhYjssumIBk2X24zHrWUnIBoECSPyTB+yl4q8HBMUIrLL/637NyrX/8F4XJO9FJo1aWgpzpcwKs/zdZEWtc1t1hsSeSgmKERepL+rQl5Z+z6OvvtXCM2mww9CUy1qNi9nkuJA7YKITYc1UodBJBkmKERkE0EQsGzxnyy2qdu9jsM9DrTx0Hnjf+fm5kIQ+P+WvAcTFCKyybad36C5rspiG6GpFm3lRS6KyLO1FOfjwMpHjcczZsxAfHw81Gq1hFERuQ4TFCKyyZZ82xIPoZml2vurpTgfNZuX9xhK02g0UKlUTFLIKzBBISKrRFHEwRrb2ioDhzg3GA8nGgTU7V7X+7X/1EbJzMzkcA95PCYoRGTVd2cv4lLQNVAGhVtspwwKh1/sWBdF5ZnayosgNNWavS6KIsrKypCXl+fCqIhcjwkKEZnobaXPhwcuQOGjRNjkeRbfGzZ5HhQ+SleE6bFsHSKrrKx0ciRE0mKCQuRFug8L2LoqpLqpFTtPdk6ODUhIRsSsxVAGDjVpowwKR8SsxQhISHZswF7I1iGy6OhoJ0dCJK1+JSgvv/wyFAoFMjMzjedaW1uRnp6OoUOHIjAwELNnz0ZVlenM/9LSUsycORMBAQEYNmwYnnvuOXR0dPQnFCKyQq1WY8yYMcZjW1eFbD6iQYfhyr4wAQnJiE5723gcoVqK4U+sZ3LiIH6xYy0OpSkUCsTFxSElJcWFURG5Xp8TlIMHD2Lt2rWYMGGCyflnnnkGX375JTZu3Ih9+/ahoqICqampxuuCIGDmzJnQ6/XIz8/He++9hw0bNuCFF17o+10QkUVqtRoqlQoajWnhL1tWhXxe2LNYWPdhHP+4cRzWcSBLQ2kKhQIAsGrVKiiV/H9Onq1PCUpzczPmzJmDf/zjHxgy5Ep3ZENDA9avX4/XXnsNd911F26++Wa8++67yM/Px/79+wEAX3/9NU6ePIkPPvgAN954I6ZPn46XXnoJq1evhl6vd8xdEZGRIAjIyMjodXdca6tCiioaUFzV5PQYyZS5obSQ8ChkZ2eb/NFH5Kn6lKCkp6dj5syZmDJlisn5wsJCtLe3m5wfPXo0RowYgYKCAgBAQUEBxo8fj8jISGObadOmobGxEUVFvddZaGtrQ2Njo8mLiGyTl5eH8vJys9ctrQrZfISl1qXS21DauAXv49e//rWEURG5jt0JyieffILDhw9jxYoVPa5ptVr4+voiNDTU5HxkZCS0Wq2xTffkpOt617XerFixAiEhIcZXXFycvWETeS1bV3v01u6r41wpIqWrh9LK6ttwoKROwoiIXMeuBKWsrAwZGRn48MMP4e/v76yYeli0aBEaGhqMr7KyMpd9bSJ3Z+tqj6523Yd6Kk8f4d46MpNdaL43jMiT2JWgFBYWorq6Gj//+c8xYMAADBgwAPv27cMbb7yBAQMGIDIyEnq9HvX19Sbvq6qqQlRUFAAgKiqqx6qeruOuNlfz8/NDcHCwyYuIbJOSkoLY2FjjBMurdV8VcvVKn5rspdCsSeMuxTKy7XglWvRc9Uiez64EZfLkyTh+/DiOHj1qfN1yyy2YM2eO8b8HDhyI3bt3G99TXFyM0tJSJCUlAQCSkpJw/PhxVFdXG9vs3LkTwcHBJj8YicgxlEolsrKyAKBHktJ9VciWLVt6XekjNNWiZvNyJikyodML2Hq89+FwIk8ywJ7GQUFBGDdunMm5wYMHY+jQocbzaWlpWLBgAcLCwhAcHIynnnoKSUlJuO222wAAU6dOxZgxY/Dwww9j5cqV0Gq1eP7555Geng4/Pz8H3RYRdZeamors7Gw8/fTTJglIbGwsVq1ahXvvvRfx8fG9rvTpUrd7HQZdnwiFjxI+vv4YuTDHFaFTLzYeKoPq5lipwyByKodXkn399ddxzz33YPbs2Zg0aRKioqJMaiwolUrk5ORAqVQiKSkJv/3tb/HII49g2bJljg6FiLpJTU3FyZMnjcdbt25FSUkJUlNTra70ATp7UtrKbdvRmJzr+/N1KKtrkToMIqeyqwelN3v37jU59vf3x+rVq7F69Wqz7xk5ciS2bt3a3y9NRHbqXtxr0qRJxmNbV/rYuk8MOZcodk6Wfea/fwadTofAwEAAnTWqBg8eLHF0RI7R7wSFiNyfrSt9bN0nhhzD0lDa54fLkTnlehdHROQ63CyQiKyu9AE6NwT0ix3rwqjIkvJLl1Hw00WpwyByGiYoRGRc6WN+iiwQNnke99yRGdZEIU/GBIWIAHROop2esbLH/i/KoHBEzFrM3YplaPsJLXRtrIlCnolzUIgIANDU2o6fAsciOu1tlGfdD6Bz/5dBo25iz4lMtegFbD/BmijkmdiDQkQAgK+OVaK13dBj/xcmJ/LGDR3JUzFBIfIigwcPhiiKEEWxx3LUzw9zPoM7OnS+1vjfubm5JnspEbkzJihEhAsXdTh4njVO3E1LcT4q/vkH4/GMGTMQHx9vUhyTyF0xQSEifH6YwwTupqU4HzWbl0NoNl1qrNFooFKpmKSQ22OCQuTlRFGEmsM7bkU0CKjbva73a//ZTykzM5PDPeTWmKAQebn9P9Wh/NJlqcMgO7SVF0FoqjV7XRRFlJWVIS8vz4VRETkWlxkTebmrJ8dyp2L5s3VPJFv3WCKSI/agEHmxFn0Hth3nLzF3Y+ueSLbusUQkR0xQiLzY9hNa6PScp+Bu/GLHQhkUbva6QqFAXFwcUlJSXBgVkWMxQSHyYqx94p4UPkqETZ7X+7X/bPi4atUqKJUsskfuiwkKkZeqqL+MgnPcDdddBSQkI2LW4h57J4VGRCE7OxupqakSRUbkGJwkS+SlNh3RwGBp+2KSvYCEZPiNnGiyd9Itv7wdqal3SBsYkQOwB4XIS31eyOEdT3D13kkntTqcqWqSMCIix2CCQuSFDpdewk+1OqnDICdRcwNB8gBMUIi8EHtPPNsXRyuMFWWJ3BUTFCIv09YhIOcYa594Mk39Zez/qU7qMIj6hQkKkZfZdbIaDZfbpQ6DnGzTEfaSkXtjgkLkZVj7xDtsO6FFazuL8JH74jJjIi9S09SG3B9rpA6DHMjc3klNrR3YdaoK90yIkSAqov5jDwqRF9lyVIMOFj/xGpu5mofcGBMUIg+g0+mgUCigUCig05lfPpzN1TteZd+PNajT6aUOg6hPmKAQeYmTFY04rWUBL2/SLojIOVYhdRhEfcIEhchLbD7K7n5vpD7M507uiQkKkRcwGER8cZR/SXujo2X1KGHVYHJDTFCIvMD+ny5C29gqdRgkkU2cLEtuiAkKkRfg8I5328LnT26ICQqRBxCEKwW5cnNzTY7bOgRsO6GVIiySiQsXW3C49JLUYRDZhQkKkZtTq9UYM2aM8XjGjBmIj4+HWq0GAOw+VY2m1g6pwiMZMOhbcfPIMKvL0InkhJVkidyYWq2GSqXqsXOtRqOBSqVCdnY2tjePkCg6IqK+Yw8KkZsSBAEZGRk9khMAxnNPZ2RgzykO7xCR+2GCQuSm8vLyUF5uvjKsKIrQlJej6cJxF0ZFROQYTFCI3FRlZaVN7YRmTo6kK1r0nI9E7oEJCpGbio6OtqmdMnCIkyMhd7LrVJXUIRDZhAkKkZtKSUlBbGwsFAqFmRYKKIPC4Rc71qVxkfyIhivLzv+5cavJMnQiuWKCQuSmlEolsrKyAKBHktJ5LCJs8jwofJQSREdy0VKcj8r1fzAef7PqGYwYOdK4DJ1IrpigELmx1NRUZGdnIyYmxuR8ZHQMImYtRkBCskSRkRy0FOejZvNyCM0XTc5XaCqgUqmYpJCsKcTe1ijKXGNjI0JCQtDQ0IDg4GCpwyGSXNf3BABs3boVh4U4/OPbCxJHRVISDQI0a9IgNNX2el2hUCA2NhYlJSVQKtnLRq5hz+9v9qAQeYDuv2BSUlLw1XFOhPR2beVFZpMToHMZellZGfLy8lwYFZHtmKAQeZjvS+pQ0cCdi72drcvLbV2uTuRqTFCIPEzOsQqpQyAZsHV5eddydZ1OB4VCwf16SDaYoBB5mK+LOLxDgF/sWCiDwi22iYwejpSUFBdFRGQfJihEHqaROxcTAIWPEmGT51lsc/vcP3KCLMkWExQiDzB48GCIoojH3z8EH19/qcMhmQhISEbErMVQBg41Oa8MCkfErMU4N3gsOgSDRNERWTZA6gCIyDEaLrfjm+JqqcMgmQlISIbfyIkoz7ofABChWopBo26CwkeJ2mY9vj1bizsShkkcJVFP7EEh8hDbT1RC38G/hqmn7tWE/ePGmRxvPqIBAJPy97m5uSyHT5KzK0F55513MGHCBAQHByM4OBhJSUnYtm2b8XprayvS09MxdOhQBAYGYvbs2aiqMp2wV1paipkzZyIgIADDhg3Dc889h44OjpkT9dfmI1y9Q/b7+mQVPv5sI8aMGWM8N2PGDMTHx7PSLEnKrgQlNjYWL7/8MgoLC3Ho0CHcdddduPfee1FUVAQAeOaZZ/Dll19i48aN2LdvHyoqKpCammp8vyAImDlzJvR6PfLz8/Hee+9hw4YNeOGFFxx7V0ReprLhMg6UXLTekOgqtcfz8NAD90Oj0Zic12g0LIdPkup3qfuwsDC88sorUKlUiIiIwEcffQSVSgUAOH36NG644QYUFBTgtttuw7Zt23DPPfegoqICkZGRAIA1a9Zg4cKFqKmpga+vr01fk6XuiUyt3XcOK7adljoMkimDvhVlr3f+XI57Jts4kZrl8MnVXFLqXhAEfPLJJ9DpdEhKSkJhYSHa29sxZcoUY5vRo0djxIgRKCgoAAAUFBRg/PjxxuQEAKZNm4bGxkZjL0xv2tra0NjYaPIiois2H+XwDtmP5fBJzuxexXP8+HEkJSWhtbUVgYGB2LRpE8aMGYOjR4/C19cXoaGhJu0jIyOh1WoBAFqt1iQ56bredc2cFStW4MUXX7Q3VCKvUKxtwqlKJu1kno+vP0YuzOlxnuXwSc7s7kFJSEjA0aNHceDAATz55JOYO3cuTp486YzYjBYtWoSGhgbjq6yszKlfj8idbD6qsd6IqBf2lsMnciW7e1B8fX1x3XXXAQBuvvlmHDx4EFlZWbj//vuh1+tRX19v0otSVVWFqKgoAEBUVBS+//57k8/rWuXT1aY3fn5+8PPzszdUIo8niiK+4PAO9VFXOXxrc1BYDp+k0O86KAaDAW1tbbj55psxcOBA7N6923ituLgYpaWlSEpKAgAkJSXh+PHjqK6+Ukxq586dCA4ONlniRkRXWNrE7eD5S9DUX5YoMnJ3lsrhKxQKAMCqVas4QZYkYVcPyqJFizB9+nSMGDECTU1N+Oijj7B3717s2LEDISEhSEtLw4IFCxAWFobg4GA89dRTSEpKwm233QYAmDp1KsaMGYOHH34YK1euhFarxfPPP4/09HT2kBD1waYjHN6h/ukqh1+3ay2E5itL1WNjY7Fq1SqTUhFErmRXglJdXY1HHnkElZWVCAkJwYQJE7Bjxw7893//NwDg9ddfh4+PD2bPno22tjZMmzYNb7/9tvH9SqUSOTk5ePLJJ5GUlITBgwdj7ty5WLZsmWPvisgL6DsM2Hqckxep/64uh79161ZMnTqVPSckqX7XQZEC66CQN9HpdAgMDAQANDc3Y/DgwQCAr4u0mPfvQilDIw/SvVZK939nRI7kkjooRCStLZwcS04iGNzu71byQExQiGSut03cmlrbsetUlYV3EfVd/jnzxduIXIUJCpGMqdXqXjdxW/rGu2jjzsXkJF/9wLlNJD0mKEQypVaroVKpet3E7bWFj6OlOF+iyMjT7TpVhct6wXpDIidigkIkQ4IgICMjA73NYe86V7d7HUQDf4mQY3SVwx+5MAeXMRBfnzS//QiRKzBBIZKhvLw8lJeXW2wjNNWirdz8JptE/cFJ2CQ1JihEMmTr5my2bvZGZK/cH2tQp9NLHQZ5MSYoRDJk6+Zstm72RmSvDoOInGPsRSHpMEEhkqGUlBTExsYa90PpjTIoHH6xY10YFXkbbqVAUmKCQiRDSqUSWVlZAGA2SQmbPA8KH5YiJ+c5UlqPczXNUodBXooJCpFMpaamIjs7GzExMSbnlUHhiJi1GAEJyRJFRt4ku/DKZG1LO2sTORoTFCIZS01NxcmTJ43HEaqlGP7EeiYn5DKbDmtgYOl7kgATFCKZ676jrH/cOA7rkEtpG1uRd5al78n1mKAQEZFF3Yd5iFyFCQoREVn0dZEWja3tUodBXoYJCpHMaRsuSx0Cebm2DgNyfqjsdWdtImdhgkIkc1//2GDcI8XH11/qcMhLZa3/oNedtdVqtYRRkSdjgkIkc5uOcPyfpNVSnI/v//GXXnfWVqlUTFLIKZigEMnYCU0DfqxioSySjmgQULd7Xe/X/rOzdmZmJod7yOGYoBDJmPowS42TtNrKiyA0mV9mLIoiysrKkJeX58KoyBswQSGSqQ7BgC9+4GZtJC1bd8y2dQduIlsxQSGSqbwztahtbpM6DPJytu6YbesO3ES2YoJCJFOfH+bkWJKeX+xYKIPCzV5XKBSIi4tDSkqKC6Mib8AEhUiGmlrbsfNkldRhEEHho0TY5Hm9X/vPTturVq0y2ZKByBGYoBDJ0NbjlWjrMEgdBhEAICAhGRGzFkMZONTkfGxsLLKzs5GamipRZOTJmKAQydDGQxzeIXkJSEhGdNrbxuN/ffw5SkpKmJyQ0zBBIZKZn2qaceiCbSsniFyp+07aFYPiOaxDTsUEhUhmuHMsuYMvjlbAYBClDoM8GBMUIhkRDCKLs5FbqGxow3fnzBdwI+ovJihEMpJ7pgbaxlapwyCyCXv7yJkGSB0AEV2RzcmxJGM+vv4YuTDHeLyjSIum1nYE+Q+UMCryVOxBIZKJ+hY9dp5i7RNyH63tBuQcY4l7cg4mKEQyseVoBfSsfUJuhsM85CxMUIhkYmNhmdQhENmt8MIllNTqpA6DPBATFCIZOFXZiBOaRqnDIOqTbCbX5ARMUIhkgJVjyZ1tOqxhTRRyOCYoRBJrFwzYfJS1T8h9VTS0Iv/cRanDIA/DBIVIYrtPVaFOp5c6DKJ+4TAPORoTFCKJcXiHPMGOoio0tbZLHQZ5ECYoRBKqbmrF3h9rpA6DqN8utwv4ijVRyIGYoBBJaNNhDQROLiQPkV1YDp1OB4VCAYVCAZ2Oy4+p75igEEloI4tckQc5dOESLrAmCjkIExQiiRwpvYSz1c1Sh0HkUFyRRo7CBIVIIp98z1UP5Hm2MEEhB2GCQiSBptZ2fHmsQuowiByu4lKL8b9zc3MhCIKE0ZA7Y4JCJIHNRzRo0fMHN3mWluJ8VK7/g/F4xowZiI+Ph1qtljAqcldMUIgk8OGBUqlDIHKoluJ81GxeDqHZtKKsRqOBSqVikkJ2Y4JC5GKHSy/htLZJ6jCIHEY0CKjbva73a2LnMvrMzEwO95BdmKAQudhH3XpPDPpWXPj7Pbjw93tg0LdKGBVR37WVF0FoqjV7XRRFlJWVIS8vz4VRkbtjgkLkQo2t7ay2SR5HaL5kU7vKSv7bJ9vZlaCsWLECv/jFLxAUFIRhw4Zh1qxZKC4uNmnT2tqK9PR0DB06FIGBgZg9ezaqqqpM2pSWlmLmzJkICAjAsGHD8Nxzz6Gjo6P/d0MkI71V1Nx0WIPL7ezmJs+iDBxiU7vo6GgnR0KexK4EZd++fUhPT8f+/fuxc+dOtLe3Y+rUqSbljJ955hl8+eWX2LhxI/bt24eKigqkpqYarwuCgJkzZ0Kv1yM/Px/vvfceNmzYgBdeeMFxd0UkUx9dNTlWNFxJVlrLTpgcE7kLv9ixUAaFm72uUCgQFxeHlJQUF0ZF7k4hds1g6oOamhoMGzYM+/btw6RJk9DQ0ICIiAh89NFHUKlUAIDTp0/jhhtuQEFBAW677TZs27YN99xzDyoqKhAZGQkAWLNmDRYuXIiamhr4+vpa/bqNjY0ICQlBQ0MDgoOD+xo+kVPpdDoEBgYCAJqbm3G6tg2z3ykwXm8pzkfdrrUmqx6UQeEImzwPAQnJLo+XqD+6VvH0pIBCAWRnZ5v8sUreyZ7f3/2ag9LQ0AAACAsLAwAUFhaivb0dU6ZMMbYZPXo0RowYgYKCzh/MBQUFGD9+vDE5AYBp06ahsbERRUVF/QmHSNa6Ly02tyRTaKpFzeblaCnOd3V4RP0SkJCMiFmLoQwcanLeNyQcGzduZHJCdhvQ1zcaDAZkZmbil7/8JcaNGwcA0Gq18PX1RWhoqEnbyMhIaLVaY5vuyUnX9a5rvWlra0NbW5vxuLGxsa9hE7lM9yWV23buQc5RAFBYXJLZpW73Ogy6PhEKH6VTYyRypICEZPiNnIjyrPsBABGqpRg06iYMHXebxJGRO+pzD0p6ejpOnDiBTz75xJHx9GrFihUICQkxvuLi4pz+NYn6Q61WY8yYMcbj+379K/z01qNoKc63uiQT6OxJaStnjyK5n+5JtX/cOCh8lCxMSH3SpwRl/vz5yMnJwZ49exAbG2s8HxUVBb1ej/r6epP2VVVViIqKMra5elVP13FXm6stWrQIDQ0NxldZGTdZI/lSq9VQqVTQaEw3TTMO35w5YNPn2Lp0k0juvjldDW0D6/yQfexKUERRxPz587Fp0yZ88803GDVqlMn1m2++GQMHDsTu3buN54qLi1FaWoqkpCQAQFJSEo4fP47q6mpjm507dyI4ONjkL87u/Pz8EBwcbPIikiNBEJCRkQFLc891RXts+ixbl24SyZ1gEPHx9+xFIfvYNQclPT0dH330EbZs2YKgoCDjnJGQkBAMGjQIISEhSEtLw4IFCxAWFobg4GA89dRTSEpKwm23dY5BTp06FWPGjMHDDz+MlStXQqvV4vnnn0d6ejr8/Pwcf4dELpSXl4fy8nKLbQyXG+EzKBiGy+bnUimDwuEXO9bR4RE5nY+vP0YuzOlx/tODZXjqruswQMn6oGQbu/6lvPPOO2hoaMAdd9yB6Oho4+vTTz81tnn99ddxzz33YPbs2Zg0aRKioqJMNolSKpXIycmBUqlEUlISfvvb3+KRRx7BsmXLHHdXRBKxtVLm4LF3WrweNnkeJ8iSR9E2tmLXqc6e896KGBJdrV91UKTCOigkV3v37sWdd1pOPgAg8sHlMFxuZh0U8iop14fj32mJPWoEDR48WOLIyFXs+f3d52XGRNRTSkoKYmNjodFozM5D6Rq+Ufgoe12SyZ4T8lTfnq1F6cUWDPWXOhJyBxwMJHIgpVKJrKys/xwpem3TffimtyWZRJ5KFIEPv78gdRjkJpigEDlYamoqsrOzMSjUdG8SZVA4ImYt5vANebXsQ+VoadMbj3Nzc02KGhJ1YYJC5ARjkqdg6Ny3jMcRqqUY/sR6Jifk9coO70XC6CslJWbMmIH4+HiTxRREAOegEDnFhu/O2zR8Y25JJpEnMrehoEajgUql4oaCZII9KEQOdkmnx+ajGusNibyIpT2ouiaUZ2ZmcriHjNiDQuRgHx8sRWu7gb0jRN1Y24NKFEWUlZUhLy8Pd9xxh+sCI9liDwqRA3UIBvy7gKsUiK5m695SthY7JM/HBIXIgbYXaVHJTdGIerB1b6no6GgnR0LuggkKkQNt+O681CEQyZJf7Fgog8LNXlcoFIiLi0NKSooLoyI5Y4JC5CDHyxtw6IJt3dhE3kbho0TY5Hm9X1N0FjVctWoVlEoWK6ROTFCIHOTd70qkDoFI1gISkhExazGUgUNNzg+PjeUSY+qBCQqRA9Q0tSHnGCf3EVkTkJCM6LS3jccRqqV4+4t8JifUAxMUIgf48MAF6AWD1GEQuYWrixi+V1AqYTQkV0xQiPpJ32HAhwf4A5aorw5duIRj5fVSh0EywwSFqJ++Ol6BmqY2qcMgchtdRQxHLsyBj68/AOBf33IOF5ligkLUD6IoYu2+n6QOg8jtfXW8ElWNrCFEVzBBIeqHb05X47S2SeowiNxeuyDi/YLzUodBMsIEhagPdDodFAoFpoyJgkHPv/qIHOGjA6W4rOdmgdSJCQoREcnCpZZ2fPQ9J5xTJyYoREQkG//I/Qn6Di7ZJyYoRH1yovxKSfvWshMQDeyWJnIEbWMrsgvLpQ6DZIAJCpGd1Go1Jif93Hhck70UmjVpaCnOlzAqIs+xZt85CAZR6jBIYkxQiOygVquhUqmgu1Rjcl5oqkXN5uVMUogcoLSuBV/8oJE6DJIYExQiGwmCgIyMDIii+b/s6nav43APkQO8uf0EFAoFFAoFdDqd1OGQBJigENkoLy8P5eWWx8aFplq0lRe5KCIiz3W2hkmJt2OCQmSjykrbdisWmi9Zb0RERBYxQSGyUdCQCJvaKQOHODkSIs/Xfag0NzcXgsChU2/DBIXIRud8YqEMCrfYRhkUDr/YsS6KiMgztRTno3L9H4zHM2bMQHx8PNRqtYRRkasxQSGyQWu7gPcPlCFs8jyL7cImz4PCR+miqIg8T0txPmo2L4fQfNHkvEajgUqlYpLiRZigENng4+9LcVGnR0BCMiJmLYYycKjJdWVQOCJmLUZAQrJEERK5P9EgoG73ut6v/Wf1XGZmJod7vAQTFCIrLusFvL33nPE4ICEZ0WlvG48jVEsx/In1TE6I+qmtvAhCU63Z66IooqysDHl5eS6MiqTCBIXIig3551HT1GZyrvswjn/cOA7rEDmArSvgbF1RR+5tgNQBEMlZY2s71uae63Hex9cfIxfmSBARkeeydQVcdHS0kyMhOWAPCpEF/8j9CfUt7VKHQeQV/GLHWl0pFxcXh5SUFBdFRFJigkJkxsXmNvzr2xKpwyDyGgofpdWVcn9Y9BKUSg6pegMmKERmrN5zDjo9VwsQuZK1lXIHxOskioxcjXNQiHpRUX8ZHxy4IHUYRF4pICEZfiMnojzrfgCdK+UGjboJCh8ljpTWY0eRFtPGRkkcJTkbe1CIrqLT6TB8SADO/G0GDPpWqcMh8kqWVsq9uqMYgsH8ruLkGdiDQnSVC7XcRZVIapZWyp2pbsbnh8vxm1viXBwVuRJ7UIiukrWr2PjfrWUnTDYtIyJ5yNp1Bm0d/N70ZExQiLp5Y/0HWJ/5a+NxTfZSaNakoaU4X8KoiOhqmvrL+HcB54l5MiYoRP+hVquR8fuHe2xSJjTVombzciYpRDKzes9ZNLWyTpGnYoJCBEAQBDw5/ymLbep2r+NwD5GMXKxvQvAgXygUCuh0nDvmaZigEAHIy8tDdWWFxTZCUy3ayotcFBERkXdjgkIE4MsC2xIPWzczIyLn696jmZubC0FgD6cnYYJCXk/X1oGt52yrd2LrZmZE5FwtxfmoXP8H4/GMGTMQHx8PtVotYVTkSExQyOu9tecsWsKut7pJmTIoHH6xY10UFRGZ01Kcj5rNy3tMaNdoNFCpVExSPAQTFPJq52t1WJ9XYtMmZWGT55lUsyQi1xMNAup2r+v9mthZXTYzM5PDPR6ACQp5tWU5J6EXDACsb1IWkJAsRYhE1E1beRGEplqz10VRRFlZGfLy8lwYFTkDS92T19pzuhrfnK42OWdpkzIikp6tE9UrKyudHAk5m909KLm5ufjVr36FmJgYKBQKbN682eS6KIp44YUXEB0djUGDBmHKlCk4c+aMSZu6ujrMmTMHwcHBCA0NRVpaGpqbm/t1I0T20HcYsCznZK/XlP6DMXJhDkYuzEHAtbcwOSGSEVsnqkdHRzs5EnI2uxMUnU6HiRMnYvXq1b1eX7lyJd544w2sWbMGBw4cwODBgzFt2jS0tl5ZJTFnzhwUFRVh586dyMnJQW5uLubNszz+T+RI678tQQk3BSRyO36xY61OaA8YEomUlBQXRUTOohC7ZhX15c0KBTZt2oRZs2YB6Ow9iYmJwbPPPos//vGPAICGhgZERkZiw4YNeOCBB3Dq1CmMGTMGBw8exC233AIA2L59O2bMmIHy8nLExMRY/bqNjY0ICQlBQ0MDgoOD+xo+eamqxlbc9epe6PScREfkjrpW8ZgTMWsx3ls2H9PHsxdFbuz5/e3QSbIlJSXQarWYMmWK8VxISAgSExNRUFAAACgoKEBoaKgxOQGAKVOmwMfHBwcOHOj1c9va2tDY2GjyIuqrFVtPMTkhcmO2TGh/KeckLusF6HQ6KBQKlsN3Qw5NULRaLQAgMjLS5HxkZKTxmlarxbBhw0yuDxgwAGFhYcY2V1uxYgVCQkKMr7i4OEeGTR7Clh9Eh87XYfNRyyXtiUj+AhKSEZ32tvE4QrUUw59Yb1xtV9HQije/OWOy3JjVZt2LWywzXrRoERoaGoyvsrIyqUMiGbL2g0jfYcBfNp1wdVhE5CTWJrSv+scH+FnCaOMxq826F4cmKFFRUQCAqqoqk/NVVVXGa1FRUaiuNl3a2dHRgbq6OmObq/n5+SE4ONjkRdSdWq3GmDFjjMe9/SB6Y/cZFFc1SREeEblYS3E+KtV/Q5XWdLkxq826D4cmKKNGjUJUVBR2795tPNfY2IgDBw4gKSkJAJCUlIT6+noUFhYa23zzzTcwGAxITEx0ZDjkJdRqNVQqFTQajcn57j+IjpXXY82+cxJFSESuxGqznsHuQm3Nzc04e/as8bikpARHjx5FWFgYRowYgczMTPzv//4vrr/+eowaNQp//etfERMTY1zpc8MNN+Duu+/GY489hjVr1qC9vR3z58/HAw88YNMKHqLuBEFARkYGeluMJooiFAoFMjIyMDrzPXQYRAitOhZhI/Jw9lSbveOOO1wXGNnF7gTl0KFDuPPOO43HCxYsAADMnTsXGzZswJ/+9CfodDrMmzcP9fX1+K//+i9s374d/v7+xvd8+OGHmD9/PiZPngwfHx/Mnj0bb7zxhgNuh7xNXl4eysvLzV4XRRHl5eVoLzwAw+Vm1O1aa7xWk70UyqBwhE2exzL2RB6E1WY9Q7/qoEiFdVCoy8cff4yHHnrIarugW+5F06EtZq9zrx0iz9FaegxVHy+22m7Pnj3sQXExyeqgELmareWsdUV7LF6v270OooHj0USewFq1WYVCgbi4OFablTkmKOTWUlJSEBsbC4VCYbaNz6BgGC5bLu4nNNWirbzI0eERkQQUPkqETba8fcqqVaugVHL+mZwxQSG3plQqkZWVBQBmk5TBY+/s9fzVbB23JiL5s1RtNm3Jm0hNTZUoMrIVExRye6mpqcjOzu6xCsxY9vp625av27pLKhG5B3PVZvfqR6GookHCyMgWnCRLHkMQBOTl5eGdrYewr7wdfrFjofBRQjQI0KxJs7jsUBkUjuFPrOeSYyIvMToqCF/M/y/4DuDf6a7ESbLklZRKJTqG3YDvlTfAf8QEY7Jhy3h02OR5TE6IvMhpbRP+b2ex1GGQBUxQyGOUXmzBgs+Oorc+QVt2PyUi7/KP3J+Qf7YWjY2Nxo1Gt23bxgqzMsEEhTxCa7uAJz4oRGNrh9k21nY/JSLvYhCB3/71TSSMvsF4jhsKygfnoJBHeG7jD9hYaL6iLBHR1VqK81GzeXmP810rArOzs7nax8E4B4U8jqUu2E++L2VyQkR24YaC8scEhWRPrVZjzJgxxuPuXbAnNA1Y8gULrBGRfezZUJCkYfdmgUSupFaroVKpeuxWrNFooFKpkPDbJWiLuUWi6IjIXXFDQfljDwrJliAIyMjI6JGcAJ1/3YiiiDOb3+IeOkRkN1sLM9q63xc5HhMUkq28vDyUl1ueW8I9dIioL7ihoPwxQSHZsrVrlXvoEJG9rBVwFEXTDQVZK8X1mKCQbNnatco9dIioLywXcFyEhsibAFieqE/OwzooJFuCICA+Ph4ajabXeSgA99Ahov4TWnUoz7ofQGcBx0GjboLCRwmFAvjN0Aq88qfHe/wMYq2UvrHn9zcTFJK1rlU8AHpNUlimnoicxdpGowqFArGxsSgpKTEOBZFlLNRGHiM1NRXr3/8IA4O4hw4RuRZrpUiLdVBI1ppa2/HZxeGIenx95w+L5ktQBg6BX+xYDusQkVOxVoq0mKCQbLV1CHjs/UM4VdkIhY8S/iMmSB0SEXkR1kqRFod4SJbaBQMyPj6K/T/VSR0KEXkp1kqRFhMUkpwgCNi7dy8+/vhj7N27F7pWPea9fwjbi7RSh0ZEXsxSrZSuVTzda6WQY3GIhySlVqvx9NNPQ6PRGM8NCh2GwDt+zwmwRCS5rlopdbvWQmi+aDzvGxyOrKwsLjF2Ii4zJsmY2wiwC1fpEJFciAahx0T9iOAAvPe7X2BsTIjU4bkN1kEh2esqwmZprx0WYSMiuQvyG4C1j9yM5GvNz1WhK1gHhWSPGwESkSdoauvA//zrIL46xqXGjsY5KCQJbgRIRJ5CLxjw1MeHUd0wGqOEMlRWViI6OhopKSmcQNsPTFBIEtwIkIg8SfPpfDz21lyTibSxsbGcSNsPHOIhp7p6CbEgCBAMIgp0ERbrCwCdc1D8Yse6KFIior5pKc5HzeblJskJAGg0GqhUKu563EfsQaE+EwQBeXl5Zrsze1tCHDN8OK751XyUhYxH2OR5qNm83Oznh02exwmyRCRrokFA3e51vV8TRSgUCmRmZuLee+/lcI+d2INCfaJWqzFy5EjceeedeOihh3DnnXciPj7e+JdC1xLi7skJAFRoNPh2zSK0FOcb6wtc3ZPCjQCJyF3Yu6Fgb73K1Dv2oJDdzNUv6erO/PTTT7FgwQKz9U0AoG73Ogy6PhEBCckYdH0iNwIkIrdk60T+82XlUKvVyMjIMFnByHkq5rEOCtnFWv0ShUKB8PBw1NTUWP2syAeXcwNAInJrraXHUPXxYqvtRk17FOe/3tDjD7eukvnZ2dlekaSwDgo5jbX6JaIo2pScAFxCTETuz9qGggDgEzgUpd990Wuvcte5zMxM43APh4E6MUEhu9hav8QWXEJMRO7O0oaCXYIm3t1jhU933eepqNVqxMfHm53f502YoJBdbK1f4jPIctcdlxATkaewNuF/YFiMTZ+zZcsWqFSqHr3U3rpcmXNQyC5dc1A0Go3ZSbDKoHAMuev3qN3ystnP4SodIvI0vW0oqPBR2jxPJSIiwuwQuUKhQGxsLEpKSqBUKq2WeZAre35/cxUP2UWpVCIrKwsqlQoKhaLXJCVs8jwEJCRDoViMut3rTJbgKYPCjdeJiDyJwkfZ68T/rnkqlpYj+wcNsTh/r/swUF1dnVesBmIPCvXKWna+6h//xuI/PYvL9Ve+oXpLPsz9RUFE5E26qs2aE3TLvWg6tMXq52RmZiIrK8ttVwPZ8/ubCQr1YGmt/qSpM/HK9mJsLCyDIDD5ICKyVUtxvtleZZ9BgV4xDMQEhfrMXBG2zuEcIO43f4HPqNskio6IyL2Z61UWDQI0a9LMDgPZU2Nqz549sh0GYh0UssjcGntBEJCRkWFhrb6Iim1rIBq8c00+EVF/dc1TGTzmdviPmGDsdba2XFkUgbvumW3T17BnNZCca66wB8WNOKK7ztLwTVhYGO68806rn8EKsEREzuHKYaAtW7a4vJeFq3g8kCP2cLC2h87MB9Ns+hxWgCUicg5L+5OJBsHiaiBbhoG6VgP97W9/w9KlS83+PpDDZFv2oLgBS/NCANtmbVvbQwfoLK5muNxoNR72oBARScPiaiCFAvOfegpvvfGG1c8JCwtDXV2dmY8xnWzrSJyDIjO2jPH1fV6I6R4O5ljbQwcADJcbWQGWiEjGLFatvXcRvrg03KbPMZecAKY1V6TEIR4ns2Voxtq8EGub83X9Q7rjjjt6nadS3axHzv6TNsU7eOydFtfih02ex6XEREQS6u8w0JAhQywmKF0cufdaXzBBsaC/k1KtzfnIzs4GAIttMjIybPpalZWVvSY6viERCLnzMfgMCrTpcwKuT4R/7FhWgCUikjFzVWu7VgOZGwYSRRG/nPUwvvxXltWvYevea87COShm2Dop1VwSY23Oh0KhwPDhnV1xltrYuu797rlPY/t75scdw+/9My5980+LpZaVQeEY/sR6YxbOImxERO7J0mqgQdcnWq25Ioc5KJImKKtXr8Yrr7wCrVaLiRMn4s0338Stt95q9X3OTlBsnZTqiCW7toiIiEBtba3Zzfl8AodCAVjczpsb+BEReRdLf2haK73/4hv/wgtPPerwmNxikuynn36KBQsWYMmSJTh8+DAmTpyIadOmobq6WqqQANg+KTU7O9tiIZwtW6zvqWCz6/7LbHICAEET77aYnACA0FQLZUCwxS3BmZwQEXkOc0XhACuTbWcths+oRFeH24NkPSiJiYn4xS9+gbfeegsAYDAYEBcXh6eeegp//vOfLb7XmT0oe/futannw1ohHFuHZmwR+eByGC43m+2uE4V21H75itXPCf/Vcxg85nYO3xAREQDzvSy//69ReP6eMQ7/erIv1KbX61FYWIhFixYZz/n4+GDKlCkoKCjo0b6trQ1tbW3G48ZG67U6+upCmcamdtYK4dTU1FitK2Lr0EzXPxhzs7ZbS4/ZFLMycAgA85OriIjIu8j594EkCUptbS0EQUBkZKTJ+cjISJw+fbpH+xUrVuDFF190SWw+g4c47LNun5GKPZ+/95+j7h1VnXNZfvfHznv659KnzLZ59NmluHFS9yy2Z0ZrEK7DC7veQH1N1VWfceWzQiOi8JffpcJH4p0siYhI/m6MC5U6BPdYZrxo0SIsWLDAeNzY2Ii4uDinfK2H7p2GxbGx0Gg0vc77sGf45oX5czH/oV/1mEgbFxeLVatWGVcDTR8fbbWNNcPXrIZKpQKgMIm7a2Lv+jVvIXVKgk2fRUREJDVJ5qDo9XoEBAQgOzsbs2bNMp6fO3cu6uvrrU4wddUqHgC9/rLvmuBrKYnpvkTLlnoqztoIMC4uzq5Eh4iIyFncYplxYmIibr31Vrz55psAOifJjhgxAvPnz5d0kmwXa7/srSUxUm205IhEh4iIyBncIkH59NNPMXfuXKxduxa33norVq1ahc8++wynT5/uMTflaq7aLNDaL3v2WBAREdnOLRIUAHjrrbeMhdpuvPFGvPHGG0hMtL72Wk67GbPHgoiIyDZuk6D0lZwSFCIiIrKNW1SSJSIiIjKHCQoRERHJDhMUIiIikh0mKERERCQ7TFCIiIhIdpigEBERkewwQSEiIiLZYYJCREREssMEhYiIiGRngNQB9EVX8dvGxkaJIyEiIiJbdf3etqWIvVsmKE1NTQA6N+YjIiIi99LU1ISQkBCLbdxyLx6DwYCKigoEBQVBoVD06TMaGxsRFxeHsrIyj93Ph/foGbzhHgHvuE/eo2fwhnsEnHOfoiiiqakJMTEx8PGxPMvELXtQfHx8EBsb65DPCg4O9uh/YADv0VN4wz0C3nGfvEfP4A33CDj+Pq31nHThJFkiIiKSHSYoREREJDtem6D4+flhyZIl8PPzkzoUp+E9egZvuEfAO+6T9+gZvOEeAenv0y0nyRIREZFn89oeFCIiIpIvJihEREQkO0xQiIiISHaYoBAREZHseGyC8re//Q3JyckICAhAaGhor21KS0sxc+ZMBAQEYNiwYXjuuefQ0dFh8XPr6uowZ84cBAcHIzQ0FGlpaWhubnbCHdhv7969UCgUvb4OHjxo9n133HFHj/ZPPPGECyO3T3x8fI94X375ZYvvaW1tRXp6OoYOHYrAwEDMnj0bVVVVLorYPufPn0daWhpGjRqFQYMG4dprr8WSJUug1+stvs8dnuPq1asRHx8Pf39/JCYm4vvvv7fYfuPGjRg9ejT8/f0xfvx4bN261UWR2m/FihX4xS9+gaCgIAwbNgyzZs1CcXGxxfds2LChxzPz9/d3UcT2W7p0aY94R48ebfE97vQMu/T2M0ahUCA9Pb3X9u7wHHNzc/GrX/0KMTExUCgU2Lx5s8l1URTxwgsvIDo6GoMGDcKUKVNw5swZq59r7/e0PTw2QdHr9bjvvvvw5JNP9npdEATMnDkTer0e+fn5eO+997Bhwwa88MILFj93zpw5KCoqws6dO5GTk4Pc3FzMmzfPGbdgt+TkZFRWVpq8fv/732PUqFG45ZZbLL73scceM3nfypUrXRR13yxbtswk3qeeespi+2eeeQZffvklNm7ciH379qGiogKpqakuitY+p0+fhsFgwNq1a1FUVITXX38da9asweLFi62+V87P8dNPP8WCBQuwZMkSHD58GBMnTsS0adNQXV3da/v8/Hw8+OCDSEtLw5EjRzBr1izMmjULJ06ccHHkttm3bx/S09Oxf/9+7Ny5E+3t7Zg6dSp0Op3F9wUHB5s8swsXLrgo4r4ZO3asSbzffvut2bbu9gy7HDx40OQed+7cCQC47777zL5H7s9Rp9Nh4sSJWL16da/XV65ciTfeeANr1qzBgQMHMHjwYEybNg2tra1mP9Pe72m7iR7u3XffFUNCQnqc37p1q+jj4yNqtVrjuXfeeUcMDg4W29raev2skydPigDEgwcPGs9t27ZNVCgUokajcXjs/aXX68WIiAhx2bJlFtvdfvvtYkZGhmuCcoCRI0eKr7/+us3t6+vrxYEDB4obN240njt16pQIQCwoKHBChI63cuVKcdSoURbbyP053nrrrWJ6errxWBAEMSYmRlyxYkWv7X/zm9+IM2fONDmXmJgoPv74406N01Gqq6tFAOK+ffvMtjH380mulixZIk6cONHm9u7+DLtkZGSI1157rWgwGHq97m7PEYC4adMm47HBYBCjoqLEV155xXiuvr5e9PPzEz/++GOzn2Pv97S9PLYHxZqCggKMHz8ekZGRxnPTpk1DY2MjioqKzL4nNDTUpDdiypQp8PHxwYEDB5wes72++OILXLx4EY8++qjVth9++CHCw8Mxbtw4LFq0CC0tLS6IsO9efvllDB06FDfddBNeeeUVi0NzhYWFaG9vx5QpU4znRo8ejREjRqCgoMAV4fZbQ0MDwsLCrLaT63PU6/UoLCw0eQY+Pj6YMmWK2WdQUFBg0h7o/B51p2cGwOpza25uxsiRIxEXF4d7773X7M8fuThz5gxiYmJwzTXXYM6cOSgtLTXb1t2fIdD5b/eDDz7A7373O4ub07rbc+yupKQEWq3W5FmFhIQgMTHR7LPqy/e0vdxys0BH0Gq1JskJAOOxVqs1+55hw4aZnBswYADCwsLMvkdK69evx7Rp06xurPjQQw9h5MiRiImJwbFjx7Bw4UIUFxdDrVa7KFL7PP300/j5z3+OsLAw5OfnY9GiRaisrMRrr73Wa3utVgtfX98ec5EiIyNl+dyudvbsWbz55pt49dVXLbaT83Osra2FIAi9fs+dPn261/eY+x51h2dmMBiQmZmJX/7ylxg3bpzZdgkJCfjXv/6FCRMmoKGhAa+++iqSk5NRVFTksA1RHSkxMREbNmxAQkICKisr8eKLLyIlJQUnTpxAUFBQj/bu/Ay7bN68GfX19fif//kfs23c7Tleret52POs+vI9bS+3SlD+/Oc/4+9//7vFNqdOnbI6acvd9OW+y8vLsWPHDnz22WdWP7/7HJrx48cjOjoakydPxrlz53Dttdf2PXA72HOPCxYsMJ6bMGECfH198fjjj2PFihWyLj3dl+eo0Whw991347777sNjjz1m8b1yeI7UKT09HSdOnLA4PwMAkpKSkJSUZDxOTk7GDTfcgLVr1+Kll15ydph2mz59uvG/J0yYgMTERIwcORKfffYZ0tLSJIzMedavX4/p06cjJibGbBt3e47uwq0SlGeffdZiFgsA11xzjU2fFRUV1WO2cdeqjqioKLPvuXryT0dHB+rq6sy+xxH6ct/vvvsuhg4div/3//6f3V8vMTERQOdf7q76xdafZ5uYmIiOjg6cP38eCQkJPa5HRUVBr9ejvr7epBelqqrKqc/tavbeY0VFBe68804kJydj3bp1dn89KZ6jOeHh4VAqlT1WTll6BlFRUXa1l4v58+cbJ9Db+9fzwIEDcdNNN+Hs2bNOis6xQkND8bOf/cxsvO76DLtcuHABu3btsrsX0t2eY9fzqKqqQnR0tPF8VVUVbrzxxl7f05fvabs5ZCaLjFmbJFtVVWU8t3btWjE4OFhsbW3t9bO6JskeOnTIeG7Hjh2ymyRrMBjEUaNGic8++2yf3v/tt9+KAMQffvjBwZE5xwcffCD6+PiIdXV1vV7vmiSbnZ1tPHf69GlZT5ItLy8Xr7/+evGBBx4QOzo6+vQZcnuOt956qzh//nzjsSAI4vDhwy1Okr3nnntMziUlJcl2gqXBYBDT09PFmJgY8ccff+zTZ3R0dIgJCQniM8884+DonKOpqUkcMmSImJWV1et1d3uGV1uyZIkYFRUltre32/U+uT9HmJkk++qrrxrPNTQ02DRJ1p7vabvjdMinyNCFCxfEI0eOiC+++KIYGBgoHjlyRDxy5IjY1NQkimLnP6Bx48aJU6dOFY8ePSpu375djIiIEBctWmT8jAMHDogJCQlieXm58dzdd98t3nTTTeKBAwfEb7/9Vrz++uvFBx980OX3Z8muXbtEAOKpU6d6XCsvLxcTEhLEAwcOiKIoimfPnhWXLVsmHjp0SCwpKRG3bNkiXnPNNeKkSZNcHbZN8vPzxddff108evSoeO7cOfGDDz4QIyIixEceecTY5up7FEVRfOKJJ8QRI0aI33zzjXjo0CExKSlJTEpKkuIWrCovLxevu+46cfLkyWJ5eblYWVlpfHVv427P8ZNPPhH9/PzEDRs2iCdPnhTnzZsnhoaGGlfSPfzww+Kf//xnY/vvvvtOHDBggPjqq6+Kp06dEpcsWSIOHDhQPH78uFS3YNGTTz4phoSEiHv37jV5Zi0tLcY2V9/jiy++KO7YsUM8d+6cWFhYKD7wwAOiv7+/WFRUJMUtWPXss8+Ke/fuFUtKSsTvvvtOnDJlihgeHi5WV1eLouj+z7A7QRDEESNGiAsXLuxxzR2fY1NTk/H3IADxtddeE48cOSJeuHBBFEVRfPnll8XQ0FBxy5Yt4rFjx8R7771XHDVqlHj58mXjZ9x1113im2++aTy29j3dXx6boMydO1cE0OO1Z88eY5vz58+L06dPFwcNGiSGh4eLzz77rEmmvGfPHhGAWFJSYjx38eJF8cEHHxQDAwPF4OBg8dFHHzUmPXLx4IMPisnJyb1eKykpMfn/UFpaKk6aNEkMCwsT/fz8xOuuu0587rnnxIaGBhdGbLvCwkIxMTFRDAkJEf39/cUbbrhBXL58uUmv19X3KIqiePnyZfEPf/iDOGTIEDEgIED89a9/bfILX07efffdXv/tdu/wdNfn+Oabb4ojRowQfX19xVtvvVXcv3+/8drtt98uzp0716T9Z599Jv7sZz8TfX19xbFjx4pfffWViyO2nbln9u677xrbXH2PmZmZxv8fkZGR4owZM8TDhw+7Pngb3X///WJ0dLTo6+srDh8+XLz//vvFs2fPGq+7+zPsbseOHSIAsbi4uMc1d3yOXb/Prn513YfBYBD/+te/ipGRkaKfn584efLkHvc+cuRIccmSJSbnLH1P95dCFEXRMYNFRERERI7htXVQiIiISL6YoBAREZHsMEEhIiIi2WGCQkRERLLDBIWIiIhkhwkKERERyQ4TFCIiIpIdJihEREQkO0xQiIiISHaYoBAREZHsMEEhIiIi2WGCQkRERLLz/wFAUp2dHJH56wAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = data.to_numpy()[\"x\"]\n", "\n", "nll = UnbinnedNLL(x, model)\n", "m = Minuit(nll, 1, 3)\n", "m.migrad()" ] } ], "metadata": { "kernelspec": { "display_name": "root", "language": "python", "name": "root" }, "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.11.0" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }