{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e0a3008f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iminuit version: 2.16.0\n"
     ]
    }
   ],
   "source": [
    "# Example of maximum-likelihood fit with iminuit version 2.x\n",
    "# pdf is a mixture of Gaussian (signal) and exponential (background),\n",
    "# truncated in [xMin,xMax].\n",
    "# G. Cowan / RHUL Physics / December 2022\n",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "from scipy.stats import truncexpon\n",
    "from scipy.stats import truncnorm\n",
    "from scipy.stats import chi2\n",
    "import iminuit\n",
    "from iminuit import Minuit\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"font.size\"] = 14\n",
    "print(f\"iminuit version: {iminuit.__version__}\")      #  should be v 2.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9da621fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select fit type, define pdf\n",
    "fitType = 'M'        # choose least squares (LS) or multinomial ML (M)\n",
    "np.random.seed(seed=1234567)        # fix random seed\n",
    "theta = 0.2                         # fraction of signal\n",
    "mu = 10.                            # mean of Gaussian\n",
    "sigma = 2.                          # std. dev. of Gaussian\n",
    "xi = 5.                             # mean of exponential\n",
    "xMin = 0.\n",
    "xMax = 20.\n",
    "\n",
    "def f(x, par):                      # fit function\n",
    "    theta   = par[0]\n",
    "    mu      = par[1]\n",
    "    sigma   = par[2]\n",
    "    xi      = par[3]\n",
    "    fs = stats.truncnorm.pdf(x, a=(xMin-mu)/sigma, b=(xMax-mu)/sigma, loc=mu, scale=sigma)\n",
    "    fb = stats.truncexpon.pdf(x, b=(xMax-xMin)/xi, loc=xMin, scale=xi)\n",
    "    return theta*fs + (1-theta)*fb\n",
    "\n",
    "class ChiSquared:                   # function to be minimized\n",
    "\n",
    "    def __init__(self, xHist, bin_edges, fitType):\n",
    "        self.setData(xHist, bin_edges)\n",
    "        self.fitType = fitType\n",
    "        \n",
    "    def setData(self, xHist, bin_edges):\n",
    "        numVal = np.sum(xHist)\n",
    "        numBins = len(xHist)\n",
    "        binSize = bin_edges[1] - bin_edges[0]\n",
    "        self.data = xHist, bin_edges, numVal, numBins, binSize\n",
    "\n",
    "    def chi2LS(self, par):           # least squares\n",
    "        xHist, bin_edges, numVal, numBins, binSize = self.data\n",
    "        xMid = bin_edges[:numBins] + 0.5*binSize\n",
    "        binProb = f(xMid, par)*binSize\n",
    "        nu = numVal*binProb\n",
    "        sigma = np.sqrt(nu)\n",
    "        z = (xHist - nu)/sigma\n",
    "        return np.sum(z**2)\n",
    "        \n",
    "    def chi2M(self, par):            # multinomial maximum likelihood\n",
    "        xHist, bin_edges, numVal, numBins, binSize = self.data\n",
    "        xMid = bin_edges[:numBins] + 0.5*binSize\n",
    "        binProb = f(xMid, par)*binSize\n",
    "        nu = numVal*binProb\n",
    "        lnL = 0.\n",
    "        for i in range(len(xHist)):\n",
    "            if xHist[i] > 0.:\n",
    "                lnL += xHist[i]*np.log(nu[i]/xHist[i])\n",
    "        return -2.*lnL\n",
    "\n",
    "    def __call__(self, par):\n",
    "        if self.fitType == 'LS':\n",
    "            return self.chi2LS(par)\n",
    "        elif self.fitType == 'M':\n",
    "            return self.chi2M(par)\n",
    "        else:\n",
    "            print(\"fitType not defined\")\n",
    "            return -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "da616747",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate data\n",
    "numVal = 200\n",
    "xData = np.empty([numVal])\n",
    "for i in range (numVal):\n",
    "    r = np.random.uniform();\n",
    "    if r < theta:\n",
    "        xData[i] = stats.truncnorm.rvs(a=(xMin-mu)/sigma, b=(xMax-mu)/sigma, loc=mu, scale=sigma)\n",
    "    else:\n",
    "        xData[i] = stats.truncexpon.rvs(b=(xMax-xMin)/xi, loc=xMin, scale=xi)\n",
    "\n",
    "# Put data values into a histogram\n",
    "numBins=40\n",
    "xHist, bin_edges = np.histogram(xData, bins=numBins, range=(xMin, xMax))\n",
    "binSize = bin_edges[1] - bin_edges[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "862e05d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize Minuit and set up fit:\n",
    "parin   = np.array([theta, mu, sigma, xi])      # initial values (here = true)\n",
    "parname = ['theta', 'mu', 'sigma', 'xi']\n",
    "parname_latex = [r'$\\theta$', r'$\\mu$', r'$\\sigma$', r'$\\xi$']\n",
    "parstep = np.array([0.1, 1., 1., 1.])           # initial setp sizes\n",
    "parfix  = [False, True, True, False]            # change to fix/free param.\n",
    "parlim  = [(0.,1), (None, None), (0., None), (0., None)]\n",
    "chisq = ChiSquared(xHist, bin_edges, fitType)\n",
    "m = Minuit(chisq, parin, name=parname)\n",
    "m.errors = parstep\n",
    "m.fixed = parfix\n",
    "m.limits = parlim\n",
    "m.errordef = 1.0                      # errors from chi2 = chi2min + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7fa0f30e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# do the fit, get errors, extract results\n",
    "m.migrad()                                        # minimize -logL\n",
    "parhat = m.values                                 # max-likelihood estimates\n",
    "sigma_parhat = m.errors                           # standard deviations\n",
    "cov = m.covariance                                # covariance matrix\n",
    "rho = m.covariance.correlation()                  # correlation coeffs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4f0ada34",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "par index, name, estimate, standard deviation:\n",
      "   0 theta       =  0.201533  +/-  0.052956\n",
      "   3 xi          =  5.148508  +/-  0.651156\n",
      "\n",
      "free par indices, covariance, correlation coeff.:\n",
      "0 0 0.002821 1.000000\n",
      "0 3 -0.018494 -0.534693\n",
      "3 0 -0.018494 -0.534693\n",
      "3 3 0.424135 1.000000\n",
      "chi2min =  35.26955011049647 ,     ndof =  37\n",
      "pval =  0.5503255503742259\n"
     ]
    }
   ],
   "source": [
    "print(r\"par index, name, estimate, standard deviation:\")\n",
    "for i in range(m.npar):\n",
    "    if not m.fixed[i]:\n",
    "        print(\"{:4d}\".format(i), \"{:<10s}\".format(m.parameters[i]), \" = \",\n",
    "         \"{:.6f}\".format(parhat[i]), \" +/- \", \"{:.6f}\".format(sigma_parhat[i]))\n",
    "         \n",
    "print()\n",
    "print(r\"free par indices, covariance, correlation coeff.:\")\n",
    "for i in range(m.npar):\n",
    "    if not m.fixed[i]:\n",
    "        for j in range(m.npar):\n",
    "            if not m.fixed[j]:\n",
    "                print(i, j, \"{:.6f}\".format(cov[i,j]), \"{:.6f}\".format(rho[i,j]))\n",
    "\n",
    "# Retrieve minimized chi-squared, etc.\n",
    "ndof = numBins - m.nfit - 1           # for fixed numVal\n",
    "chi2min = chisq(parhat)\n",
    "print (\"chi2min = \", chi2min, \",     ndof = \", ndof)\n",
    "pval = chi2.sf(chi2min, ndof)\n",
    "print (\"pval = \", pval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b3e39fab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD/CAYAAADxL6FlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABKGklEQVR4nO3dd1yV5fvA8c8FiCC4BUVQgTQHgqgoOdNy61dzpZZlWmml7aGV7aGWVjZ+pVmZVqbZUMsstRyZIwfuLAcpYm4FAZn374/nYKiMwzyM6/163Z1znnmdE57rPPdzDzHGoJRSSjk5OgCllFLFgyYEpZRSgCYEpZRSNpoQlFJKAZoQlFJK2WhCUEopBYCLowMoSDVq1DD+/v6ODkMppRxuy5Ytp4wxXrnZp1QlBH9/fzZv3uzoMJRSyuFE5J/c7qNVRkoppQBNCEoppWw0ISillAIKMSGISB0R+VVE9orIbhF50La8mogsF5G/bY9Vs9i/h4jsE5H9IjKhsOJUSillKcwrhBTgUWNMY+A6YKyINAEmACuNMQ2AlbbXlxERZ+A9oCfQBBhm21cppVQhKbSEYIw5ZozZanseC+wFfIF+wKe2zT4Fbspk99bAfmPMQWNMEvClbT+llFKFpEjuIYiIP9Ac2AjUNMYcAytpAN6Z7OILHMnwOsq2TCmlVCEp9IQgIp7A18BDxpgYe3fLZFmmEzeIyGgR2Swim0+ePJnXMJVSqswr1IQgIuWwksHnxphvbIuPi4iPbb0PcCKTXaOAOhle+wHRmZ3DGDPTGBNmjAnz8spVpzyllFIZFGYrIwE+AvYaY97IsGoxMML2fASwKJPd/wAaiEiAiLgCQ237KaWUKiSFeYXQDrgNuEFEImylFzAZ6CoifwNdba8RkdoishTAGJMCjAN+wroZvcAYs7sQY1VKqTKv0MYyMsb8Rub3AgBuzGT7aKBXhtdLgaWFE51SSqkraU9lpZRSgCYEpRxKRLjtttsuvU5JScHLy4s+ffrk6XiLFy9m8uTJBRVernXq1OmqEYc3b97MAw88AMDzzz/P1KlT7T6ev78/p06dAqBt27YArFq1Ks+fT3YK67glSaka/lqpksbDw4Ndu3aRkJCAu7s7y5cvx9c3711u+vbtS9++fQswwvwLCwsjLCws38f5/fffCyAalR29QlDKwXr27MkPP/wAwLx58xg2bNildZs2baJt27Y0b96ctm3bsm/fPgDeeOMNRo0aBcDOnTtp2rQp8fHxzJ49m3HjxgFwxx13cO+999K5c2cCAwNZvXo1o0aNonHjxtxxxx2XzuHp6Xnp+cKFCy+ts3f/nGT1y/vDDz+kZ8+eJCQk8Nlnn9G6dWtCQ0MZM2YMqampV22fMc4LFy4waNAgGjVqxK233ooxVjellStX0rx5c4KDgxk1ahSJiYnZLl+2bBmNGjWiffv2fPPNN1eds6wpVQnh4Mk4hsxYn2X5YuNhR4eo1FWGDh3Kl19+ycWLF9mxYwfh4eGX1jVq1Ig1a9awbds2XnzxRZ566ikAHnroIfbv38+3337LyJEjmTFjBhUqVLjq2GfPnuWXX37hzTff5H//+x8PP/wwu3fvZufOnUREROQYW373z8q7777LkiVL+O6774iMjGT+/PmsW7eOiIgInJ2d+fzzz7Pdf9u2bbz11lvs2bOHgwcPsm7dOi5evMgdd9zB/Pnz2blzJykpKbz//vvZLr/77rtZsmQJa9eu5d9//83z+yktykyV0Z5jVifpW8LrOjgSpS4XEhJCZGQk8+bNo1evXpetO3/+PCNGjODvv/9GREhOTgbAycmJ2bNnExISwpgxY2jXrl2mx/7f//6HiBAcHEzNmjUJDg4GICgoiMjISEJDQ7ONLb/7Z2bu3Ln4+fnx3XffUa5cOVauXMmWLVto1aoVAAkJCXh7ZzaizX9at26Nn58fAKGhoURGRlKxYkUCAgK49tprARgxYgTvvfcenTt3znR5p06dCAgIoEGDBgAMHz6cmTNn5vr9lCalKiEEenkwf0ybTNcNmbG+iKNRyn59+/blscceY9WqVZw+ffrS8meeeYbOnTvz7bffEhkZSadOnS6t+/vvv/H09CQ6OtNO/ACUL18esBJI+vP01ykpKYB1YzvdxYsXc71/bjVt2pSIiAiioqIICAjAGMOIESOYNGmS3cfIGIuzszMpKSmXqo2ulNVyuPy9q1JWZaRUSTVq1CieffbZS7/A050/f/7STebZs2dftvzBBx9kzZo1nD59moULF+b53DVr1mTv3r2kpaXx7bff5vk49mrevDkzZsygb9++REdHc+ONN7Jw4UJOnLBGsTlz5gz//JPr6YBp1KgRkZGR7N+/H7CuRK6//vpslx86dIgDBw4A1v2bsk4TglLFgJ+fHw8++OBVy5944gmefPJJ2rVrd9mN1ocffpj77ruPa6+9lo8++ogJEyZc+kLNrcmTJ9OnTx9uuOEGfHx88vwe0vXu3Rs/Pz/8/PwYPHhwptu0b9+eqVOn0rt3b7y9vXn55Zfp1q0bISEhdO3alWPHjuX6vG5ubnzyyScMHjyY4OBgnJycuOeee7JdPnPmTHr37k379u2pV69eft96iSfZXU6VNGFhYebKNtDp0quMsqpSUkqp0kREthhjctXeV68QlFJKAZoQlFJK2WhCUEopBWhCUEopZaMJQSmlFKAJQSmllE2h9VQWkY+BPsAJY0xT27L5QEPbJlWAc8aY0Ez2jQRigVQgJbdNp5RSSuVeYQ5dMRt4F5iTvsAYMyT9uYhMA85ns39nY8ypQotOKaXUZQpzCs01IuKf2TqxBhC5GbihsM6vlFIqdxx1D6EDcNwY83cW6w3ws4hsEZHRRRiXUkqVWY4a7XQYkN1IUu2MMdEi4g0sF5E/jTFrMtvQljBGA9Stq0NbK6VUXhX5FYKIuAADgPlZbWOMibY9ngC+BVpns+1MY0yYMSbMy8uroMNVSqkywxFVRl2AP40xUZmtFBEPEamY/hzoBuwqwviUUqpMKrSEICLzgPVAQxGJEpE7bauGckV1kYjUFpGltpc1gd9EZDuwCfjBGLOssOJUSillKbSEYIwZZozxMcaUM8b4GWM+si2/wxjzwRXbRhtjetmeHzTGNLOVIGPMK4UVo1LFQWpqKg8++CBBQUEEBwdz8ODBfB1v2bJlNGzYkPr16zN58uSr1h85coTOnTvTuHFjgoKCmD59ul37jxo1Cm9vb5o2bXrZ9v7+/gQHBxMaGkpYWO66DOUUa07bZXbuixcv0rp1a5o1a0ZQUBDPPfdcns5ZJhljSk1p2bKlycrNH/xubv7g9yzXK+UoL730knnrrbeMMcbMnDnTPProo3k+VkpKigkMDDQHDhwwiYmJJiQkxOzevfuybaKjo82WLVuMMcbExMSYBg0aXNomu/1Xr15ttmzZYoKCgi47Xr169czJkyezjevXX381I0aMyHWsOW2X2bnT0tJMbGysMcaYpKQk07p1a7N+/fpcnbM0ADabXH6H6tAVSjlQXFwc33777aXZ0gICAi5N9ZgXmzZton79+gQGBuLq6srQoUNZtGjRZdv4+PjQokULACpWrEjjxo05evRojvt37NiRatWq5Tm2vMSam+3SiQienp4AJCcnk5ycfGnu5Nweq6xxVLNTpRSwYsUKjhw5QmhoKGDNJ9ylS5ertuvQoQOxsbFXLZ86depl2x89epQ6depceu3n58fGjRuzPH9kZCTbtm0jPDw8T/uD9QXcrVs3RIQxY8YwevR/XYfCw8NJTEzkwoULnDlz5tL7nDJlCrGxsXadK7uYsjp3amoqLVu2ZP/+/YwdOzZf768s0YSglANFRETw4osvcs899wBw1113ERISctV2a9eutet4JpMpcdN/HV/pwoULDBw4kLfeeotKlSrlev9069ato3bt2pw4cYKuXbvSqFEjOnbsCHDpy3bVqlXMnj2b2bNnX9rvq6++sutc2cWU1bmdnZ2JiIjg3Llz9O/fn127dtG0adM8vb+yRKuMlHKgs2fPUqFCBQBSUlL4+eef+d///nfVdh06dCA0NPSqsmLFisu28/Pz48iRI5deR0VFUbt27auOl5yczMCBA7n11lsZMGBArvfPKH29t7c3/fv3Z9OmTXa8c/vPld12OZ27SpUqdOrUiWXLluX5/ZUpub3pUJyL3lRWJc17771n7r33XmOMMa+99pq555578nW85ORkExAQYA4ePHjppumuXbsu2yYtLc3cdttt5sEHH8z1/ocOHbrspvKFCxdMTEzMpedt2rQxP/74Y4HFmt12WZ37xIkT5uzZs8YYY+Lj40379u3NkiVLcnXO0oA83FR2+Jd4QRZNCKqkOXPmjAkPDzfXXHONGT58uImPj8/3MX/44QfToEEDExgYaF5++eXL1vXs2dOsXr3aACY4ONg0a9bMNGvWzPzwww857j906FBTq1Yt4+LiYnx9fc2sWbPMgQMHTEhIiAkJCTFNmjS56nytW7e+dI6MZdmyZXbFevTo0Sy3y+rc27dvN6GhoSY4ONgEBQWZF154we7PpzTJS0IQa7/SISwszGzevDnTdUNmrAdg/pg2RRmSUko5hIhsMbmcS0bvISillAI0ISillLIpVQkhLa30VH8ppVRRK1UJ4cipC44OQSmlSqxSlRDOp7pz/KwmBaWUyotSlRDE2YV7P8y8lZFSSqnslaqE4Jwcx2bTnL+jTjs6FKWUKnFKVULwreICrp7c+8lOR4eilFIlTqlKCJU9yuNz9nf+qtCGP/486uhwlFKqRCnMKTQ/FpETIrIrw7LnReSoiETYSq8s9u0hIvtEZL+ITMjNed+9JRCAsV8cylf8SilV1hTmFcJsoEcmy980xoTaytIrV4qIM/Ae0BNoAgwTkSb2nrRVI18aJazn36pt+Xrtn3kMXSmlyp7CnFN5DXAmD7u2BvYba27lJOBLoF9uDjBrdHO4eI6nf4rXzmpKKWUnR9xDGCciO2xVSlUzWe8LHMnwOsq2LFMiMlpENovI5pMnTwJQ17syN1TYSYJXC6Z+rc1QlVLKHkWdEN4HrgFCgWPAtEy2yWz6oix/5htjZhpjwowxYV5eXpeWvzu6DZyL5P1dVbmYlJK/qJVSqgwo0oRgjDlujEk1xqQBH2JVD10pCqiT4bUfEJ3bc3m6uzKq/r+kVavPQ7N+z1vAShWB1atX4+3tjbOzMwEBAUybltnvJPv5+/sTHBxMaGgoYWGZj348atQovL29adq0aabrU1NTad68OX369Lm07M033yQoKIimTZsybNgwLl68aPf5srJs2TIaNmxI/fr1mTx5cpbbnTt3jkGDBtGoUSMaN27M+vXrsz13Vu/vyJEjdO7cmcaNGxMUFMT06dNzFW+pl9OECcBrQCWgHLASOAUMt2eyBcAf2JXhtU+G5w8DX2ayjwtwEAgAXIHtQJA957tygpzU1DTTYGKEqfPKCXPTO2t1ghxVLM2fP9/cd999l2b5yq969eqZkydPZrvN6tWrzZYtWy6b/SyjadOmmWHDhpnevXsbY4yJiooy/v7+lybwGTx4sPnkk0/sPt+vv/5qRowYcdmylJQUExgYaA4cOHBp9rLdu3dnuv/tt99uPvzwQ2OMMYmJiZc+q6zOndX7i46ONlu2bDHGGBMTE2MaNGiQ5TlLOvIwQY49VwjdjDExQB+sX+/XAo/ntJOIzAPWAw1FJEpE7gReE5GdIrID6GxLCohIbRFZaktQKcA44CdgL7DAGLPbjjiv4uQkvNzFFdyr82d0bF4OoVSh+/TTT+nSpQuVK1cusnN27NiRatWqZbouKiqKH374gbvuuuuy5SkpKSQkJJCSkkJ8fHy+5yLetGkT9evXJzAwEFdXV4YOHcqiRYuu2i4mJoY1a9Zw5513AuDq6kqVKlWyPXZW78/Hx4cWLVoAULFiRRo3bszRo9pnKZ2LHduUsz32AuYZY86IZFbNfzljzLBMFn+UxbbRtuOnv14KXNUkNS9uvr4xb/yylmPOlTl34WJBHFKpAnX//fczaNAgnJ2def/997nlllsuW9+hQwdiY6/+QTN16lS6dOly1XIRoVu3bogIY8aMYfTo0bmK56GHHuK111677Jy+vr489thj1K1bF3d3d7p160a3bt1yPF94eDiJiYlcuHCBM2fOEBoaCsCUKVOIjY2lTp3/aof9/PzYuHHjVfEcPHgQLy8vRo4cyfbt22nZsiXTp0/Hw8MjX+81MjKSbdu2ER4enqvPp1TL6RICmAz8CWzDSg5ewMbcXooURclqTuVdB4+buk+vNYHjl+XyokupwrV3714THBxsfvnlF5OWllYgx0yfh/j48eMmJCTErF69OtPtDh06dFWVypIlS8y9995rjLGqedKrjM6cOWM6d+5sTpw4YZKSkky/fv3M3Llz7T5fZlVGCxYsMHfeeeel13PmzDHjxo27at8//vjDODs7mw0bNhhjjHnggQfMxIkTczx3Zu8vXWxsrGnRooX5+uuvM11fGlAYVUbGmAlAGyDMGJMMxJPLfgGOFhTgTSWJJcW1Em98o81QVfExY8YMHnnkETp37kxWV94dOnQgNDT0qrJixYpMt0+vyvH29qZ///5s2rTJ7njWrVvH4sWL8ff3Z+jQofzyyy8MHz6cFStWEBAQgJeXF+XKlWPAgAH8/vvv+Tqfn58fR47818I8Kioq02ooPz8//Pz8Lv2SHzRoEFu3bs3zuZOTkxk4cCC33norAwYMsCvWsiLHhCAiFYCxWE1GAWoDuWtKUAw0rF0JUi7y9o6qxF9MdnQ4SgFw8eJF/v3332y3Wbt2LREREVeVzKqL4uLiLlX1xMXF8fPPP2fZkigzkyZNIioqisjISL788ktuuOEGPvvsM+rWrcuGDRuIj4/HGMPKlStp3Lix3efr1KkTs2fPvmxZq1at+Pvvvzl06BBJSUl8+eWX9O3b96p9a9WqRZ06ddi3bx8AK1eupEmTJnl6r8YY7rzzTho3bswjjzxi9+dSVthzD+ETYAvQ1vY6CvgK+L6wgioMzk5OlCeBxLQTXPfKXzT2q3LVNv1CfbklvG7RB6fKrMcff5y7776bzz77DB8fH+bMmYOPj0+ej3f8+HH69+8PWDeBb7nlFnr0sEaQ6dWrF7NmzaJ27doMGzaMVatWcerUKfz8/HjhhRcu3bTNTHh4OIMGDaJFixa4uLjQvHlzRo8ezdGjR7M8X/p+iYmJVx1vypQpdO/enXfffZfu3buTmprKqFGjCAoKurRNxnjfeecdbr31VpKSkggMDOSTTz7J9r1m9f7WrVvH3LlzLzVVBXj11Vfp1SvTYdXKHLGqmrLZQGSzMSZMRLYZY5rblm03xjQrkghzISwszGzenHmV0BcbD7Mo4ih/HIol1aUCwV6peLi5Xlq/51gMTXwqMX9Mm6IKV6nLjBgxgptvvpnevXs7OhRVCojIFmNMrmpz7LlCSBIRd2y9hUXkGuDqlF/M3RJel1vC67Jmxz8M/7km/57awuaX/vvyHzJjvQOjU2Xd999/T1xcXKbVQEoVFXv6ITwHLAPqiMjnWJ3TnijUqApRx5B6hKZs4GSNdrz/fYSjw1EKgD59+rBw4ULKly/v6FBUGWZPK6PlwADgDmAeVmujVYUbVuH6dGw4nD/Ca3946A1mpZSyyTIhiEgj22MLoB7WYHTRQF3bshKrakV3RjeIJq1aA+58b52jw1FKqWIhu3sIjwCjyXxEUgPcUCgRFZEnh7Tmy2c3sq5yKzbsOZLzDkopVcplmRCMMaNFxAmYaIwpdT+jnZyEj4f7MXCx4a4vjhFUJ+d9lFKqNMv2HoKxhqmeWkSxFLlWjXzpWG4zsV6tiTwe4+hwlFLKoexpZfSziAwUe0a0K4FmjW2P8+k/+feiG4nJOpGOUqrssqcfwiOAB5AiIhexZjQzxphKhRpZEXFzdeGVTmlM+NWVvdHnHB1OsZbeuS872ttbqZLLnmanFY0xTsYYV2NMJdvrUpEM0g27oQnuqee46FyFWT9ud3Q4xdaiiKPsOZZ11dqeYzE5JgylVPGV4xWCiKw0xtyY07KSrolvRbYcTeTljTUZ3OEilT3dHB1SsZTd8B7a21upki27fghuIlINqCEiVUWkmq34Y414mi0R+VhETojIrgzLXheRP0Vkh4h8KyJVstg30jazWoSIFMl41eVcnPHzSMJUDWTY2xuK4pRKKVWsZFdlNAZrlNNGtsf0sgh4z45jzwZ6XLFsOdDUGBMC/AU8mc3+nY0xobkdnCk//GpUxOfMWna5teerNXuL6rRKKVUsZJkQjDHTjTEBwGPGmEBjTICtNDPGvJvTgY0xa4AzVyz72VhzJgNsAPzyE3xhmH9vMCScZvwKJx3WQilVpthzU/kdEWkrIreIyO3ppQDOPQr4MavTYjV33SIi2U6SKiKjRWSziGw+efJkvoOqV6sKdwdGklq9IcPeKnX98ZRSKkv2zJg2F6tzWnugla3kqxpHRJ4GUoDPs9iknTGmBdATGCsiHbM6ljFmpjEmzBgT5uXllZ+wLnlmWDjep39jm2s7FqzeUyDHVEqp4s6efghhQBOT00w6dhKREUAf4MasjmmMibY9nhCRb4HWwJqCOL+9Ft7XlA6zTjH+l3L0bJVIxQo6LLFSqnSzp6fyLqBWQZxMRHoA44G+xpj4LLbxEJGK6c+BbrYYilS9WlW4v1EUadUaMPjN34v69EopVeTsuUKoAewRkU1kmCnNGHP1bNgZiMg8oBNWs9UorIl2ngTKA8ttI2FsMMbcIyK1gVnGmF5ATeBb23oX4AtjzLLcvrGC8PigVnz73Fr2VOnIuE/WczIp++21l65SqiSzJyE8n5cDG2OGZbL4oyy2jQZ62Z4fBIrNfM0L729Gm/+LZsnJi1Qsl0hQ7cqZbpfeg1cTglKqpMoxIRhjVotIPaCBMWaFiFQAnAs/tOKhdo1KPNs6khfXlScxKVF76SqlSi17WhndDSwEZtgW+QLfFWJMxc6dPULwSDtPUrnKPPeZfvErpUone24qjwXaATEAxpi/Ae/CDKo4CvKrjCTH88mRa9lx4F9Hh6OUUgXOnoSQaIy5dDtVRFywOo6VKU5OwrU1nMDFnZs/OUpKapqjQ1JKqQJlT0JYLSJPAe4i0hX4ClhSuGEVT1U93fif5xYSvFpyyxurHR2OUkoVKHsSwgTgJLATa8C7pcDEwgyqOHtnTHuqn1rPBucOzFm+09HhKKVUgbFnLKM0Y8yHxpjBxphBtudlrsoonZOTsPiBJsiFYzyzvgqHT5x3dEhKKVUg7LlCUFfw86rMy21jMJ4+9H17N2lpZTY/KqVKEU0IeXRblyDamXWc9WrLne8W6TBLSilVKLKbMW2u7fHBogunZPns4Y5UOrmJlcltOHEuztHhKKVUvmR3hdDS1kN51BVTaFazTa1Z5jk7CUvGXYvEn+DgeRcuJqXkvJNSShVT2SWED4BlXD2F5hagSOY5Lgn8a1VhcvtYcHZl59F47Z+glCqxsptC821jTGPg4yum0AwwxgQWYYzF3rDOjakqMaS6VmLw66scHY5SSuWJPc1O7xWRZiIyzlZCiiKwkqahb2Vck2PYUr4TU+ZvdHQ4SimVa/YMbvcA1lSX3rbyuYjcX9iBlUQhdTxxOfc37x1qyOqISEeHo5RSuWJPs9O7gHBjzLPGmGeB64C7CzesksnF2YnPh3gChju+SeLE2QuODkkppexmT0IQIDXD61Tbsux3EvlYRE6IyK4My6qJyHIR+dv2WDWLfXuIyD4R2S8iE+yIsdho08SXRxodIrVqfbpN3U6q3mRWSpUQ9iSET4CNIvK8iDwPbCCLmc+uMBvoccWyCcBKY0wDYKXt9WVExBl4D+gJNAGGiUgTO85XbDx0Uws6yjrOerdjwGu/OjocpZSyiz03ld8ARgJngLPASGPMW3bst8a2T0b9gE9tzz8Fbspk19bAfmPMQduw21/a9itR5jzQnlpnfyeiwo08M2edo8NRSqkc2TV0hTFmq60Z6nRjzLZ8nK+mMeaY7ZjHyHyiHV/gSIbXUbZlmRKR0SKyWUQ2nzx5Mh+hFSwnJ2H54y1xPb2bT/9txsLVexwdklJKZas4jmWU2f2JLEePM8bMNMaEGWPCvLy8CjGs3KvkUZ5Fd9ZEkmJ5ZE1ltv111NEhKaVUloo6IRwXER8A2+OJTLaJAupkeO0HRBdBbIWiSb0avH1jHLhVYeCcsxw/E+vokJRSKlPZJgQRcRaRFQV4vsXACNvzEcCiTLb5A2ggIgEi4goMte1XYvVrU5+xAftIqd6ELlN3kpSsYx4ppYqfbBOCMSYViBeRyrk9sIjMA9YDDUUkSkTuBCYDXUXkb6Cr7TUiUltEltrOmQKMA34C9gILjDG7c3v+4mb8oBZ0K7+BmJpt6f7yKsrwHENKqWLKxY5tLgI7RWQ5cGmMZ2PMA9ntZIwZlsWqGzPZNhroleH1UqypOkuVD+9tS5dX17G/WheGT13O5493dXRISil1iT33EH4AngHWcPmIpyqXRGDZ+DZUO7uZtS43Mn7WKkeHpJRSl+R4hWCM+VRE3IG6xph9RRBTqVbOxYlVj4fQ+rW9zEu9Dr9vNnL/gHBHh6WUUnYNbvc/IAJrbgREJFRESvRNXker7OHK8rF1cY7/l9f2NWDhqp2ODkkppey6h/A8Vu/hVQDGmAgRCSjEmMqEut4VWTgsgQELknn4Ny+2H9/EX+dSc94xG/1CfbklvG4BRVj6fLHxMIsisu8Lop+hKsvsuYeQYow5f8UybSJTAFo28Obj3kmIS3nmRMSxM+pcno+151hMjl92Zd2iiKPsORaT5Xr9DFVZZ88Vwi4RuQVwFpEGwAPA74UbVtlxY2gdpsUd4LEfyxGXnMaL3f1p6O+T6+MMmbG+EKIrfZr4VGL+mDaZrtPPUJV19lwh3A8EAYnAPCAGeKgQYypzBre7hnqeiVDOnR4z/iXyaPEZk0kpVXbYM9ppvDHmaaz+A52NMU8bYy4Wfmhli0/VCvi6x5PqFcwNbx8i6vhpR4eklCpj7Gll1EpEdgI7sDqobReRloUfWtlTp7oHo+vuJaVWGJ2m7ePE6XOODkkpVYbYU2X0EXCfMcbfGOMPjMWaNEcVgomDgrndZzfJvm1pN3kHJ06ddXRISqkywp6EEGuMWZv+whjzG6BDdhail4cGM6TGTpL8OtJu8g7+PaHVR0qpwpdlQhCRFiLSAtgkIjNEpJOIXC8i/4etT4IqPK/dGsxQr10k1bme9lN2cuz4KUeHpJQq5bJrdjrtitfPZXiu/RCKwJRbmuL05V6+oBPtX1vFyocb4u+X+yapSilljywTgjGmc1EGUhrsORaTbVv2vPSCnTS0Mc5f7WMunej0xlqW3ZdEo/r1CuX8OfXk3XMshiY+lewPvoiV9PiVcrQcO6aJSBXgdsA/4/Y5DX9d1vQLzXLaZ4BLPWTzMizCy4Mb4rF4P++bdnT7YBPf3ZFIi6bXFvj503vyZvWl2cSnUo7ncaSSHr9SjmZPT+WlwAZgJ5BWuOGUXLeE1832yza/vWCf7Fsfz+WHed2E0ffTCOYO2ErnNi0K/PzZ9eQtCUp6/Eo5kj0Jwc0Y80hBnVBEGgLzMywKBJ41xryVYZtOWNNrHrIt+sYY82JBxVBS3d+1Lp5u0Ty/KZjhi/fzxulfGdJHa/aUUgXDnoQwV0TuBr7HGr4CAGPMmbyc0DanQihYczYDR4FvM9l0rTGmT17OUZqN7FCbahVP88CKAB7Z5M6xk9/w0MgBjg5LKVUK2NMPIQl4HWt+5PTZ0jYX0PlvBA4YY/4poOOVCf1Cq/NF/zRcK9Zg6pE2PDLpQ52jWSmVb/YkhEeA+raeygG2ElhA5x+KNWBeZtrYhsn4UUSCsjqAiIwWkc0isvnkybIzKFy7azz5/jY33N3L85XcTL9xU0hMTMx5R6WUyoI9CWE3EF/QJxYRV6Av8FUmq7cC9YwxzYB3gO+yOo4xZqYxJswYE+bl5VXQYRZrQbVcWTWmKtVdE4kIfIS2Iydz+rT2alZK5Y09CSEViLD1Vn47vRTAuXsCW40xx69cYYyJMcZcsD1fCpQTkRoFcM5Sx7eSsPo+bwLdz3Eq/DnC7pnFjh06JadSKvfsSQjfAa9gTYqzJUPJr2FkUV0kIrVERGzPW9vi1J++WajsBstGe9Om2mnS2o+n65Q/WLDwG0eHpZQqYXJsZWSM+bSgTyoiFYCuwJgMy+6xne8DYBBwr4ikAAnAUFMEd02z6+lbEL1cC/P4bi7wxfDqPPdzLHMYxbhflvNHxEtMfuEpnJ2d83zcolQW5jwuC+9RlVz29FQ+RCZjF+XnxrIxJh6ofsWyDzI8fxd4N6/Hz4ucerDmt5drYR8fwEngpe4VaVwrmadMZz4/7svGm0aw4MOp+TpuUcmpp3F+ensXF2XhPaqSy55+CGEZnrsBg4FqhROO4+TU07e4H/+yczUrh381GPnNNRyq9DYtbhpNcP87qVKlSpGcPz/KwpzHZeE9qpLJnik0T2coR209im8o/NBUfrStAz+PKI+/VwXKDfmSXZEniYyMJDk52dGhKaWKKXum0GyRoYTZ6vorFkFsKp/qVYGlt7lxfT1w8WrA0QQ32l9/I/v373d0aEqpYsieVkbTMpRJQEvg5sIMShWciuXhk5tc8K0IzpVqcbT9W7S4vg8zZ87U3s1KqcvY08pIR08r4ZydoE4l8HSFIy6hlHtwAw9MHciCBQv48MMPCQgIcHSIxdb58+eJjo4mJiaG2NhYEhIS8PDwwNPTk0qVKlGvXj3c3d0dHaZSBcKeVkblgYFcPR9CmR99tKSp6gbv3+zE6CWVMWOXs2PlqzQNDuaVl19m3Lhxjg7PodLS0oiLiyMmJoa77rqLPXv2sH//fnIaDkVE8Pf3p3HjxrRq1Ypu3brRunVrXFzsaa+hVPFiz1/tIuA8Vmc0HSynhLumKiweKjz1i/AdE/Fq1oNHJ/bkk08+ofbwKVSqVNnRIRaZqKgoli1bxrJly1ixYgVuvSYAsPunbwkJCeGmm26iQYMG+Pn5UblyZSpWrIibmxvx8fFcuHCB8+fPs3//fvbu3cuePXtYtmwZL7zwApUqVaJ79+6MHDmSbt26lZh+IErZkxD8jDE9Cj0SVWQ8XOGt7hDuC8+tbkmTSYc5O+92jm+LoGbNmkRF1cHPz8/RYRaKv/76i6+//pqFCxeydetWAHx9fRk0aBD76zaiUqVKLPr8FLaO8rly5swZfvnlF3766Se+++47vvrqK/z8/Bg5ciTjxo3D29u7oN+OUgXKnoTwu4gEG2N0gJxSRARuCYaQmsK4H925MHQB1Y7/yIl/dtGgQQMeeughxo8fXyL6LuTk8OHDfPnll3zxxRds374dgPDwcKZMmUKvXr0ICgpCRBgyYz17jsUwdOaGLI+VXS/iatWqMWjQIAYNGsR7773HkiVL+Oijj3j55ZeZNm0aY8eOJbluH8qVK1co71Op/LKnlVF7YIuI7BORHSKyU0R2FHZgqmg09YYfhsGARkKcc1WqXtuO3sNGM3nyZAIDA3n11Vc5f/68o8PMtbNnz/Lhhx9y/fXXU69ePcaPH4+7uztvvfUWhw8fZsOGDTzxxBM0bdr00tVAv1DfbIcP2XMsJsdhJ9K5uroycOBAli5dyt69exkwYADTpk1jw4YNREZGkpCQUCDvU6mCZM8VQs9Cj0I5lIcrTOsGEbvh0Dkn/mw1nTdufpif376Pp59+milTpjBu3DgefPDBYl3tkZaWxpkzZxg4cCDff/89SUlJNGrUiJdeeolbbrmFwMDsR1sprHmxGzZsyNy5c5k4cSKD/u83/vnnH5o2bcrbb79N796983RMpQqDPT2V/8msFEVwqmjVqAAhNaFRDXhrnz+B45ayakMEXbt2ZdKkSdSpU4fbbruN9evXF5s+DElJSSxbtoyRI0fy+++/s3v3btatW8d9993H5s2b2bNnDxMnTswxGRSFhg0b0qRJE5o1a4arqyt9+vRh0KBBObZkUqqo2FNlpMqQ8s4wfyA80RaWHYDHdjbj3tcXsmfPHkaPHs2iRYto27YtoaGhTJkyhcjIyCKPMTY2lu+++44RI0ZQs2ZNevbsyTfffIOXlxchISFERUXx5ptv0rJlyzzdHC5sVapUYfv27UyaNIklS5YQHBzM0qVLHR2WUpoQ1NWcnWBsK1g0BKq5wajFMOtII16Z+g7R0dF88MEHVKhQgQkTJhAQEMB1113HSy+9RGxsDJkMjJtvycnJbNy4kSNHjrB9+3aqV69O//79Wbx4Mf369WPJkiUcP36chg0bUrVq1RLRB8DV1ZUJEybwxx9/4O3tTe/evbnvvvtIS0tzdGiqDCv+/3KUwzT1hsVDYfomeH8zrD0Mr97gyZgxYxgzZgyHDh1iwYIFLFy4kOeeew7voa/i4uJCz+9eokWLFrRs2ZKgoCDq1q1rd2/elJQUNm/ezM6dO9m5cyfbtm1j48aNJCQkUHPYJDw8PHj44Yfp2bMn7dq1K/EtdkJCQti0aRMTJ05k2rRpXDO6CUFNspxCXKlCpQlBZau8i1V91D0QHl9hXS30awjPdYSAgADGjx/P+PHjOXXqFDd/8Dtnzp4hens0y5cvJzU19dJxvL29qV27Np6ennh4eFChQgWSk5NJSEggISGBkydPcqH1XaSlpdFq3JMAuLm50bRpU0aPHk379u35/JgXrq6uTBlzv6M+jkLh5ubG1KlT6dixI/ct3MeWLVv46ackunfv7ujQVBnjkIQgIpFALNZ8zSnGmLAr1gswHegFxAN3GGO2FnWc6j/NasH3w6wrhXc2wZp/YGIHGNjY6tNQo0YNvLy88PLyYv607SQkJLBjxw7++usvDh8+zD///EN0dDRxcXGcOXOGqKgoXF1dcXd3x93dnebNm3Ogdm3Kly/PewsXEhwczDXXXHNZL9+vSvlcAX379mXWwV/YvXs3PXv2ZPLkyTz++OPF8j6IKp0ceYXQ2RhzKot1PYEGthIOvG97VA7k6gwPhkOPa+DJX+DR5fDVXnj1BmtIjIzc3d0JDw8nPNz+/23pzToHDsx88piywN3dnRYtWhAWM5jx48ezb98+3n//fVxdXR0dmioDiutN5X7AHGPZAFQRER9HB6UsDWvAwsEw6QbYcxJ6fA6v/w5pxaMlaonn5OTEvHnzeOaZZ/j444/p1q0bp0+fdnRYqgxw1BWCAX4WEQPMMMbMvGK9L3Akw+so27JjRRSfyoGTbeiLroHwym/w7h9Q/hz4VwFjrGqk0mjPsZh8TXOZ3XzKGTk5OfHiiy/SsGFDRo0aRbt27Vi2bBn+/v55PrdSOXHUFUI7Y0wLrKqhsSLS8Yr1mX2dZPr7U0RGi8hmEdmsHXyKnpeHNVDegkHgLLDvNNz+HfxVCn/Q5jS0hT2a+FSiX6iv3dvfeuutrFixguPHj9OmTRsiIiLydX6lsuOQKwRjTLTt8YSIfAu0BtZk2CQKqJPhtR8QncWxZgIzAcLCwrTSwkHCfSG4Jhy/ABH/WtVItzeDh8KhipujoysYOQ1tUVg6dOjAb7/9Ro8ePejYsSPffvstN954Y5HHoUq/Ir9CEBEPEamY/hzoBuy6YrPFwO1iuQ44b4zR6qJiToBanrBqBAxrCp9uh46z4aNtkJSa094qO0FBQaxfv5569erRq1cvvv76a0eHpEohR1QZ1QR+E5HtwCbgB2PMMhG5R0TusW2zFDgI7Ac+BO5zQJwqj6pXgFdugKXDrLGRXlwDXebCD39b9xdU3vj5+bFmzRrCwsK4+eab+eijjxwdkiplirzKyBhzEGiWyfIPMjw3wNiijEsVvMZe8Fl/WP0PvLIW7lsKzWrC+HbQrk7O+6urVa1alZ9//plBgwZx1113cebMGR5//HFHh6VKieLa7FSVItfXgx9vgde7wMk4uOUbuO1b2HHc0ZGVTB4eHixatIghQ4bwxBNP8PTTTxeb0WdVyaZDV6gi4ewENwdB34Ywdwe89wf870ur2eoj10ETL0dHWLK4urry+eefU7lyZV599VViYmKYPn06Tk76G0/lnSYEVaTcXODuFjA0CD7ZDh9ugZ5fQM/6EJ8MFUr2WHVFytnZmQ8++IBKlSoxdepUYmJi+Oijj0rEaK+qeNK/HOUQFcvDA61hRDOYtRU+iYDYE9Zw2ztPQHDxnZjtcnFx8M8/EBkJFy78t9zNDerVA39/qFy50E4vIrz22mtUrlyZZ555htjYWObNm0f58uUL7Zyq9NKEoApUXnvyNgJ2pMUQk1iJPvOs+w5jW0Hr2pf3es7u+Pb2As4zY2DPHli1yirr1sExO1pDV6sGbdrA9ddDp07QogVkGLTvSnn6DL1upN0zjfh9/346jnySXz58CQ8Pj9wdQ5V5mhBUgclND9wrOTtB8zqV6N7UlwRXq+/CzQuhpQ+MaWnda8jp+LntBWy3Q4fgs8+s8tdf1rI6deCGGyAoyLoKqFcPqlT5b5+MVw5//glr18IPP1jrfH3hlltg+HAICbnsVPmJ39fXl5Mp5fknErp37873339PlYwxKZUDKU2tE8LCwszmzZsdHUaxlP6Lc/6YrEcStWebonIxBebvhplbISoGAqtY9x4GNLbuQxQ6Y2DlSnjtNVi+3FrWqRMMGwZdu1pJILcDNh07Br/+Cl9+CT/+CCkpEBYGjz8OAwdme9VgryEz1nPq1CnWvjCQJk2a8NNPP1GzZs18H1eVPCKy5cqpBXKiTRJUseTmYt1fWD0C3ukBFVytIbfbfAxT18PxuEI6sTGwcKH1Rd21K+zaBS+/bP3a//VXGD0aAgLyNnqfj491ZbB4sZUc3nkHYmJgyBBo2BBmzIDk5Hy/hRo1avD999/z999/0759e4fMe61KJk0IqlhzcbKaqn4/FOYNsKqQ3t0E7T6GB5fB1mMF2Pv5t9+suv7Bg60bxB9+aFUXPf001C3gMYxq1IBx46x7El9/DdWrwz33QNOmsGhRvt9Ut27dWLFiBadPn6Zdu3bs3r27gAJXpZkmBFUiiEDbOjDrf9ZYScNDYMUh6L8A+n4JX+2xqpny5PBhq8qmQwc4cgQ+/tj6or7rLijs1jrOzjBgAGzYAEuWgJMT3HSTVT21Y0e+Dt2mTRvWrFmDMYYOHTqwbt26AglZlV6aEFSJ418Fnr8eNt4JL3eGhBR4bDm0mgUvrIa/z9h5oJQUePNNaNIEli2Dl16ybhqPHFkg9fm5IgJ9+sDOnfD++1ZCatkSJkyA+Pg8H7Zp06asW7eOGjVq0KVLFxYtWlSAQavSRhOCKrE8XeG2EFg+3JqPoVM9qxd0l7kwcIF11RCfVZX89u0QHg6PPGI1B929GyZOBEc31XRxsaqO/vwTbr8dpkyxqpFWrMjzIQMCAli3bh0hISEMGDCAmTOvnI9KKYsmBFXiiVjzMbzTEzbcCU+2gzMJ/101TFgJm6Nt1fIpKTBpErRqBUePwvz58P33Vquh4qR6dfjoI6u/g6urdYN73DirOWseeHl58csvv9CjRw/GjBnDU089RVpaWsHGrEo8TQiqVKlRAe4Jg19uh68GQfdr4Ls/YeBX0PmjJN4Z8xFHps6w6ul374abby7e831efz1s2wYPPQTvvQfNm1v3G/IgfVC8MWPGMGnSJIYNG0ZCQkLBxqtKNO2YpkolEWjta5UXrzf8OOd3Fu5MY2rIGKaGjKF1bcNN0UKvClDV3dHR5sDd3brX0bcv3HEHtG8PL74I48dneq8jp57OfUc+SYMGDXj88cc5fPgwixYtwtvbGivki42HWRRxNNtw+oX6OmTmOFX49ApBlW4xMXjeeRuDR7dn/vZnWNfjGI+3gdMJwlO/QNgsuGMRfPMnxCY6OtgcdO5stTwaNMhqCtu1q1XtlUFO8z7vORbD4ohoHn30Ub7++mu2b99OWFgYW7ZsAWBRxFH2HIvJdv+cEoYquYr8CkFE6gBzgFpAGjDTGDP9im06AYuAQ7ZF3xhjXizCMFVpsGWL1ekrMtJqQfTkk/g5OzMOa5yk3Sdh8V+wZB88HAnlnaFjPejdALoEWAPwFTuVK8O8edC9u3VPoVkzmDMHevUCcp73OeOVQ//+/fntt9/o378/7du3Z9asWUAgTXwqZdlbPS/jVKmSwxFXCCnAo8aYxsB1wFgRaZLJdmuNMaG2oslA2c8YePddaNsWkpJg9WqrBVGG6hURaOoNT7WHdaNg4WC4NdgaafWhn6DFh9aVw5e74HTeW30WDhGraezWreDnB717W81T89DLuUWLFmzevJnw8HCGDx/OgQMHdLKdMqzIE4Ix5pgxZqvteSywFyiEEclUmXT+vHWj+P77rSqVbdugXbtsd3ESaFUbnrse1o+Crwdbw2b8fQbGr7SqlQZ9BTO2wKGzRfQ+7NGwIaxfD2PGWM1TO3eGqKhcH8bLy4vly5fzwAMPEBUVRUREBIcPHy6EgFVx59B7CCLiDzQHNmayuo2IbBeRH0UkqGgjUyXS1q3W0NLffmsNSrd4sdV8MxecBMJqw8QO8Nsd8MMwuL8VxCXBq79BpznQ+VN4aQ38fgSSUwvnrdjN3R0++AC++MLqWxEaanWyy6Vy5coxffp0mjRpQlxcHKGhoSxZsqTg41XFmsMSgoh4Al8DDxljrryLtRWoZ4xpBrwDfJfNcUaLyGYR2Xzy5MlCi1cVY8ZYX4pt2kBiolVF9Pjj1jAQ+ZBerfRIG/jxVvhtpNVDuk5lmLMDhn0DoTNh9PfwxU44mvW92MI3bBhs3gy1a0PPntZN55Tcj+Xh5eVFy5Yt8ff3p2/fvtx3333E5bHvgyp5HJIQRKQcVjL43BjzzZXrjTExxpgLtudLgXIiUiOzYxljZhpjwowxYV5eOjFvmRMTY40geu+91vwEERE5VhHlVZ1KMDIU5twEEaNhZh/oey3sPG6NxNr2E7hhDjy/GlYehAtJhRJG1ho2tPoo3HknvPoqdOkC0dG5Poy7uzu///47jz76KB988AHNmzdnQx77PqiSxRGtjAT4CNhrjHkji21qAceNMUZEWmMlrtNFGKYqCbZts+4XHDwIr7xi3VgtoknmPVytTm/dr7EuUP4+A2v+gTWHYd4ua0pQFydoVhPa1bEG5mteqwjmcqhQAWbNgo4drSQZGmpN7NOtW64O4+bmxtSpU+nTpw8jRoygXbt2PPLII6QF9sPJqYjHeVJFxhEd09oBtwE7RSTCtuwpoC6AMeYDYBBwr4ikAAnAUKNNH1S69Cqihx+2hpFetcoaqdRBRODa6la5q4U16urmaOsew7oj8O4f8PYmq1lrCx+4zhfC/Qo5Qdx+uzU8x+DB0KMHPPkkvPBCrg/TqVMnduzYwWOPPcbUqVOpN9KXBg0aFELAqjjQGdPKiCEz1uc453D6+uIwY1qWTp+2hqX+7jvri27OHCjmVYXnE+GPo7A+yip7ToIByjlBSE2rhVOr2tZcDwXeazo+Hh54wBoXqU0bhgx5hT1nk7L8O8jub2DNmjWMmLONVM9aVDYXCAwMxM3NLdchaU/n7BVUb/G8zJimQ1eUEfbM1VtocxIXlNWr4dZb4cQJmDbNGt+niKqI8qNyeegSaBWwEsTmaNh4FP6ItuaP/sDqKMw1Va2riJY+1hVEg2rWfNN5ll6F1KULjBlDv3lvQ5+RWW6e3d9Ax44decHZl//74Q8OHz7N6dOn8POrQ926dXG2c7jw9F7QmhCylt5bPLukDYXzGeoVgir+EhPh2Wfh9dehfn1rTuIWLRwdVYG5mALbj1tJYusx2HIMzl601nmUs64imteyHpvVBB/PPI7Hd+iQdQN+wwYYPtyawrNKlTzFHBUVxVNPPcXcuXPx8vJiwoQJ3Hvvvbi7Z3+JU5zm7S6ucvqM7P0M9QpBlT47d1pfXjt2wN13wxtvgKeno6MqUG4u1vDd4bYf5sbAoXMQ8S9s/dd6nLkVUmyjVXtVsJrDBttKU287k0RAAKxZY92Af/ll64pr9myrdVYu+fn5MWfOHMaNG8fTTz/No48+yuuvv86TTz7JXXfdRYUKFXJ9TOV4mhBU8ZScbF0RvPACVK1qTS/Zp4+joyoSIhBY1SoDGlvLLqbA3lPWlcTO49YQG6v/gTTbBX5VNwjygibppYa1f7kra3LKlYPnn7fGPrrtNrjxRmtMpEmT8pRoW7duzfLly1mzZg3PPvssDz74IC+++CL33XcfY8eOpWbNmvn6LFTR0oSgip+tW2HUKKvn7aBB8H//V+xvHBc2Nxer2qh5rf+WxSdbSWL3CWugvl0n4dPtkGjrPe3qDPWrQqMa0LAGNKxuFR9PkNatrWa7Tz5pVR0tXgwzZ1qD5uVBx44d+fXXX/ntt9+YNm0aL7/8Mq+99ho333wzY8aMoW3btkhxnndCAZoQVHESG2uN8//mm1YC+OYb6N/f0VEVWxXKWTefW/r8tyw5FQ6ehd2nYN8pK2H8HmUN752uoivUrwYNqlWgwYjp1O82hvovPIhvz144D7/VujLLwy97EaFDhw506NCBffv2MX36dD777DPmzp1LkyZNGDVqFInlWlK+fHEcRlaB3lRWxYExsGCBNb9xdLTV0/b1162qIlUgzl2Ev05bZd9pqyPd/jNwMsNIruVNMgH//kng6QMENq3NNV1aElDdmYAqUCX3rUsBuHDhAgsWLGDGjBls2rSJmsMmUaVKFR5uJvTt2xcfH5+cD1LG6E1lVXZt2QKPPWZ1LmvRAr7+Gq67ztFRlTpV3P6bQS6jswmw/ywcOAP7z5bjwJEA9kZW5afUWqSucL5s/4AqUK8y+Nse69qKV4Wsb2h7enoyatQoRo0axV9//cXts7dy/MRx7rnnIe655x5at25N37596d69O82bN7e7+aoqHJoQlGMcOmTNUfDFF1Zv4//7Pxg9OtMpIVXhqeoOrdytjnEWTzAeJH3zHYdffY/IRHcOtevDoRsHEulSg83RsGif1bEunbsL+FWyxnqqU9l69EsvFa1kIgLXXnst9eqdpl69eix5aCeLFy9m0aJFTJw4kYkTJ1K1alVuuOEGOnXqRLt27QgJCdEEUcQ0IaiiFRkJkyfDxx+Di4s1KucTT0ClrHtQqyImguvA/tTv04v6778PLz0FM++xbvBPnEhiUDOiYuHwefjnvPUYFQNHzsOm6KsH9atQDmpXBN+KcPAcxFyI4YHlUN6zM7Vu70y/4UmcP3eOs2fP8sfZs/z2+wXiZj0JB9bRunVrWrVqRVhY2KVRWHO6Oe3oeaEdff780ISgisaff1r3BebMsXoX33mnlQz8/BwdmcpK+fJWb/A77oCpU63WSAsXUr5fP6558kmuCQ+/ahdjICYRjsRYJTrWGhY8KtZ6Hufsy0WBA5dNNOSKi3jjWsWbKtUhLv48Na8Nwy9pBIf3bOadBatJ/OAzUmOPU9HDneDgYIKDg2ncuDGNGjWiUaNG1KlTBydbr3VH9vQtDufPD00IqvCkpcHPP8P06dakLW5u1gicTzyhiaAkqVLF6sj26KPw9tvw1luwaBGEh1sJY+BAq38DVtVQZTerNPXO7GB1uZhSl2OxcOwC/HvBShT/xsFx2+s/963nfBqcq9kFag6jRuf0fQ2uKRc4HneCyFOHWfhHFKm/7CA1djlOCWfwruhE3eoVOF2nDRXLl2OIVxz16tWjTp061KhR41LCKIp5oUvqvNSaEFTBi4qCTz+1esHu3w+1alnNSceMAe9MvyVUSVC1Kjz3nNUabPZsKzkMG2b9/739dmue50aNcjyMmwsEVLVKZobMsO5RvH0rnIj7rxyPE07GVeRkfEWOxwVy4kIap+Ih2Vj3GRKBvwFOW1+4D/zZkLQtJ0m98BfEr8edBDydk0iuVIlyTnDva1/hU9Wdul6e+NeqQkDt6tT09srTgH2lhSYEVTBOnbJGIF2wAFassOoOrr/e+gK5+WZwdXV0hKqgVKxozVk9diz8+KPVoW3aNGva0vBwGDLEumqom/cqEQFqeVol6y2cMQZik+B0PJyKh1MJMGUxJCan0SwghWPnPTgdX4GY5GuIN+7EunjC2T+4CCwt3wbigX+sYtLSSEs4CwlROKdcwDUtHncS8XBJoWK5NCq7CdU9nKnuWQ6vSuXxqVoBn2oV8KtRidpelXErX/L/xjUhqLwxBvbsgZ9+gh9+sMbFSU2Fa66xWg+NGGE9V6WXkxP07m2Vf/+Fzz+HuXOtK4hHHoHWra3hRnr0gJYtC2VkWhGoVN4q6Vccsz0AnPjg5quvRo2BwTOsDnwTel8k8t9zHD4ZS/TZi5yKTeZMquG8EeJTXUigBrFOFThfzhPjWglxcoI0IMZWojIcNy0Nk3QeSbqAU3wCYlJp8vQm3CQZd+cUPFzS8HCFSuWF/ecu4uIkvPHNZqp7lqd6JTe8KrvjXcUDryoeBf4Z5YYmBGUfY6wbw7/9ZpVff4UjR6x1TZpYs5UNGgTNmuVxKE5VotWqZd1jePRRq5rw66+t8txz1ki1NWpA587Qvr1VQkKsVmZFTAScBZxdoI2/G238awG1ctwvzcDZuGSOnDjP0dMXOHYmjuPnL3L6QjJn4lI4l2i4kAxxqc6cIAUjQqJzReKd3Ukr5wGunkg520iwaeshDaYfyaLP2Kk0SEuj3ktROKVcxDntIi4miXIk4yopxKbE4iSG3q+uwr0ceLgKFVyFiuWd8XRzoZK7C5Xcy+Xp83FIQhCRHsB0wBmYZYyZfMV6sa3vhXVRd4cxZmuRB1pWxcXBvn2wd6813s3WrdbjuXPWem9va4ayZ56xxr7JR9WAKoXq14fx461y8iQsX241Kli9Gr76ytqmQgVres8WLaB5c+tHRePGULmyQ0PPipNAdc9yVPesQWhgptO7X/JfT+LGly2PT0zh3zOxjP0ykeSUNG5ttJMzF5I4E5fE+fgULiSmEZto2JuWQKoRqiRHkmicSTLlSBZXksSDGGc3SIsEcWanexskfTrTJFuJzd/7dMScys7Ae0BXrIuuP0RksTFmT4bNegINbCUceN/2qPLLGLhwAY4dg6NHrXLkiDUv8aFD1q+7f/75b/vy5a1f/UOHWlMyduhg/YPXqwBlDy8vaw6GW26xXh85AuvWwfr11o+M2bPh3Xf/297HB/o9A+5u8OJy68eGr69Vate2EkYJ/durUN6FQJ+qVKpgjeU0sntwptvlPHRFNADz7nYi7mIip87HczomgTOxF4lNSOJ8XBLn45N5JA8xOuIKoTWw3xhzEEBEvgT6ARkTQj9gjm0e5Q0iUkVEfIwxx4o+3EJgzH8lLe2/kpp6dUlJsUpyslWSkqySmAgXL1qPCQnWVInx8dav+9hYq8TEWL/qz561ysmTVklMvDomb28IDIR27aw+Ao0bW+Xaay81KVQq3+rUsX5cDB1qvU5LgwMHrKvR9ALW3+tbz1v/RjIqV86qfvLygmrVrCaxVataHRsrVrRK3DXW/Yr5860rkfLlreLmZjVuiI8DcbI6Sbq4WMXZ+eri5GQVkf8ei1EycnISKlYoT8UK5QnwubrJVklJCL7AkQyvo7j6139m2/gC2SeEbdusPwh72TuwX8btrtznynXprzM+XlmKgqen9VlUrWqV2rWtX/peXlbx8bGW+fpafQI8HHszS5VRTk7QoIFV+va1ltnm/x7y3lpISoRE2w+gpCTbD6MkSEr+78dSSgqcToETqUAae7xTaHLiEDz/ZObnHDaJPd4BDHlqU/5il0v/ucwe7wCanIyEhzJMPJQxkQx4nj1e/gwZvjaz3dlTw58mpyLh0a6Z73/Ts+ypUY8ht/2Wv/gzUeSjnYrIYKC7MeYu2+vbgNbGmPszbPMDMMkY85vt9UrgCWPMlkyONxoYbXvZFNhVyG+hOKsBnHJ0EA6mn4F+BqCfAUBDY0wufiE75gohCqiT4bUfEJ2HbQAwxswEZgKIyObcDvdampT19w/6GYB+BqCfAVifQW73KfiGwTn7A2ggIgEi4goMBRZfsc1i4HaxXAecLzX3D5RSqpgq8isEY0yKiIwDfsJqdvqxMWa3iNxjW/8BsBSryel+rGanI4s6TqWUKmsc0g/BGLMU60s/47IPMjw3wNg8HHpmPkMr6cr6+wf9DEA/A9DPAPLwGZSqKTSVUkrlnSPuISillCqGSkVCEJEeIrJPRPaLyARHx+MoIuIsIttE5HtHx+IIIvKwiOwWkV0iMk9ESv04xiLysYicEJFdGZa9LiJ/isgOEflWRKo4MMRCl9lnYFt+v+17YbeIvOao+AqbiNQRkV9FZK/tvT5oW15NRJaLyN+2xywGHP9PiU8IGYbC6Ak0AYaJSBPHRuUwDwJ7HR2EI4iIL/AAEGaMaYrVYGGoY6MqErOBHlcsWw40NcaEAH8BWfTQKjVmc8VnICKdsUY8CDHGBAFTHRBXUUkBHjXGNAauA8bavgMnACuNMQ2AlbbX2SrxCYEMQ2EYY5KA9KEwyhQR8QN6A7McHYsDuQDuIuICVCCLviuliTFmDXDmimU/G2NSbC83YPXjKbUy+wyAe4HJxphE2zYnijywImKMOZY++KcxJhbrR6Ev1vfgp7bNPgVuyulYpSEhZDXMRVnzFvAE1ojtZY4x5ijWr8DDWEOcnDfG/OzYqIqFUcCPjg7CAa4FOojIRhFZLSKtHB1QURARf6A5sBGomd5/y/aY43SFpSEhZDbaVJlqOiUifYATmQ3tUVbY6kf7AQFAbcBDRIY7NirHEpGnsaoTPnd0LA7gAlTFqkJ5HFhgG1a/1BIRT+Br4CFjTExejlEaEoLdw1yUYu2AviISiVVldoOIfObYkIpcF+CQMeakMSYZ+AZo6+CYHEZERgB9gFtN2WxbHgV8YyybsK6cs5/IoAQTkXJYyeBzY8w3tsXHRcTHtt4HyLHarDQkBHuGwijVjDFPGmP8jDH+WO//F2NMWft1fBi4TkQq2H4J3kjZvcHeAxgP9DXGxDs6Hgf5DrgBQESuBVwppYPd2f7ePwL2GmPeyLBqMTDC9nwEsCinY5X4KTSzGgrDwWGpImaM2SgiC4GtWNUk2ygDvVVFZB7QCaghIlHAc1itisoDy221JBuMMfc4LMhClsVn8DHwsa0pahIwohRfKbUDbgN2ikiEbdlTwGSsqrI7sX4wDc7pQNpTWSmlFFA6qoyUUkoVAE0ISimlAE0ISimlbDQhKKWUAjQhKKWUstGEoJRSCtCEoJRSykYTglIFwDYefVfb85dF5G1Hx6RUbpX4nspKFRPPAS+KiDfWaJN9HRyPUrmmPZWVKiAishrwBDrZxqVXqkTRKiOlCoCIBAM+QKImA1VSaUJQKp9sQwt/jjUfQ5yIdHdwSErliSYEpfJBRCpgzb3wqDFmL/AS8LxDg1Iqj/QeglJKKUCvEJRSStloQlBKKQVoQlBKKWWjCUEppRSgCUEppZSNJgSllFKAJgSllFI2mhCUUkoB8P+QWP3adOF9FgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot fit\n",
    "fig, ax = plt.subplots(1,1)\n",
    "plt.gcf().subplots_adjust(bottom=0.15)\n",
    "plt.gcf().subplots_adjust(left=0.15)\n",
    "yMin = 0.\n",
    "yMax = np.max(xHist)*1.1\n",
    "plt.xlim(xMin, xMax)\n",
    "plt.ylim(yMin, yMax)\n",
    "plt.xticks(np.arange(xMin, xMax+1, 4.0))\n",
    "xCurve = np.linspace(xMin, xMax, 100)\n",
    "yCurve = f(xCurve, parhat)*numVal*binSize\n",
    "sig_parhat = np.copy(parhat)\n",
    "sig_parhat[0] = 1.\n",
    "bkg_parhat = np.copy(parhat)\n",
    "bkg_parhat[0] = 0.\n",
    "sigCurve = f(xCurve, sig_parhat)*numVal*binSize*parhat[0]\n",
    "bkgCurve =f(xCurve, bkg_parhat)*numVal*binSize*(1.-parhat[0])\n",
    "plt.plot(xCurve, yCurve, color='black')\n",
    "plt.plot(xCurve, sigCurve, color='red')\n",
    "plt.plot(xCurve, bkgCurve, color='dodgerblue')\n",
    "\n",
    "# Plot data histogram\n",
    "binLo, binHi = bin_edges[:-1], bin_edges[1:]\n",
    "xPlot = np.array([binLo, binHi]).T.flatten()\n",
    "yPlot = np.array([xHist, xHist]).T.flatten()\n",
    "plt.plot(xPlot, yPlot)\n",
    "plt.xlabel(r\"$x$\")\n",
    "plt.ylabel(r\"number of entries\")\n",
    "y_fitval = 0.8\n",
    "delta_y_fitval = 0.08\n",
    "if fitType == 'LS':\n",
    "    plt.figtext(0.55, y_fitval, r\"Least Squares\")\n",
    "elif fitType == 'M':\n",
    "    plt.figtext(0.55, y_fitval, r\"Maximum Likelihood\")\n",
    "for i in range(len(parin)):\n",
    "    if not parfix[i]:\n",
    "        y_fitval -= delta_y_fitval\n",
    "        plt.figtext(0.55, y_fitval, parname_latex[i] + ' = ' + f'{parhat[i]:.4f}' + r'$\\pm$' + f'{sigma_parhat[i]:.4f}')\n",
    "plt.show()\n",
    "#plt.savefig(\"histFit.pdf\", format='pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffcbd460",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}