starter_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "\n",
    "# $K^* \\mu \\mu$ angular analysis - Getting started"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from iminuit import Minuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Loading the data\n",
    "Let us load the data from each bin in separate data frames, and put them in a list called `bins`.\n",
    "\n",
    "The following assumes that the csv\n",
    "file for bin i is called `toy_data_bin_i.csv` and is located in the folder `../kstarmumu_toy_data/`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>mB</th>\n",
       "      <th>ctl</th>\n",
       "      <th>ctk</th>\n",
       "      <th>phi</th>\n",
       "      <th>mkpi</th>\n",
       "      <th>eff</th>\n",
       "      <th>year</th>\n",
       "      <th>pol</th>\n",
       "      <th>id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>75.956925</td>\n",
       "      <td>-0.986918</td>\n",
       "      <td>0.116878</td>\n",
       "      <td>-2.855914</td>\n",
       "      <td>1.035573</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2012</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>121.280030</td>\n",
       "      <td>-0.657315</td>\n",
       "      <td>0.514766</td>\n",
       "      <td>0.317335</td>\n",
       "      <td>0.987207</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2012</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>90.329200</td>\n",
       "      <td>-0.840377</td>\n",
       "      <td>0.309197</td>\n",
       "      <td>0.939818</td>\n",
       "      <td>0.927182</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2012</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>84.977505</td>\n",
       "      <td>-0.278589</td>\n",
       "      <td>-0.883990</td>\n",
       "      <td>2.403850</td>\n",
       "      <td>0.800599</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2012</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>106.125464</td>\n",
       "      <td>0.292535</td>\n",
       "      <td>0.422031</td>\n",
       "      <td>-0.377728</td>\n",
       "      <td>0.925444</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2012</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0          mB       ctl       ctk       phi      mkpi  eff  year  \\\n",
       "0           0   75.956925 -0.986918  0.116878 -2.855914  1.035573  1.0  2012   \n",
       "1           1  121.280030 -0.657315  0.514766  0.317335  0.987207  1.0  2012   \n",
       "2           2   90.329200 -0.840377  0.309197  0.939818  0.927182  1.0  2012   \n",
       "3           3   84.977505 -0.278589 -0.883990  2.403850  0.800599  1.0  2012   \n",
       "4           4  106.125464  0.292535  0.422031 -0.377728  0.925444  1.0  2012   \n",
       "\n",
       "   pol  id  \n",
       "0  1.0   1  \n",
       "1  1.0   1  \n",
       "2 -1.0   1  \n",
       "3 -1.0   1  \n",
       "4  1.0   1  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_path = 'kstarmumu_toy_data/'\n",
    "files = [f'{data_path}toy_data_bin_{i}.csv' for i in range(7)]\n",
    "bins = [pd.read_csv(file) for file in files]\n",
    "\n",
    "bins[0].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Let us have a look at the distribution of the variable $cos(\\theta_l)$, which will be used below, in one of the bins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdyklEQVR4nO3df5wcdZ3n8dfbwWBkBALoCIGFcMb1ormLMBJ/7OKMBzLKw8TVuEQxZ9AY/BGXk+w+CKuyLHor6iPnqYunOVRcz2VEVtkI0SiQwXNX2CT8MCYYGCK3JGLwB0QHEAh87o+qhqLT3VPd09XTmXo/H495TFfVt6rfU5P0Z77141uKCMzMrLyeMdkBzMxscrkQmJmVnAuBmVnJuRCYmZWcC4GZWcm5EJiZldwBkx2gWUcccUQcd9xxLa374IMPctBBB7U3UBs4V/O6NZtzNce5mjORXJs3b/51RDy35sKI2K++TjzxxGjVhg0bWl63SM7VvG7N5lzNca7mTCQXsCnqfK760JCZWcm5EJiZlZwLgZlZyRVaCCQNSdouaVTSqhrLPy3p1vTrDkkPFJnHzMz2VdhVQ5J6gEuAU4GdwEZJayNiW6VNRHww0/4DwEuLymNmZrUV2SM4CRiNiB0R8SgwDCxs0P6twOUF5jEzsxoUBQ1DLWkRMBQRy9LpJcD8iFhRo+2xwI3A0RHxeI3ly4HlAH19fScODw+3lGlsbIze3t6W1i2SczWvW7M5V3OcqzkTyTU4OLg5IvprLqx3XelEv4BFwKWZ6SXA39dpex7wuTzb9X0EndOtuSK6N5tzNce5mlPUfQRF3lm8CzgmM310Oq+WxcD7C8xiZh103KprcrVbOXcvS9O2d198epGRrIEizxFsBGZLmiVpGsmH/drqRpJeBMwAflxgFjMzq6OwQhARe4EVwHrgduCKiNgq6SJJCzJNFwPDadfFzMw6rNBB5yJiHbCuat4FVdMXFpnBrN0aHfbIHurI8mEP62a+s9jMrOT2u2GozWzi8p7MrXCPZmpzj8DMrORcCMzMSs6FwMys5FwIzMxKzoXAzKzkXAjMzErOhcDMrORcCMzMSs6FwMys5FwIzMxKzoXAzKzkXAjMzErOg85Z16s1QFq94Z7BA6SZNcs9AjOzknOPwMz2Sx5Ku33cIzAzKzkXAjOzkvOhIeu4Zrv0ZlYs9wjMzEqu0B6BpCHgM0APcGlEXFyjzZ8DFwIB3BYRbysyk1k191Cs7AorBJJ6gEuAU4GdwEZJayNiW6bNbOB84FURcb+k5xWVx8zMaivy0NBJwGhE7IiIR4FhYGFVm3cDl0TE/QARcV+BeczMrAZFRDEblhYBQxGxLJ1eAsyPiBWZNlcBdwCvIjl8dGFEfK/GtpYDywH6+vpOHB4ebinT2NgYvb29La1bpLLl2rJrz4S30Tcddj9ce9ncmYd0PE9FvVzNZmrWeD9Do/2VR1H7NJur6N9bM9ufiv8nBwcHN0dEf61lk33V0AHAbGAAOBr4oaS5EfFAtlFErAHWAPT398fAwEBLbzYyMkKr6xapbLnqDQ3RjJVz97J6S+1/vnefOdDxPBX1cjWbqVnj/QyN9lceRe3TbK6if2/NbL9s/yeLLAS7gGMy00en87J2AjdFxGPAzyXdQVIYNhaYy6Y4n/xtP+/Tqa3IcwQbgdmSZkmaBiwG1la1uYqkN4CkI4AXAjsKzGRmZlUKKwQRsRdYAawHbgeuiIitki6StCBtth74jaRtwAbgryLiN0VlMjOzfRV6jiAi1gHrquZdkHkdwLnpl5mZTQLfWWxmVnIuBGZmJedCYGZWci4EZmYlN9k3lJlZDb5u3zrJPQIzs5Jzj8DMuoJ7QZPHPQIzs5JzITAzKzkXAjOzkhu3EEg6SNIz0tcvlLRA0jOLj2ZmZp2Q52TxD4E/lTQD+D7JqKJnAGcWGcwmR/aE3cq5e3ON+X73xacXGcnMCpbn0JAi4iHgTcDnI+ItwIuLjWVmZp2SqxBIegVJD6Dy52FPcZHMzKyT8hSC/wacD3w7fZ7A8STPDjAzsylg3HMEEXEDcIOkZ6fTO4C/KDqYmZl1Rp6rhl6RPkHsZ+n0f5b0+cKTmZlZR+Q5NPQ/gdOA3wBExG3AyQVmMjOzDsp1Q1lE3FM16/ECspiZ2STIcx/BPZJeCUR6I9k5JA+jNzPbbzQzqN3KuXsZKC5K18nTI3gP8H5gJrALmAe8L8/GJQ1J2i5pVNKqGsuXSvqVpFvTr2VNZDczszbI0yP444h42l3Ekl4F/EujlST1AJcApwI7gY2S1kbEtqqm34iIFU1kNjOzNsrTI/hcznnVTgJGI2JHRDwKDAMLmwlnZmbFq9sjSO8mfiXwXEnnZhYdTL47i2cC2ZPMO4H5Ndq9WdLJwB3AB2ucmLYu5weKmO3fFBG1F0ivBgZIzhF8IbPo98B3IuLOhhuWFgFDEbEsnV4CzM8eBpJ0ODAWEY9IOhs4IyJeU2Nby4HlAH19fScODw/n/wkzxsbG6O3tbWndInVTri279jz5um867H54EsM00K3Z6uWaO/OQpraT/T20w/62vyZb33R43mHN/c46YSKfFYODg5sjor/WsrqF4MkG0rER8f+afdO0R3FhRJyWTp8PEBEfr9O+B/htRDTc+/39/bFp06Zm4wAwMjLCwMBAS+sWqZtyVY8+unpLdz7NtFuz1cvV7Ait7e5l7W/7a7KtnLuXD5zZfUeyJ/JZIaluIcjzG3hI0qdIRhx9VmVmrb/cq2wEZkuaRXK10WLgbVXBjoyIe9PJBfiyVDOzjstzsvjrJMNLzAL+Frib5EO+oYjYC6wA1pN8wF+RDlp3kaQFabO/kLRV0m0k4xctbfonMDOzCcnTIzg8Ir4k6ZzMAHTjFgKAiFgHrKuad0Hm9fkkI5uamdkkyVMIHku/3yvpdOAXwGHFRTIzs07KUwg+JukQYCXJ/QMHAx8sNJWZmXVMnucRXJ2+3AMMFhvHbGryvRbWzRrdUPY5oO61pRHhh9OYmU0Bja4a2gRsJrlk9ATgzvRrHjCt8GRmZtYRdXsEEfFVAEnvBf4kvRwUSV8A/m9n4pmZWdHy3Ecwg+QEcUVvOs/MzKaAPFcNXQzcImkDIJLHVF5YZCgzM+ucPFcNfUXSd3lq5NDzIuKXxcYyM7NOqXtoSNKL0u8nAEeRDCl9D3BUOs/MzKaARj2ClcC7gdU1lgUw3qBzZma2H2h01dC70+++iczMbAprdEPZmxqtGBHfan8cazff0Wpm42l0aOgN6ffnkTyy8vp0ehD4V8CFwMxsCmh0aOgsAEnfB+ZUHiAj6Ujgso6kMzOzwuW5oeyYzFPEAHYDf1RQHjMz67A8N5RdJ2k9cHk6fQZwbXGRzMysk/LcULYiPXH8p+msNRHx7WJjmZlZp+TpEVSuEPLJYTOzKWjcQpD2Bj5BcvWQ0q+IiIMbrmhmth9r9tLruy8+vaAkxcvTI/gk8IaIuL3oMGZm1nl5rhra3WoRkDQkabukUUmrGrR7s6SQ1N/K+5iZWevy9Ag2SfoGcBXwSGXmeHcWS+oBLgFOBXYCGyWtjYhtVe2eA5wD3NRc9KnhuFXXsHLuXpbm7Ibuz91PM+tOeQrBwcBDwGsz84LxTx6fBIxGxA4AScPAQmBbVbuPkpyD+Ks8gc3MrL3yXD56VovbnkkybHXFTp56pgHw5BDXx0TENZJcCMzMJoEionED6VnAu4AXkzzIHoCIeOc46y0ChiJiWTq9BJgfESvS6WeQjF+0NCLuljQC/GVEbKqxreXAcoC+vr4Th4eHc/+AWWNjY/T29ra0blG27NpD33TY/XC+9nNnHtL09lvVTK5O69ZsztWcqZSr2f+brZjIZ9jg4ODmiKh5HjbPoaGvAT8DTgMuAs4E8pw83gUck5k+Op1X8RzgJcCIJIDnA2slLaguBhGxBlgD0N/fHwMDAznefl8jIyO0um5RlqbnCFZvyXVLB3efOdD09lvVTK5O69ZsztWcqZSr2f+brSjqMyzPVUMviIiPAA9GxFeB06k6xFPHRmC2pFmSpgGLgbWVhRGxJyKOiIjjIuI44EZgnyJgZmbFylMIHku/PyDpJcAhJDeXNRQRe4EVwHqSHsQVEbFV0kWSFrQa2MzM2itP32eNpBnAh0n+ou8FLsiz8YhYB6yrmldz3YgYyLPNsvODZsys3fJcNXRp+vKHwPHFxjEzs04b99CQpL+TdGhmeoakjxWayszMOibPOYLXRcQDlYmIuB94fWGJzMyso/IUgh5JB1YmJE0HDmzQ3szM9iN5ThZ/neQpZV9Jp88CvlpcJDMz66Q8J4s/Iek24JR01kcjYn2xsczMrFPyPqHse8D3Cs5iZmaToPvu7TYz2w/tz080y3Oy2MzMprC6hUDSden3T3QujpmZdVqjQ0NHSnolsCB9qIyyCyPi5kKTmZlZRzQqBBcAHyEZPvp/VC0L4DVFhTIzs86pWwgi4krgSkkfiYiPdjCTmZl1UJ77CD6aDht9cjprJCKuLjaWmZl1Sp5B5z4OnEPy0PltwDmS/q7oYGZm1hl57iM4HZgXEU8ASPoqcAvw10UGMzOzzsh7H8GhmdfFP6HZzMw6Jk+P4OPALZI2kFxCejKwqtBUZmbWMXlOFl8uaQR4WTrrvIj4ZaGpzMysY/IOOncvyfOKzcxsivFYQ2ZmJVdoIZA0JGm7pFFJ+5xXkPQeSVsk3SrpR5LmFJnHzMz21bAQSOqR9LNWNiypB7gEeB0wB3hrjQ/6f4yIuRExD/gk+w5lYWZmBWtYCCLicWC7pD9qYdsnAaMRsSMiHgWGgYVV2/9dZvIgkjGMzMysg/KcLJ4BbJX0b8CDlZkRsWCc9WYC92SmdwLzqxtJej9wLjCNKTCQXbMPpzAzm2yKaPxHuKRX15ofETeMs94iYCgilqXTS4D5EbGiTvu3AadFxDtqLFsOLAfo6+s7cXh4uGHmesbGxujt7W1p3by27NrT9Dp902H3wwWEmaBuzQXdm825mlPmXHNnNn9v7kQ+wwYHBzdHRH+tZXnuI7hB0rHA7Ii4VtKzgZ4c77sLOCYzfXQ6r55h4H/VybAGWAPQ398fAwMDOd5+XyMjI7S6bl5LW+gRrJy7l9Vbuu+pod2aC7o3m3M1p8y57j5zoOl1ivoMyzPo3LuBK4EvprNmAlfl2PZGYLakWZKmAYupuhdB0uzM5OnAnTm2a2ZmbZSn5L2f5MTvTQARcaek5423UkTslbQCWE/Sg/hyRGyVdBGwKSLWAisknQI8BtwP7HNYyMzMipWnEDwSEY9KyZMqJR1Azqt7ImIdsK5q3gWZ1+fkj2pmZkXIc0PZDZL+Gpgu6VTgm8B3io1lZmadkqcQrAJ+BWwBzib5C//DRYYyM7POyXPV0BPpw2huIjkktD3Gu+bUzMwaauWeo8uGDiogSY5CIOl04AvAXSTPI5gl6eyI+G4hiczMrKPynCxeDQxGxCiApP8AXAO4EJiZTQF5zhH8vlIEUjuA3xeUx8zMOqxuj0DSm9KXmyStA64gOUfwFpKbxczMbApodGjoDZnXu4HKmEO/AqYXlsjMzDqqbiGIiLM6GcTMzCZHnquGZgEfAI7Lts8xDPV+z0NKm1kZ5Llq6CrgSyR3Ez9RaBozM+u4PIXgDxHx2cKTmJnZpMhTCD4j6W+A7wOPVGZGxM2FpTIzs47JUwjmAktIHiNZOTQUTIHHSpqZWb5C8Bbg+PQB9GZmNsXkubP4p8ChBecwM7NJkqdHcCjwM0kbefo5gil/+aiZWRnkKQR/U3gKMzObNHmeR3BDJ4KYmdnkyHNn8e956hnF04BnAg9GxMFFBjMzs87I0yN4TuW1kifYLwReXmQoMzPrnDxXDT0pElcBp+VpL2lI0nZJo5JW1Vh+rqRtkn4i6TpJxzaTx8zMJi7PoaE3ZSafAfQDf8ixXg9wCXAqsBPYKGltRGzLNLsF6I+IhyS9F/gkcEYT+c3MbILyXDWUfS7BXuBuksND4zkJGI2IHQCShtP1niwEEbEh0/5G4O05tmtmZm2U5xxBq88lmAnck5neCcxv0P5d+DnIZmYdp4iovUC6oMF6EREfbbhhaREwFBHL0uklwPyIWFGj7duBFcCrI+KRGsuXA8sB+vr6ThweHm701nWNjY3R29ubu/2WXXtaep9m9U2H3Q935K2a0q25oHuzOVdznKs5sw7paeozLGtwcHBzRPTXWtaoR/BgjXkHkfzlfjjQsBAAu4BjMtNHp/OeRtIpwIeoUwQAImINsAagv78/BgYGxnnr2kZGRmhm3aUdejDNyrl7Wb0lz1G6zurWXNC92ZyrOc7VnMuGDmrqMyyvRo+qXF15Lek5wDnAWcAwsLreehkbgdnpE852AYuBt2UbSHop8EWSnsN9Tac3M7MJa3j5qKTDJH0M+AlJ0TghIs7L86EdEXtJDvesB24HroiIrZIuklQZp+hTQC/wTUm3Slo7kR/GzMyaV7dHIOlTwJtIDsnMjYixZjceEeuAdVXzLsi8PqXZbZqZWXs16hGsBI4CPgz8QtLv0q/fS/pdZ+KZmVnRGp0jaOquYzMz2z9132nxAm3ZtadjVwKZme0v/Fe/mVnJuRCYmZWcC4GZWcm5EJiZlZwLgZlZybkQmJmVnAuBmVnJuRCYmZWcC4GZWcm5EJiZlZwLgZlZybkQmJmVnAuBmVnJuRCYmZWcC4GZWcm5EJiZlZwLgZlZybkQmJmVXKGFQNKQpO2SRiWtqrH8ZEk3S9oraVGRWczMrLbCCoGkHuAS4HXAHOCtkuZUNft3YCnwj0XlMDOzxop8eP1JwGhE7ACQNAwsBLZVGkTE3emyJwrMYWZmDSgiitlwcqhnKCKWpdNLgPkRsaJG28uAqyPiyjrbWg4sB+jr6ztxeHi4pUz3/XYPux9uadVC9U3HuZrUrdmcqznO1ZxZh/TQ29vb0rqDg4ObI6K/1rIiewRtExFrgDUA/f39MTAw0NJ2Pvf1f2b1lu77kVfO3etcTerWbM7VHOdqzmVDB9Hq518jRZ4s3gUck5k+Op1nZmZdpMhCsBGYLWmWpGnAYmBtge9nZmYtKKwQRMReYAWwHrgduCIitkq6SNICAEkvk7QTeAvwRUlbi8pjZma1FXoQLCLWAeuq5l2Qeb2R5JCRmZlNEt9ZbGZWci4EZmYl50JgZlZyLgRmZiXnQmBmVnIuBGZmJedCYGZWci4EZmYl50JgZlZyLgRmZiXnQmBmVnIuBGZmJedCYGZWci4EZmYl50JgZlZyLgRmZiXnQmBmVnIuBGZmJedCYGZWci4EZmYl50JgZlZyhRYCSUOStksalbSqxvIDJX0jXX6TpOOKzGNmZvsqrBBI6gEuAV4HzAHeKmlOVbN3AfdHxAuATwOfKCqPmZnVVmSP4CRgNCJ2RMSjwDCwsKrNQuCr6esrgf8iSQVmMjOzKoqIYjYsLQKGImJZOr0EmB8RKzJtfpq22ZlO35W2+XXVtpYDy9PJPwa2txjrCODX47bqPOdqXrdmc67mOFdzJpLr2Ih4bq0FB7Sep3MiYg2wZqLbkbQpIvrbEKmtnKt53ZrNuZrjXM0pKleRh4Z2Acdkpo9O59VsI+kA4BDgNwVmMjOzKkUWgo3AbEmzJE0DFgNrq9qsBd6Rvl4EXB9FHasyM7OaCjs0FBF7Ja0A1gM9wJcjYquki4BNEbEW+BLwNUmjwG9JikWRJnx4qSDO1bxuzeZczXGu5hSSq7CTxWZmtn/wncVmZiXnQmBmVnJTrhBIeoukrZKekFT3Mqt6w1+kJ7dvSud/Iz3R3Y5ch0n6gaQ70+8zarQZlHRr5usPkt6YLrtM0s8zy+Z1Klfa7vHMe6/NzJ/M/TVP0o/T3/dPJJ2RWdbW/TWR4VIknZ/O3y7ptInkaCHXuZK2pfvnOknHZpbV/J12MNtSSb/KZFiWWfaO9Hd/p6R3VK9bcK5PZzLdIemBzLJC9pmkL0u6T8m9VbWWS9Jn08w/kXRCZtnE91VETKkv4D+S3HQ2AvTXadMD3AUcD0wDbgPmpMuuABanr78AvLdNuT4JrEpfrwI+MU77w0hOoD87nb4MWFTA/sqVCxirM3/S9hfwQmB2+voo4F7g0Hbvr0b/XjJt3gd8IX29GPhG+npO2v5AYFa6nZ4O5hrM/Bt6byVXo99pB7MtBf6+xrqHATvS7zPS1zM6lauq/QdILnQpdJ8BJwMnAD+ts/z1wHcBAS8HbmrnvppyPYKIuD0ixrvzuObwF5IEvIZkuAtIhr94Y5uiZYfTyLPdRcB3I+KhNr1/Pc3metJk76+IuCMi7kxf/wK4D6h55+QETWS4lIXAcEQ8EhE/B0bT7XUkV0RsyPwbupHkfp5OyLPP6jkN+EFE/DYi7gd+AAxNUq63Ape36b3riogfkvzhV89C4B8icSNwqKQjadO+mnKFIKeZwD2Z6Z3pvMOBByJib9X8duiLiHvT178E+sZpv5h9/wH+97Rb+GlJB3Y417MkbZJ0Y+VwFV20vySdRPIX3l2Z2e3aX/X+vdRsk+6PPST7J8+6RebKehfJX5UVtX6n7ZI325vT39GVkio3oHbFPksPo80Crs/MLnKfNVIvd1v21X4xxEQ1SdcCz6+x6EMR8c+dzlPRKFd2IiJCUt3rdtNKP5fkHoyK80k+EKeRXEt8HnBRB3MdGxG7JB0PXC9pC8mHXcvavL++BrwjIp5IZ7e8v6YiSW8H+oFXZ2bv8zuNiLtqb6EQ3wEuj4hHJJ1N0qN6TQfffzyLgSsj4vHMvMneZ4XYLwtBRJwywU3UG/7iNyRdrgPSv+pqDYvRUi5JuyUdGRH3ph9c9zXY1J8D346IxzLbrvx1/IikrwB/2clcEbEr/b5D0gjwUuCfmOT9Jelg4BqSPwJuzGy75f1VQzPDpezU04dLybNukbmQdApJcX11RDxSmV/nd9quD7Vxs0VEdjiZS0nOC1XWHahad6RTuTIWA+/Pzih4nzVSL3db9lVZDw3VHP4ikrMvG0iOz0My/EW7ehjZ4TTG2+4+xyXTD8PKcfk3AjWvLigil6QZlUMrko4AXgVsm+z9lf7uvk1y7PTKqmXt3F8TGS5lLbBYyVVFs4DZwL9NIEtTuSS9FPgisCAi7svMr/k7bVOuvNmOzEwuAG5PX68HXptmnAG8lqf3jgvNlWZ7EcnJ1x9n5hW9zxpZC/zX9OqhlwN70j922rOvijgDPplfwJ+RHCd7BNgNrE/nHwWsy7R7PXAHSTX/UGb+8ST/UUeBbwIHtinX4cB1wJ3AtcBh6fx+4NJMu+NIqvwzqta/HthC8oH2f4DeTuUCXpm+923p93d1w/4C3g48Btya+ZpXxP6q9e+F5FDTgvT1s9KffzTdH8dn1v1Qut524HVt/vc+Xq5r0/8Hlf2zdrzfaQezfRzYmmbYALwos+470305CpzVyVzp9IXAxVXrFbbPSP7wuzf997yT5HzOe4D3pMtF8qCvu9L37s+sO+F95SEmzMxKrqyHhszMLOVCYGZWci4EZmYl50JgZlZyLgRmZiXnQmBmVnIuBGZtJGm6pBsk9aTTPZI+o2So7C3p0ATZ9tMk/TC9G9lsUrgQmLXXO4FvxVPj05wP7IiIFwOfJRmu+kmRjIB5HXAGZpPEhcAsQ9JRkv5J0i2SfibpJEkvknS9koeRXJsOL1B5IMjmdPTMH6WbOJN0OAxJBwF/FhGfSZf9HHhBjbe9Kl3PbFL4zmKzVHp4ZjPJsANXS3o28EzgX4EzI+JWSecBBwMXAzeRDGnxqKRDgYeAf4+I56fbWwj8b+AX6VscBlwbEe+set8e4JcRUcSzFMzG5R6B2VPeCNweEVcDRPJAlyHgRxFxa9pmG/A84HFgOrBaUn9EPAAcATyQ2d484IKImBcR84Dvk4z3g6TLKo3Sw0iPSnpOIT+V2ThcCMyeMo/kKV5Zc0gG+aqYSzLy6kPAS4B/AdZIeh/wMMngcxUzSHoJld7Ga4HvpD2NB6ve50DgD+35Mcya40Jg9pRfAi+uTEh6LslIsHPS6eOBJcA/SJodEQ9GxDBwNfCsSB4V2COpUgzuIHm+LMAHgWsieVzlCcDNmfc5HPh1ZJ4/YdZJLgRmT7kM6Esv9bwVeAXJk8+OSp/INgy8M5IHqnxI0nZJN5M8zvDz6Ta+D/xJ+vpy4ARJo8B/As5N57+MZFz8ikGSh+uYTQqfLDZrI0knAB+MiCUN2lwKnF25xFTSt4BVEXFHh2KaPY17BGZtFBE3AxsqN5TVabMsUwSmAVe5CNhkco/AzKzk3CMwMys5FwIzs5JzITAzKzkXAjOzknMhMDMrORcCM7OScyEwMys5FwIzs5L7/075JgtFb1GHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(bins[3]['ctl'], bins=25, density=True)\n",
    "plt.xlabel(r'$cos(\\theta_l)$')\n",
    "plt.ylabel(r'Number of candidates')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Fitting\n",
    "Let us assume a flat acceptance distribution for the following.\n",
    "We have\n",
    "$\\frac{d^2\\Gamma_P}{dq^2dcos\\theta_l} = \\frac{3}{8}[\\frac{3}{2} - \\frac{1}{2}F_L + \\frac{1}{2}cos2\\theta_l(1-3F_L) + \\frac{8}{3}A_{FB}cos\\theta_l]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def d2gamma_p_d2q2_dcostheta(fl, afb, cos_theta_l):\n",
    "    \"\"\"\n",
    "    Returns the pdf defined above\n",
    "    :param fl: f_l observable\n",
    "    :param afb: a_fb observable\n",
    "    :param cos_theta_l: cos(theta_l)\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    ctl = cos_theta_l\n",
    "    c2tl = 2 * ctl ** 2 - 1\n",
    "    acceptance = 0.5  # acceptance \"function\"\n",
    "    scalar_array = 3/8 * (3/2 - 1/2 * fl + 1/2 * c2tl * (1 - 3 * fl) + 8/3 * afb * ctl) * acceptance\n",
    "    normalised_scalar_array = scalar_array * 2  # normalising scalar array to account for the non-unity acceptance function\n",
    "    return normalised_scalar_array\n",
    "\n",
    "def log_likelihood(fl, afb, _bin):\n",
    "    \"\"\"\n",
    "    Returns the negative log-likelihood of the pdf defined above\n",
    "    :param fl: f_l observable\n",
    "    :param afb: a_fb observable\n",
    "    :param _bin: number of the bin to fit\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    _bin = bins[int(_bin)]\n",
    "    ctl = _bin['ctl']\n",
    "    normalised_scalar_array = d2gamma_p_d2q2_dcostheta(fl=fl, afb=afb, cos_theta_l=ctl)\n",
    "    return - np.sum(np.log(normalised_scalar_array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Let us check how the negative log-likelihood function defined above looks like for a given bin and a fixed $A_{FB}$ or a fixed $F_L$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mark/.local/lib/python3.6/site-packages/pandas/core/series.py:726: RuntimeWarning: invalid value encountered in log\n",
      "  result = getattr(ufunc, method)(*inputs, **kwargs)\n",
      "/home/mark/.local/lib/python3.6/site-packages/pandas/core/series.py:726: RuntimeWarning: invalid value encountered in log\n",
      "  result = getattr(ufunc, method)(*inputs, **kwargs)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAADQCAYAAAAK/RswAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABEM0lEQVR4nO3deXhU5fXA8e/JvickQIAQ9n0XEXBDcEOtFbXaaq1ad63a2lb9aWtrW2trq9bW1rVq1dbivmDFBZeIGyL7voQ9gZBA9oTs5/fH3MAACSRkZu7M5HyeZ55M3nvnznkTcjlz7/ueV1QVY4wxxphwEuF2AMYYY4wxvmYJjjHGGGPCjiU4xhhjjAk7luAYY4wxJuxYgmOMMcaYsGMJjjHGGGPCjiU4xhhjjAk7luCYIyIiUSKyWkTWB0Es6SLyhohUicgWEfn+IfatPODRKCJ/D2S8xpj2EZEuIqIH/O3mi4j48T3bfF45kv2N/0W5HYAJWdcD3YE0EUlU1SoXY3kEqAMygXHAOyKyVFVXHrijqiY1PxeRJKAAeCVAcRpjjsw4oFhVMwL4nm0+rxzh/sbP7AqOaTcRSQXuBm4CGoFRLsaSCHwH+JWqVqrq58As4NI2vPw7QCHwmR9DNMZ03DhgVaDerL3nlQ6eh4yfWIJjjsSvgI2qOhNYB4zxxUFF5H8iUtrK43+tvGwI0KCq67zalgIj2/CWlwPPq61XYkywO4ojTHACdF7pyHnI+IndojLtIiIDgJuB6U7TSrwSHBF5GcgC6oFtqnqpVxtADXCWqtYeeGxVPfsIQkoCyg9oKwOSD9OPvsBJwFVH8J7GmMAaBwwSke85329T1dEi8jpwxwGJxX4CdF45ovOQ8S9LcEx7/Rn4UFVznO9XAqd6bR8CTFDVhpbaRORNPJ9qFvkonkog5YC2FKDiMK+7FPhcVTf5KA5jjB+ISCwwHDhWVRccsHkgsMEPb9ve88qRnoeMH1mCY9pMRE7Ac5+5XEQKnOZ4oMnZHgNEeCc33m3OferuQG4rx38XOLGVt/9MVc9soX0dECUig1W1eUbXWDyJ16FcBtx3mH2MMe4bheccs8y7UUTi8NwWajzUiwN0XjnS85DxI7HhB6YtnOmYXwMLgd96beoLzHO+pgEfse+P+jt4bk01tw0Dfq2qT/o4thcBBa7Gcyl7NnBca7MXROQ4YA7QQ1XtE5YxQUxErgZuUNWjD2gfB9ymqpf46X3be15p1/7G/2yQsWmrS/BMf7xNVQuaH8B8PJdhx+D5pPUPVZ3qPHZ7t+FJcG7xQ2w/wnMlqRCYiedkuPekIiLvisgvvPa/HHjdkhtjQsI44MBbU+A5t/hzZlV7zyuH3N8Enl3BMT4jIn8AvlHVN1pqE5GBeGYtHe9akMaYsCAifwS+VtU33Y7FBCdLcIzPiMgs4KequuGAtiygFE/NnFtVdVnLRzDGmLZxzi092Dd76UJVLXExJBNkLMExxhhjTNixMTjGGGOMCTuW4BhjQp6IPCMihSKywqttnIjME5ElIrJARCY67SIiD4tIrogsE5HxXq+5XETWO4/L3eiLMcY3LMExxoSDZ4EzDmj7M/BbVR0H/Nr5HuBMYLDzuBZ4DDyrQeNZY20SMBG4W0S6+DtwY4x/hHWhv65du2q/fv0Ou19VVRWJiYn+DygArC/BqTP2ZeHChbtUtVsAQkJV54pIvwOb2VddNhXY7jyfwb41yOaJSJqI9ASmAnNUtRhARObgSZpmHuq97TwT2qwvwak9fWntXBPWCU6/fv1YsKCl8gn7y8nJYerUqf4PKACsL8GpM/ZFRLb4P5pDugV4X0QewHO1+jinPQvY5rVfntPWWvtBRORaPFd/yMzM5IEHHjhsMJWVlSQlJbWvB0HK+hKcOmtfpk2b1uK5JqwTHGNMp3YDnrIFr4nId4Gn2X/dtCPmVON+EmDChAnaloSvMya5ocD6Epx80Rcbg2OMCVeXA687z1/BM64GIB/I9tqvt9PWWrsxJgRZgmOMCVfbgZOc5ycDzYsgzgIuc2ZTTQbKVHUH8D5wuoh0cQYXn+60GWNCkN2iMsaEPBGZiWeQcFcRycMzG+oa4G8iEgXU4IyZwbMI4ll4VrWvBq4AUNViEbkH+MbZ73fNA46NMaHHEpwA+8ucdTQ2NXHb9GFuh2JM2FDVi1vZdPSBDc7sqRtbOc4zwDM+DM2YsFFV28AFj3/FfeePZmx2mtvhHJbdogqwhz9azyOfbDj8jsYYY0wQWby1lNU7yvnTe2vcDqVNLMExxhhjzGHVNTYCEBsVGqlDaEQZhuoamtwOwRhjjGmz2nrP/1uxUZEuR9I2luAEkPfK7YUVNS5GYowxxrRPrfPBPMau4JgD1Xpdtdm6u9rFSIwxxpj2qW2wW1R7iUi2iHwiIqtEZKWI/MRpTxeROc6KvXOaF7QL91V+y2vq9z7fUFTpYiTGGGNM++ypcxKcaEtwABqAn6vqCGAycKOIjADuAD5S1cHAR873EOar/Jbvadj7fENRlYuRGGOMMe1T5SQ4MZE2BgdV3aGqi5znFcBqPIvXzQCec3Z7DjjXeb53lV9VnQc0r/I7HWeVX1UtAZpX+Q0pFXYFxxhjTIjyvgsRCgJW6E9E+gFHAV8DmU5pdIACINN57vNVfnNycg4bW2VlZZv266jlRZ4rON3ihZXbdvvlPQPVl0CwvgSncOqLMabtSqrqAKhvDI1ZwAFJcEQkCXgNuEVVy0Vk7zZVVRHRVl/cTsG8ym/lsu2wcDEnDu/F64vyGXPMcaQnxvj0PWw12eBkfTHGhLriKs8VnFBJcPw+UkhEovEkNy+oavPKvjudW084Xwud9rBe5bd5DM6Uwd0AWLqt1MVojDHGmLYrrqoFoM4SHM+sKOBpYLWq/sVr0yygeSbU5cBbXu1hu8pv8xic4wZmECGw2BIcY4wxIaKkuvkKjs9uuviVv29RHQ9cCiwXkSVO2y+A+4CXReQqYAvwXWdbWK/yW15TT4RAt+RYhmQms3hridshGWOMMW2yu9JzBac+RCrx+zXBUdXPAWll8ykt7B/Wq/xW1DSQHBeNiDCxfzqvLMijtqExZMpeG2OM6Zyq6xoor/EMs7AxOOYgVbWNJMV6csopg7uxp76RhZvtKo4xHSUiz4hIoYisOKD9ZhFZ4xQa/bNX+51OQdG1IjLdq/0Mpy1XRO7AGANAQdm+5YVsDI45SHVdAwkxnqs1xw7MIDpS+HR9kctRGRMWnuWA2lgiMg1Pba2xqjoSeMBpHwFcBIx0XvOoiESKSCTwCJ6CoyOAi519jen0dnglOHYFxxykqq6RBOcKTmJsFEf37cLcdbtcjsqY0Keqc4EDx+XdANynqrXOPs2zNWcAL6pqrapuwjPmb6LzyFXVjapaB7zo7GtMp9ec4GSlxdsgY3Ow6toGEmP2jbc5aUh3/vTeGvJL95CVFu9iZMaEpSHAiSJyL1AD3Kqq3+ApEjrPaz/vwqEHFhSd1NKBg7mgaCBYX4KTP/vy5QZPkb/UyFp2ldT6/Wfmi75YghNAVXWNpCXsK+x31uge/Om9NbyzbDvXThnoYmTGhKUoIB3POnjH4Jm5OcAXBw7mgqKBYH0JTv7sy5yS5aQnFjCod1eW5ZX6/Wfmi77YLaoAqq5rIDF23xWcvhmJjO2dyttLdxziVcaYI5QHvO6sbTcfaAK6EuYFRY3xh4KyGnqkxJEUF0VlbcPhXxAELMEJoKraRhJi9r9odvaYXizPL2PTLltd3BgfexOYBiAiQ4AYYBeegqIXiUisiPQHBgPz8dTZGiwi/UUkBs9A5FluBG5MsNlaXE1Wl3iSY6OoqLEExxyguq6BpNj9a96cM64XkRHCi/O3uhSVMaFPRGYCXwFDRSTPKSL6DDDAmTr+InC5czVnJfAysAp4D7hRVRtVtQG4CU+V9NXAy86+xnRqTU3KluJq+ndNJCk2itqGJupCoNifjcEJkNzCCqrrGslOT9ivPTMljtNHZPLygm389LQhxEVb0T9j2ktVL25l0w9a2f9e4N4W2mfjqahujHFsL9tDXUMT/TISqW1oBKCqtoGYKN8uFu1rdgUnQL7e5JnBOm1o94O2XTq5LyXV9cxebmNxjDHGBJfNu6oB6Nc1YW+x2lAYh2MJToAUV3qm2GWmxB207diBGQzolshzX23Bs1qFMcYYExw27faMEe3fNZHkOE+CU+4sHh3M/L2a+EHl00XkJRFZ4jw2Ny/CKSL9RGSP17bHvV5ztIgsd8qnP+ysUh5SiqvrSI6NIibq4B+5iHDF8f1Zuq2ULzfsdiE6Y4wxpmWbiqqIi44gMzlu70SZPXWNLkd1eP6+gvMsB5RPV9Xvqeo4VR0HvAa87rV5Q/M2Vb3eq/0x4Bo8sx0GH3jMUFBSVUeXxNbvV353Qm96pMTxt4/WBzAqY4wx5tDWFJQzNDOZiAgh3ilWu6e+kyc4rZRPB8C5CvNdYOahjiEiPYEUVZ3nrDb+PHCuj0P1u+Lq+kMmOLFRkVx30gDmbypm3ka7imOMMcZ9qsrqHeUM75kCQLwzEaY6BK7guDmL6kRgp6p6X7LoLyKLgXLgLlX9DE8J9TyvfbzLqh8kWEuoby3YQ0qMHPJ9shqVtFjhly/N567JcRzJnTgrOx6crC/GmFC0s7yWkur6fQmOcwWnJgSu4LiZ4FzM/ldvdgB9VHW3iBwNvCkiI9t70GAtoV4/72MGZaczdeq4Q+5XmrqN219dRnmXIcwY12oe1yorOx6crC/GmFC0ekc5wN4EJyEmdK7guDKLSkSigPOBl5rbnJV9dzvPFwIb8CyWl4+nZHqzkCyfXlp96DE4zS4Y35uRvVL407trQmIQlzHGmPC1yklwhvVMBvbdogqF/5/cmiZ+KrBGVffeehKRbiIS6TwfgGcw8UZV3QGUi8hkZ9zOZcBbbgR9pGrqG6mqayS9DQlORITw67NHsL2shoc+XBeA6IwxxpiWrdpRTnZ6PClx0QA2yLhZK+XTwbPGy4GDi6cAy5xp468C16tq8wDlHwFPAbl4ruy868+4fa202lMvoEtC26o+ThqQwcUT+/DUZxtZuq3Uj5EZY4wxrVuZX8YI5/YUQExkBBESGldw/DoGp7Xy6ar6wxbaXsMzbbyl/RcAo3waXAAVVdQC0CUhus2vufOsYXy8Zie3v7qMt28+ocX6OcYYY4y/7KqsZfPuai6e2Gdvm4iQEBNlY3CMx2e5RQCMyU5r82tS4qL5w3mjWbuzgj+/t8ZPkRljjDEtW7SlBICj+3bZrz0+JpLK2k5eydh4rN5RQd+MBLLS4tv1ulOGZ3LZsX156vNNfLR6p5+iM8YYYw62cGsJMZERjMpK3a99YLdE1hRUuBRV21mCEwAlVXVtGmDckl+cNZwRPVP4+StLySup9nFkxoSHlpaF8dr2cxFREenqfC/Oki+5IrJMRMZ77Xu5iKx3HpcHsg/GBJtFW0oYlZVCnDNzqtm47C6s2l5OfWOTS5G1jSU4AVBSXUd6GwcYHyguOpJHLhlPY5Ny9XMLqAiBBc6MccGztLCEi4hkA6cDW72az2Tfsi/X4lkKBhFJB+4GJgETgbtFZP9r88Z0ErUNjSzLK2N8n4P/BLLS4mhoUkqq6lyIrO0swQmA0up60o4wwQHPCq6PXjKe9YWV/HjmYhqbbMVxY7wdYlmYh4DbAe8/mhnA8+oxD0hzloSZDsxR1WJVLQHmEILr3hnjCwu3lFDb0MSxAzMO2paRFAvA7iBPcNysZNxpFFfVtWsGVUtOHNyN354zkrveXMFv317Jb88ZeURLORjTWYjIDCBfVZce8LeSBWzz+r55+ZfW2ls6dlAuCRMo1pfg5Mu+vLqujkiBuvxV5Oxcvd+2rcWeGVSffPkNOzMiW3p5h/miL5bg+FlJVR176hvplhzb4WP9YHJfthZX8+TcjSTGRnH79KGW5BjTAhFJAH6B5/aUzwXrkjCBYn0JTr7sy0MrPmd83wjOPPW4g7b1LqzgvvlzyRo4jKlHsKRQW/iiL3aLys++3OBZGXxCv3SfHO/OM4dxyaQ+PJazgX98nOuTYxoThgYC/YGlIrIZzxIvi0SkB56lXrK99m1e/qW1dmM6ldLqOpbll3H8oK4tbm8uWmtjcDq5NQXlRAiMPmCa3ZESEe6ZMYrzx2fx4Jx1PPjBWlQ7z5icDUWV9LvjHb7I3eV2KCaIqepyVe2uqv1UtR+e203jVbUAmAVc5symmgyUOUvCvA+cLiJdnMHFpzttxnQqc9fvQhVOHNxygpMU57n5UxXkxf7sFpWfbSuupldavE8rEUdECPdfMJbYqAj+/nEuu6vquGfGKCIjwv921adrPUUT31m+o9VPF6bzcZaFmQp0FZE84G5VfbqV3WcDZ+FZ+qUauAJAVYtF5B7gG2e/33ktF2NMp/HBygK6JsUwLrvlSYSxUZFERwqVtQ0Bjqx9LMHxs63F1WR3SfD5cSMjhD+cN5r0xBge+WQDJVV1/OW743z+PsGmbI9nmnxynP3TNfu0tiyM1/Z+Xs8VuLGV/Z4BnvFpcMaEkNqGRnLWFnH2mJ6H/NCcGBtFVZAnOHaLyo8qaxtYsb2c4V4LlfmSiHDb9GH86uwRvLeygAse/5Jde4K78FJH5RZWAgT9H5YxxoSieRuLqaxt4PSRmYfcLzEmKuiv4Ph7NfGDqouKyG9EJF9EljiPs7y23elUF10rItO92s9w2nJF5A5/xuxLX2/cTV1DE6eO6O7X97nqhP48c/kxbC2u5rdf7uErZ2BzuGlqUuZt9PStsLzW5WiMMSb8fLCygISYSI4beOghAEl2Bafl6qLAQ6o6znnMBhCREcBFwEjnNY+KSKSIRAKP4Kk+OgK42Nk3qDU2Ke+vLCBCYFw7Ftk8UtOGdeetG48nKUa45Kl5/OWDtTQEeRnt9lq8rXRvYanSPVbR2RhjfKmhsYkPVu3kpCHdDlqe4UBJcZ38Cs4hqou2ZAbwoqrWquomPAMAJzqPXFXdqKp1wIvOvkHtsZxcXl6Qx4S+6STEBGa8yIBuSfz62HjOO6o3D3+cy4VPfMWW3VUBee9AePmbbcRFRzB5QDrlluAYY4xPfbFhN0UVtZwzttdh902MjaKy1mZRteQmEbkMWAD83CmLngXM89rHu4rogdVFJ7V24GCoMKqqPPfZHtLjhB8OrAlolczGmiq+3V3oPjaWZ1eWctqDOZw7KJrT+0UTFWKzrLx/LxV1yuuLqjk+K4r6mnKKShtDqvqoVUs1xgS7NxfnkxIXxcnDDz+sIiMxhsVbSiivqSclrmOV+v3FjQTnMeAePGvD3AM8CFzpq4MHQ4XRRVtLKHr/Sx68cCxnHd3bZ8dti+a+TAUuK9vDr99aycurdrKsPI4/nj+ao1pYOC1Yef9efjNrJY26mbsuPI7/zNvKsuK8kKo+atVSjTHBrKq2gfdWFHDuUVnERh1++YWzx/TkjcX5rMgr47ggLdkR8FlUqrpTVRtVtQn4J55bUBBG1UWXbC0F4MQh7v7Se6bG88/LJvD4D46mtLqe8x79kp+8uJhtxdWuxtVeS7aV8p95W7hoYh8GdU8mxbn322SLjhpjjE+8v7KAPfWNnD++bUsvdE+OA4K72F/AExxn1d5m5wHNM6xmAReJSKyI9AcGA/PxFN0aLCL9RSQGz0DkWYGMub3WFJTTNSlm7z8At50xqgdzfjaFG6cN5L0VBZzy4Kfc+84qdlUG/0yksup6bp65iMyUOP7vjGEAJMdFowqVdcE9wM0YY0LFi/O30TcjgQl923aVPyHWc5WnOojPw369RdVSdVFgqoiMw3OLajNwHYCqrhSRl4FVQANwo6o2Ose5CU/J9EjgGVVd6c+4O2pNQQXDevin9s2RSo6L5rbpw/jB5L785YN1PPX5Jv49bwsXHdOHa6YMICst3u0QD1LboFzx7HwKymp48drJpMZ77vM2F/kr3xO8936NMSZUrCkoZ/7mYn5x1rA2L+CcENOc4ATvFRy/JjitVBdtrXw6qnovcG8L7bPxlFcPeo1NytqCCi6d3NftUFrUMzWe+y8cy/VTB/J4zgb+M28L/5m3hbNG9+T7k/owqX96UKxQXlBWw33za9hSUc0j3x/P0X33LVbavDJ7YUUtvf1QJdoYYzqTF+ZtJSYqgguPzj78zo7m2cHBXAvH6t372MrtZdQ2NDEyK7iu4BxoYLck7r9wLD89bQhPfbaJVxduY9bS7QzslsjFE/twzthedE9x5xbbh6t2cucby6mobuKJSydw2oj9K2r2TPVcbdpRWgN93IjQGGPCQ2VtA28szufsMT3pkhjT5tc1X8HZ01mv4HRGH60uJELgpCH+rV7sK73S4vn1t0dw2/Sh/G/Zdv47fyu/f2c1985ezTH90vnW6J5MH9mDHqn+T3ZWbS/noQ/XMWfVToZmJvPjMREHJTeemD2x7Cjb4/eYjDEmnL04fyuVtQ1cdmy/dr0uOjKCmKiIoB5kbAmOjy3aWsKwHimktyMTDgbxMZFcOCGbCydks35nBf9btoPZy3dw96yV3D1rJYO7J3H8oK6cMKgrE/p1IS3BN/0rr6nno9U7efmbPL7auJvkuChumz6Ua04cwJefz23xNanx0STFRrFld2jNBjPGmGBS19DEU59tYvKA9COquJ8QE8mezjrIuDNakV/GGaN6uB1GhwzOTOanpyXz09OGsG5nBR+vKeSL3F3MnL+VZ7/cDEDvLvGMzkpleM8U+mYkkJ2eQO8u8aTFxxATdfDkvKYmpbymnh1lNazbWcHaggq+3lTMkm2lNDYpfdITuPX0IVx6bL+9g4lbIyIM75nMqh3l/ui+CUEi8gxwNlCoqqOctvuBbwN1wAbgClUtdbbdCVwFNAI/VtX3nfYzgL/hmdDwlKreF+CuGBMwby3Jp6C8hj9dMOaIXp8YExV+V3BEZJiqrvF1MKGurqGJkup6eqQE34ykIzUkM5khmclcf9JAauobWby1lCXbSlmRX8by/DLeXVFw0GvioyNJiY8iUoQmhYYmz8+l0atuTWSEMCorlRtOGshJQ7txdJ8uRLSj0vLorDRe+HoLlbUNJMVanh6u2nGueRb4B/C8V9sc4E5VbRCRPwF3Av93wLp3vYAPRWSI85pHgNPwVEz/RkRmqeoq3/TGmODR1KQ8/ukGRvRMYcrgI6vZlhofTbGzPmAwOuz/DCJyG3A+cKWqrnaa80XkelV93K/RhZiSas8vOj0ptG5PtVVcdCTHDszg2IEZe9tq6hvJK6lmy+5q8kv3UFZdT9meespr6lGFCBEiIoT0xGgyEmPplhzL4Mwk+ndNbFO1zNZ8a0xPnvliE3//eD13njncF90zLuvIuUZV54pIvwPaPvD6dh5wgfN877p3wCYRaV73Dpx175x4mte9swTHhJ13VxSwoaiKv1007ohnzvbuEs/mIF7vsC0ffQcBtwBFzQ2qWiEi3wYswfGyu9KT4HQNsfE3HREXHcmg7skM6p4c0Pcd3yeNrLR4nvh0IzPGZjGiV3DPWjNt4s9zzZXAS85zn6x7Z0yoamhs4sE5axmSmcTZYw6/sGZrstMT+Gz9LlQ1KMqLHKgtCc7HwAlA8ycqRKQrcLy/ggpVzZfqQm2AcSgSEf520TguePwrPl1XZAlOePDLuUZEfomneOgLHYpu/2O6vqivm6wvwamtfZmbV8/GojpuPiqWz+Z+esTvV7O7nj31jbz9QQ4psb5NcHzxezlsgqOqLzmXjjeIyGbgA+BrYG2H3jkMLc0rBXCtfkxnM6FfOv0yElieX+p2KMYH/HGuEZEf4hl8fIqqNg8CO9T6dm1a9y4YFvV1k/UlOLWlLzX1jfzigRzGZqfxs+8e16ErL42rd/LC6gVkDx/n84WcffF7adNaVKp6P56SanfjmV1wK1DRoXcOM9tL93D/+2tJT4yhX4ZV1w2UwZnJrN9Z6XYYxkd8ea5xZkTdDpyjqt41BcJm3Ttj2uvpzzexvayG26cP7fBtpeZK8ttKgrMmWZunn6jqHjzLJcwGEJGp/gkpNK0p8ExZ/umpg4PyXmS4GpKZxCdrCqmua9hbOtyEtiM517Sy7t2dQCwwx/mbnKeq14fTunfGtMeOsj384+NcTh+RyfGDjmzmlLfsdM+M4byS4KxJ1q7VxJ3BfgCoak4b9n9GRApFZIVX2/0iskZElonIGyKS5rT3E5E9IrLEeTzu9ZqjRWS5iOSKyMMShBlEc9G5M0f3PMyexpeOG9iVhibli9zdbodifKi95xpVvVhVe6pqtKr2VtWnVXWQqmar6jjncb3X/veq6kBVHaqq73q1z1bVIc62g9bFMyaU3fvOappU+dXZI3xyvISYKJJjoygsr/XJ8XytXQkOLSyEeRjPAmcc0DYHGKWqY4B1eD5lNdvQ0skIeAy4Bs+l5MEtHNN1W3ZXkxgTSYYNMA6oY/qlkxQbxcdrdrodivEtSy6M8aGvNuzmf8t2cMPUgWSn+24YRbfkWIoqwiPBadeVE1WdCxQf0PaBqjbXdp6HZyBf628o0hNIUdV5ziDB54Fz2xNHIGzZXUXfjES7PRVgMVERTBnSlY9WF7JvDKkJA/aHZIyP1NQ3cteby+ndJZ7rTxro02MHc4LT3kELvv4fxLs2BUB/EVkMlAN3qepneOpT5Hnt412z4iBuTN9sUuXL3GpGd410fbphZ5zy2It6Civq+P0LH3Ji70Mv8+CWzvh76SDLVo3xkYc+XMeGoiqev3IicdFHXmC1Jd2SY1meX+bTY/qKa6MyW6hNsQPoo6q7ReRo4E0RGdne47oxffO1hXnUNi5l4vB+TJ067IiP4wudbcojwLENjcwt/JxnVlYyfPhwvjM+K+iupHXG34sxxn2Ltpbwz7kbuXhiNlOGdPP58bPTE3hvRQENjU1ERbb3ppB/uRKNV22KS5prU6hqrarudp4vxLM43hA8dSi8b2O1WpvCLd9s9tyFu3HaIJcj6ZxioyK5ZFJfVOHWV5Yy5K53g3p9FGOMCYSa+kZue2UpPVLi+MVZ/lnSpn/XRBqalGv/vXC/9QaDQXsTnA6P5GytNoWIdBORSOf5ADyDiTeq6g6gXEQmO7OnLgPe6mgcvvDVht1M+P2HvPjNNqYN7UaiLfromsuO7ctnt08DoL5Ref6rze4GZDrKRo0b00H3/G8VG4qq+NMFY0iO88/t+0HdkwD4eE0hG4qCqyZZuxIcVT2tPfs7tSm+AoaKSJ6IXIVnxd9kPLUpvKeDTwGWicgS4FXgelVtHqD8I+ApIBfPlZ13CQIfrt7Jrspaxman8bsZo9wOp1MTEbLTE3ji0qMB2FgUvAvAmcNr77nGGLO/d5bt4IWvt3LdlAGcONj3t6aaHZWdxrVTBgCwdXdw1cPx6yUHVb24heanW9n3NeC1VrYtAIIqg6ipb2TB5mLG9k7lrRttWa5gMX1kD04Y1JUtxcH1h2ZMICzeWkJxVR2nDM90OxTjom3F1dzx2jKO6pPGrdOH+vW9RITrTxrIk3M3cvXzC9h837f8+n7tEVwjgkJEbUMj5z36JUvzyvjWGCvsF2wGdU9ibUG5jcMxncqOsj2c9+iXXPXcAiuZ0IlV1zVw/X8WgsDDFx1FdAAG/nZJ2Hf7q2xPvd/fr60swWmnzbuquOJf37B6RzmXTOrDtVN8W1PAdNz3J/Whpr6JVxduczsUYwLmlheX7H2+rTg41wYy/qWq3PbKMlbtKOfhi4/yaUG/QxER/nXFMQCsLQieZSotwWmHnLWFTH0ghy837CZC4DY/X/ozR2ZIZjIjeqaQs7bI7VCMCZhV28v3Pl+8rcTFSIxbZm2o553lO7jzzGFMG9o9oO89vEcKsG9dxmBgCU47vLrQU2/w5pMHsfw300lLsGUZglX/bonsKKtxOwxjAuKL3F1U1DZw55nDiIuOYPHWUrdDMgE2e/kO3sit5/yjsrjmxAEBf//MlFjSEqJZvcOu4ISc2oZGFm8tZfrITH5++lCbEh7keqXGsb10j41FMJ3CJU99DcCIXikc0y+duevt6mVn8uWGXdzy4hIGpUXwh/NHu1LoVEQY2SuFmfO3smBz8eFfEACW4ByGqrIiv4yhd71Hfukezjuq1VUiTBDpkRpPbUMTK/KD53KpMf7SLTkWgEn9Mzh9RCYbi6rILQyumiTGP1bkl3Ht8wvp1zWBW8bH+XwphvYY2SsVgJ+9vNS1GLxZgnMYL32zjbP//jkAN0wdyBmjbNZUKDh1uOf+8w0vLKShscnlaIy/icgzIlIoIiu82tJFZI6IrHe+dnHaRUQeFpFcEVkmIuO9XnO5s/96Ebncjb4cqe9NyCYmKoJTR3imiM9evsPliIy/bSyq5If/mk9qfDTPXTmRpBh3l6hpXsgzJT447nBYgnMIC7eUcMfrywEYnZXK7TaoOGT0zUjk1OHdySvZw1H3zOGpzzba7arw9ixwxgFtdwAfqepg4CPne4Az8VRKH4xnYd7HwJMQAXcDk4CJwN3NSVEwU1VKqurISPKMCeyZGs/xgzJ46ZttNAVZ6XzjOxuKKrnoyXmownNXTqRnarzbIZGeGMMlk/qwrTg4hgdYgnMIT87dAMB954/m7ZtPCLoFHM2h/eackcRGRVBR08Dv31nNH99d43ZIxk9UdS5w4I3/GcBzzvPngHO92p9Xj3lAmoj0BKYDc1S1WFVLgDkcnDQFneKqOhqalPTEfZMeLjqmD/mle/gsd5eLkRl/yS2s5OIn59GkysxrJ+9dLiEYDOuZQtmeevJK3C9VEBzXkYLMhqJKvv33z6mua+TyY/ty0cQ+bodkjkDvLgms+t0ZbCyq5Ldvr+LZLzZz88mD/LYmiwk6mc5adgAFQHN53yzAu0hSntPWWvtBRORaPFd/yMzMJCcn57DBVFZWtmm/9pq9yVPQsnbnJnJytgIQ16QkR8NDby9Ej47z+Xv6qy9uCLW+bKto4oEFNajC/02MY/vqhWxf7dkWDH1pLGsE4O6Zn3HZyNgjPo4v+mIJjpemJuWmmYuYvbxgb9uVJ/R3MSLTUZERwuDMZG6dPpRzH/mC57/aYqu+d0KqqiLis2vmqvok8CTAhAkTdOrUqYd9TU5ODm3Zr73+V7SUbslF3HjBKfu1X9W0jr9+uJ7uQ8YzoleKT9/TX31xQyj1Zd7G3fz5+QUkxMbwn6smMTgzeb/twdCXpiblL4s/4ONtDfzknEmMzU47ouP4oi92i8rLXz9av19yk3vvmfTNSHQxIuMr47LTOG1EJo9+khtUpcSNX+10bj3hfC102vOBbK/9ejttrbUHtcKKWnqlHnyV5orj+5McG8XDH613ISrja+8u38Flz8yne3Isr//o+IOSm2ARESH887IJACzNK3U3Fn+/QajMbsgrqebhj9Zz4uCufHrbVN6+6QSiArCGhwmcH588mKq6Rptd0nnMAprPFZcDb3m1X+acbyYDZc6trPeB00Wki3NOOt1pC2qF5TV0Sz44wUmNj+bKE/rz3soC5m8Kjrokpv1Ulac/38SP/ruIUb1SePX648hKc39A8aFM7J9Ot+RY3lnm7rk2EP+DP0sIzG7Y4izz/qOpg+ibkcjo3qm+PLwJAqOyUkhLiGaZy58qjO+JyEzgK2CoiOSJyFXAfcBpIrIeONX5HmA2sBHIBf4J/AhAVYuBe4BvnMfvnLagpaoUlNfQPaXlsQ7XnTSAXqlx3D1rJY02oyrk1NQ3cusry7jnf6s4fUQmL1w9mS6JwV9BX0S4+oT+fL2pmE27qlyLw+8JTijMbmhS5YOVnltTvbsEd2ZsjpyIMKxHMquCqJS48Q1VvVhVe6pqtKr2VtWnVXW3qp6iqoNV9dTmZMU5v9yoqgNVdbSqLvA6zjOqOsh5/Mu9HrXNmoIKSqvrGZPV8geyhJgofvmtEazeUc7zX20ObHCmQ3aW13DRk/N4bVEePzllMI9dcjTxMe4V8WuvM0b1AOCa5xfwydrCw+ztH24NMg6a2Q0VdcrczVW8stFzBWfdkq/ZEBG608GDYRS9r/ijL92ljlnb6rngofe4dkwsMZGB+V3b78X4Wll1Pf/8bCMAJw9rfWHFs0b34KQh3fjTe2s4cXC3oJpSbFo2d10RP3t5CdV1jTz+g/EhWWC2j7OSeW5hJVf86xs23/etgMfg+iwqN2c3qCoj736f6jrPf3I3TRvEKSeHdjG/YBhF7yv+6MugsdV8+ciXLNhZS/zmBJ6/cmJA6hvZ78X42i/eWM47y3cwrEcy3VNanwouItx/wRim/3Uut7y0mNduOI7YqNC5EtCZ1DU08eAHa3li7kaGZCbx32vGMyRIBxMfjvd5Ncal8axujaINitkNjU3Kpcf2BeD2M4Zyq1UqDnu9uySw4K5T+eVZw/ls/S4+Wu3OpVNjOmp5fhkA931nzGH37Z4Sx33fGcOK/HJ+9eaKoKgya/a3eVcVFz7xFU/M3cglk/ow66YTQja5afbEpUeTFBtFoyp1DYFfMsetBCcoZjdERUZw55nDeXhaAj+aarVROpNLj+3L0MxkfvTCIv44ezWvLNhmJ30TMlSVkuo6fjC5D+PaWGdk+sge3HzyIF5ekMfTn2/yb4CmzRqblKc+28gZf5vLxqJKHr1kPPeeN9rVRTN9ZfrIHtx/wRgam9SVyR1+v0XlzG6YCnQVkTw8s6HuA152ZjpsAb7r7D4bOAvP7IZq4ArwzG4QkebZDeDj2Q0psaE75sYcmbjoSGZeO5kZj3zOE3M94xjKaxq4ygo7mhBQUF5DRU1Duz/h//TUIeQWVvKH2avpmRrPt8aE3tiOcJJbWMFtry5j8dZSThnWnXvPG02PFmoahbLJAzIQgQse/4p7zh3FpZP7Buy9/Z7gqOrFrWw65cAG9XyEvrGV4zwDPOPD0Ewnl54Yw8xrJrNoayk/fWkJf5y9mu9P7BNSMxVM57RuZyUAg7u3L8GJiBAe/O5Yip6u5ScvLiY6Ujh9ZA9/hGgOobqugUc+yeWfczeREBvJX783jhnjeoXleoddEmOIFKFBlV+9uYLzjsoiKTYww3+tkp3p1Hp3SeCcsb247/zRNDTp3nENxgSrzbuquPyZ+QAM79n+MRoJMVH864pjGJWVyo3/XeR6MbbORFWZtXQ7Jz/wKY98soGzx/Rkzk9P4tyjssIyuWn25wvG7C1OuHRbacDe1xIcY4BTh2cSGxXBvbNXU1Pf6HY4xrQqx6kpct1JA0hLOLKib8lx0Tx35UTG9k7jxv8u4p9zN9oYND9buq2Ui56cx49nLiYjKYZXrz+Wv3xvHN2Sj3xBylBx/vjezP7xiUBgl2+wBMcYPJdRH7hwLEu3ldpSDiZoNTYpM+dvIzs9njvPHN6hY6XGR/Ofqydx1uge3Dt7NXe9uYLaBkvufW3dzgqu+/cCZjzyBesLK7n3vFHMuukEJvRLdzu0gEpNiKZrUixbnVUDAsESHGMcZ4/pSb+MBH728lJO+NPHrCkodzskY/azaGsJa3dW8JNThvjkeHHRkfzj4vFcd9IAXvh6K9957Eu27HavtH442VBUyc9eWsL0v87ly9zd/Oy0Icy9fRqXTOpLZAgXk+2IrLQ4Vm4v57p/L9i7eoA/WYJjjENE+Mf3x9M3I4G8kj38+s2VbodkzH7mrisiMkKYPjLz8Du3UUSEcOeZw3ny0qPZuruasx/+nJnzt9Jka1cdkcVbS7ju3ws49S+f8s7yHVx74gDm3j6NH58yOGCDa4NVt+RYlueX8f7Kncycv9Xv79e5f9rGHGBUViqf3jaNhz9az1/mrKOgrCbspm2a0LS2oIJHczYwvGcyyXHRPj/+6SN78E7PFG59ZSl3vr6ct5bk88fzx9C/a6LP3yvcNDQ28dGaQv71xSbmbSwmJS6Km6YN4vLj+tE1KfzH2LTVSUO786FTXPWTtUUUltccsgp3R9kVHGNacNZoT32QyX/8iFcWbDvM3sb43/S/zqWxSRnbO81v75GdnsDMayZz3/mjWbm9nOkPzeXed1ZRWl3nt/cMZTvLa/jbh+s58c+fcN2/F7JldzV3fWs4X955Cj8/faglNwe4dHJfVv1uOmN7exaHvfP15X59P0twjGnBoO5J3Hyyp7r1ba8u47p/L6CqtsHlqMyREJGfishKEVkhIjNFJE5E+ovI1yKSKyIviUiMs2+s832us72fy+EDsL10z97nF0/s49f3iogQLprYhw9/dhIzxvXiqc83MeXPn/BoTi7lNfV+fe9QUNvQyHsrdnDdvxdw3H0f89CH6xjUPYknLj2az26fxtUnDuj0t6IOJSEmir98bxwAH60p9OtYR0twjGnFz08fyme3TyM7PZ73V+7k5pmLWVtQ4XZYph1EJAv4MTBBVUcBkcBFwJ+Ah1R1EFACXOW85CqgxGl/yNnPdcfd9zHgSW5GZaUG5D0zU+K4/8KxzP7xiYzv24U/v7eW4//4MX98dzUlNYFfV8hNjU3KVxt2c8dryzjm9x9y/X8WsXBLCVce349Pbp3Kv6+axPSRPYhyaVHJUDOwWxKv3XAsAGf89TO/jfeyNNOYQ8hOT+Djn0/lV2+u4MVvtjF3XRE3TB3Ixl1VXD9lIKN7B+Y/G9MhUUC8iNQDCcAO4GTg+87254DfAI8BM5znAK8C/xARUReLxJRV77tq0suF8WDDe6bw7BUTWZ5XxhNzN/DPuRsR4P2ihXzvmGxOHNw1LP9jr65r4LP1u/hw1U4+XlPI7qo6EmIiOWNkD2YclcXxAzPCst+BMrLXvnNnfukestMTfP4eluAYcxjRkRH84bzRHNMvnd/9bxV//zgXgG3F1fz+3FGM8eOYCNMxqpovIg8AW4E9wAfAQqBUVZvvOeYBWc7zLGCb89oGESkDMoBd3scVkWuBawEyMzPJyck5bCyVlZVt2u9Aq3fvq02TXLmVnJz8dh/DVy7oBVPS4nkvt5ov1hfw3soC0mKFiT0iGZ8ZxeC0iJCbAt38e2lSZWt5E6uLm1i9u5HVxY3UN0F8FIztFsl3B8UyrlsksVGl6PZSPt/uduQHO9J/Y2576I3P+faA6P2qOfuiL64kOCIyFHjJq2kA8GsgDbgGKHLaf6Gqs53X3Inn8nEj8GNV7fBq4sa0VUSE8J2je/Odo3vT1KTc8tISZi3dzjn/+IJLJvVhWM8U4qIi2FlewzVTBhAbZetZBQMR6YLnqkx/oBR4BTijo8dV1SeBJwEmTJigU6dOPexrcnJyaMt+3kqr67j1L3NJT4xhzk+nkBEkg1a75+Tw+AlT+HhNIa8uzCNnfREfbKkhPTGGk4d15/hBGUwekEHP1Hi3Q21VdV0DK7eX837OQnYXJ/H1pmLK9niulg3olsglk7tx2ohMJvZPJzpErtQcyb8xNy2dVM/Y337A6+vrOX/KeE4Y3HXvNl/0xZUER1XXAuMARCQSyAfewLN6+EOq+oD3/iIyAs9985FAL+BDERmiqlZ20wRcRITwx/NHc95RWfzm7ZW88PX+9Rwe/3Qjk/qn84/vj7eFO913KrBJVYsAROR14HggTUSinKs4vfGcg3C+ZgN5IhIFpAK7Ax+2x/+W7WBXZS2v/+i4oElumsVERXDGqB6cMaoHVbUNfLquiA9WFjBn1U5eXZgHQN+MBCb0TWdkrxRG9kphRK8Uv0xxPxRVZXtZDRsKK8ktrGT1jnKW5ZWxvrCC5qEf2enlTB+ZyXEDu3LswAwy/Th12eyTGr/v30J+qe8rHAfDLapTgA2quuUQi43NAF5U1Vpgk4jkAhOBrwIUozH7SYyNYtqw7kwb1p36xiZKquqorG3g0ZwNvLUkn4/WFHLygzlMH9mD6roGJiZ2rkGZQWQrMFlEEvDcojoFWAB8AlwAvAhcDrzl7D/L+f4rZ/vHbo6/mb+pmJ6pcRyVneZWCG2SGBvFWaN7ctbonjQ2KWsKypm3sZivNuzm03VFvLYob+++vbvE0zcjgT7pCfRJT6R3l3gykmLISIwlIymGLgkxbb7NVdvQSPmeBspr6impqmNHWQ0FZTWer+V72FpczcaiKqrr9n0WTk+MYUzvVKaP6sHY3qlUbF3JudNP9vnPxLTN3NumMeX+T9heWuPzYwdDgnMRMNPr+5tE5DI8J6Gfq2oJnvvi87z28b5nvp9A3hsPRtYXd53dDc46NZ4X19TxxfYa/vPVZhoUXgaWFM7h1L6B/fTqD6H0e1HVr0XkVWAR0AAsxnNr6R3gRRH5vdP2tPOSp4F/Ox+iivGcnwJOVXngg7XMWrqd4wdlhNRK05ERwsheqYzslcpVJ/QHoLC8hpXby1m5vYz1hZVs2V3N+yt3Ulx1cH0dEYiLiiQ2OoLYqAhioyKJihAampTGJqWhqYmGRqWytoHahpY/OCTERNIjNY6stHi+d0w6A7slMah7EgO7JdE1KWa/n2fOztX++UGYNumTkUBMVAR/+2g9w3okc6ZTg8wXXE1wnNoT5wB3Ok2PAfcA6nx9ELiyPccM1L3xYGV9CQ4nT9v3/N9fbeZXb63k9Q2N/PqSU4iJCo37+a0Jtd+Lqt4N3H1A80Y8V4EP3LcGuDAQcR3Kkm2lPPLJBpLjorh0cl+3w+mw7ilxdE+JY9qw7vu1V9TUs720ht1VtRRX1VFcVcfuyjr21DdSW99IbUMTNfWN1Dcp0RFCZEQE0ZFCZISQFBtFSnw0KXGer2kJMfRMjaNHahzJsVEhlRR2dtdPGcDDH+fy2qK88ElwgDOBRaq6E6D5K4CI/BP4n/Nt833xZt73zI0Japce24+SvFz+srCWbz38GZERwn+unmRVTk2rVuSXAfDeLVPISgvegbodlRwXzdAe0UCy26EYF/3s9KG8vWxHq1fkjpTbHyUvxuv2lIh4p27nASuc57OAi5wqo/2BwcD8gEVpTAeN7hrJD4/rR31jE2sKKvjD7NW8/M02PlhZwKyl29lYVOl2iCZIbN5VxWM5G+ieHOtK3Rtj3DC2dypLt5XuV/epo1y7giMiicBpwHVezX8WkXF4blFtbt6mqitF5GVgFZ776DfaDCoTSkSE35wzEtURzHjkC15flM/ri/ZdhOyeHMvTlx/DqKwUu7TeyT2ak8v2shpmXjPZ/i2YTuOy4/rx5pLt3PH6Mh69ZLxPjulagqOqVXgKaHm3XXqI/e8F7vV3XMb4k4gw66YTqK5roLiqjtLqehZvK+VXb67g2//4nClDunHPjJH0zbAVnDurL3J3860xPTl2YMbhdzYmTByVnUZKXBTvrijglQV5dD/8Sw7L7VtUxnRKCTFR9O6SwKisVC6d3Jf/Xj2Jwd2TmLuuiJPuz+GH/5rPR6t3UlNvFyo7k7ySavJL9zC8h41JMZ2LiPDqDccBsCy/1CfHdHuQsTEGOG5QV+b87CQWbC7mV2+tJGdtETlri0iJi+LT26bRJTHG7RCNn9XUN3LFv74BPOs/GdPZDMlMZliPZHaW13rWNeggS3CMCSIT+qUz+8cnsHhbKQ9+sJYvcncz/a+eUv0nDOrKqSMyyUqL98vCdMZdN/13EesLK7nomGymDfXFBXpjQk+35Fg27apiV7cmGpu0Q2ub2S0qY4KMiDC+TxdeuHoyD144luMGZrCrspanPt/ERU/O47tPfMWeOrt1FU6Kq+r4aE0h1580kPu+M4aIEFuw0hhfyUyJI7ewkls/3cOmXR2bXWpXcIwJYs0LfNY1NPHpuiLmrivi3/O28PriPKYO7U6v1DibaRPiKmsbmPSHD1GFb4/1XZEzY0LRxRP78M3mYgYl1pGV1rEr1ZbgGBMCYqIiOG1EJicO7srri/L45RueElHfGt2T608aSHZ6PKnx0ZbshKDXFuZR36icO64XI3uluh2OMa46um8XPr1tGjk5OR1erNgSHGNCSFx0JLN/ciLrd1Yyc/5W3lm+g3eW7wA86+/0y0jkT98Zw8DuiTQ0KSkBXrnZHN624mr++dlGzhzVk7HZqdw9ayU9U+P460VHuR2aMWHFEhxjQkzfjET6ZiRyyvDurCmoYGtxNduKq1m9o4LXFuXx7X98vnffu789guwuCZwwuCtx0R37NGQ6RlX5fP0unv1yEx+uLuT5r7YwsJun3tHxg7q6HJ0x4ccSHGNClIgwvGfKflOKrz6xP8vzyyiqqOX+99fy27dXAXD+UVlcfeIARvSy6cduWbm7iQfe/xqACX27sGBLCRuKqoiLjuD3545yOTpjwo8lOMaEEe+E5weT+7JwSzG3v7qc1xfn8/rifK6bMoCx2WkM7JZE34wEu6oTIKrK7E11e7+/48xhREVG8PKCbXx/Yh/7PRjjB5bgGBOmUuOjOXlYJvPu7MaaggpufWUpT8zduHe7CGR3SWBAt0T6ZSTSPSWW7slxpCdGk5YQw/g+XVyM3ndEJA14ChiFZ527K4G1wEtAPzzr3n1XVUvEM0r7b8BZQDXwQ1Vd1NEYdlfVsWp3E2eP6ckVx/fj6L7pAIzLTuvooY0xrbAEx5gwFxUZwaisVN67ZQpVtQ1s2lXFhqJKNhbt+7pwcwkVtQ37ve61G47d+x9xiPsb8J6qXiAiMUAC8AvgI1W9T0TuAO4A/g84ExjsPCYBjzlfO2R76R4AzhnbK1x+psYEPVcTHBHZDFQAjUCDqk4QkXQC+MnKmM4kMTaKUVmpjMo6eDrynrpGCitqyC/Zww0vLOLKZxcwvGcyvbsk0LtLPL27JNAjJY6i6iYXIj8yIpIKTAF+CKCqdUCdiMwApjq7PQfk4ElwZgDPq6oC80QkTUR6quqOjsTRnOBkdYnvyGGMMe0QDFdwpqnqLq/v7yCAn6yMMR7xMZF7Z2g9e8UxPPvlZvJL9vD5+l3srKhB1bPflN5RXOhuqO3RHygC/iUiY4GFwE+ATK+kpQDIdJ5nAdu8Xp/ntO2X4IjItcC1AJmZmeTk5BwyiE831wOwacUiitaFfq2iysrKw/Y5VFhfgpMv+hIMCc6BAvrJyhhzsKP6dOEorzE4tQ2NbC+toaiilnUrlrgXWPtFAeOBm1X1axH5G54PTXupqoqItuegqvok8CTAhAkTdOrUqYfcf1hZDZnvf863TpsaFsUYc3JyOFyfQ4X1JTj5oi9uJzgKfOCcXJ5wThoB/WQFlvUGK+tLcEqT6lDqSx6Qp6pfO9+/iifB2dn8AUlEegKFzvZ8INvr9b2dtg7pkRrHuO5RYZHcGBMq3E5wTlDVfBHpDswRkTXeGwPxyQos6w1W1pfgFEp9UdUCEdkmIkNVdS1wCrDKeVwO3Od8fct5ySzgJhF5Ec8t8DK7SmxMaHI1wVHVfOdroYi8AUwkwJ+sjDFh72bgBWcG1UbgCiACeFlErgK2AN919p2NZyJDLp7JDFcEPlxjjC+4luCISCIQoaoVzvPTgd/h+QRln6yMMT6hqkuACS1sOqWFfRW40d8xGWP8z80rOJnAG8496Sjgv6r6noh8g32yMsYYY0wHiGq7hriEFBEpwpMkHU5XYNdh9woN1pfg1Bn70ldVu/k7GLfZeSbkWV+CU3v60uK5JqwTnLYSkQWq2tIl7JBjfQlO1hcTTj8360twsr7sL8JXwRhjjDHGBAtLcIwxxhgTdizB8XjS7QB8yPoSnKwvJpx+btaX4GR98WJjcIwxxhgTduwKjjHGGGPCjiU4xhhjjAk7nTLBEZELRWSliDSJSKvT0ETkDBFZKyK5InJHa/u5SUTSRWSOiKx3vnZpZb9GEVniPGYFOs5DOdzPWURiReQlZ/vXItLPhTDbpA19+aGIFHn9Lq52I87DEZFnRKRQRFa0sl1E5GGnn8tEZHygYwx2dp6x84y/hMt5Bvx8rlHVTvcAhgNDgRxgQiv7RAIbgAFADLAUGOF27C3E+WfgDuf5HcCfWtmv0u1Yj/TnDPwIeNx5fhHwkttxd6AvPwT+4XasbejLFGA8sKKV7WcB7wICTAa+djvmYHvYeSZ4HnaeCd6HP881nfIKjqquVs/KwocyEchV1Y2qWge8CMzwf3TtNgN4znn+HHCue6Eckbb8nL37+CpwijhrfASZUPk3c1iqOhcoPsQuM4Dn1WMekOYsjmscdp4JKnaeCVL+PNd0ygSnjbKAbV7f5zltwSZT9y06WoBnja+WxInIAhGZJyLnBia0NmnLz3nvPqraAJQBGQGJrn3a+m/mO86l1ldFJDswoflcqPx9BLtQ+TnaeSZ4dKbzDHTgb8TNxTb9SkQ+BHq0sOmXqvpWC+1B61B98f5GVVVEWpv331dV80VkAPCxiCxX1Q2+jtUc1tvATFWtFZHr8HxiPNnlmMwRsvPMQew8ExzsPEMYJziqemoHD5EPeGe9vZ22gDtUX0Rkp4j0VNUdzmW7wlaOke983SgiOcBReO7juq0tP+fmffJEJApIBXYHJrx2OWxfVNU77qfwjG0IRUHz9+EmO88cdAw7z/hfZzrPQAf+RuwWVeu+AQaLSH8RicEz6CyoZgU4ZgGXO88vBw761CgiXUQk1nneFTgeWBWwCA+tLT9n7z5eAHyszuizIHPYvhxw7/gcYHUA4/OlWcBlzgyHyUCZ1y0M03Z2ngkMO8+EriM/17g9gtqNB3Aenvt4tcBO4H2nvRcw22u/s4B1eD6B/NLtuFvpSwbwEbAe+BBId9onAE85z48DluMZbb8cuMrtuA/ow0E/Z+B3wDnO8zjgFSAXmA8McDvmDvTlj8BK53fxCTDM7Zhb6cdMYAdQ7/ytXAVcD1zvbBfgEaefy2llllBnfth5xs4zLvYlJM4zTqx+O9fYUg3GGGOMCTt2i8oYY4wxYccSHGOMMcaEHUtwjDHGGBN2LMExxhhjTNixBMcYY4wxYccSHGOMMcaEHUtwjDHGGBN2LMExQUNErhORHSKyxHn8x+2YjDHhQ0T+LiJb3I7DBEbYrkVlQtJo4C5VfdrtQIwx4UVE+gHTgBgRSVbVCpdDMn5mV3BMMBkDLHE7CGNMWPot8Hs862ONdDkWEwCW4JhgMhL4l3N76kO3gzHGhAcRGQmMAl7Cs/DkKHcjMoFgt6hMUBCRbKBAVccc0P62qn7bpbCMMeHh98CvVVVFZDVeV3BEZCGeFboBbgE+B+YBw4Efq+rKAMdqfMQSHBMsRuNZ/XYv5575ZjeCMcaEBxGZBJwBHCUij+BZNXy5sy0b+EZVr/f6fr6q3iQiPwV6c8B5yYQOS3BMsBjDwSeS8cAiF2IxxoSPPwDfVtUPAUQkE1jsbDsaGC4ij+O5fZUKDHG+76mqD7kRsPENG4NjgsVoPIP/vB0NLHQhFmNMGBCRU4GY5uQGQFV3Akkiko7nHHOLql6vqp843//cuaLTICKJrgRufEJU1e0YjGmRiLwLbAEagFdU9VOXQzLGhBEReQ34nqo2ON/PxnPOaQJKVPUuN+MzHWMJjjHGGGPCjt2iMsYYY0zYsQTHGGOMMWHHEhxjjDHGhB1LcIwxxhgTdizBMcYYY0zYsQTHGGOMMWHHEhxjjDHGhB1LcIwxxhgTdv4f7YgcWjIVFPkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_test_bin = 1\n",
    "_test_afb = 0.7\n",
    "_test_fl = 0.0\n",
    "\n",
    "x = np.linspace(-1, 1, 500)\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3))\n",
    "ax1.plot(x, [log_likelihood(fl=i, afb=_test_afb, _bin=_test_bin) for i in x])\n",
    "ax1.set_title(r'$A_{FB}$ = ' + str(_test_afb))\n",
    "ax1.set_xlabel(r'$F_L$')\n",
    "ax1.set_ylabel(r'$-\\mathcal{L}$')\n",
    "ax1.grid()\n",
    "ax2.plot(x, [log_likelihood(fl=_test_fl, afb=i, _bin=_test_bin) for i in x])\n",
    "ax2.set_title(r'$F_{L}$ = ' + str(_test_fl))\n",
    "ax2.set_xlabel(r'$A_{FB}$')\n",
    "ax2.set_ylabel(r'$-\\mathcal{L}$')\n",
    "ax2.grid()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Note the warning about an inalid value. For some values of the fit parameters you can get a negative PDF which leads to invalid logarithms.\n",
    "\n",
    "### Fitting\n",
    "We then use `iminuit` to optimize $A_{FB}$ and $F_L$. The function `m.migrad()` does the gradient descent minimisation of the negative log likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bin 0: 0.293 pm 0.037, -0.099 pm 0.021. Function minimum considered valid: True\n",
      "Bin 1: 0.652 pm 0.038, -0.187 pm 0.023. Function minimum considered valid: True\n",
      "Bin 2: 0.738 pm 0.041, -0.097 pm 0.023. Function minimum considered valid: True\n",
      "Bin 3: 0.722 pm 0.033, 0.072 pm 0.019. Function minimum considered valid: True\n",
      "Bin 4: 0.651 pm 0.028, 0.215 pm 0.017. Function minimum considered valid: True\n",
      "Bin 5: 0.299 pm 0.034, 0.383 pm 0.02. Function minimum considered valid: True\n",
      "Bin 6: 0.239 pm 0.042, 0.332 pm 0.024. Function minimum considered valid: True\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ht616/Documents/kstarmumu_course/ic-teach-kstmumu-public/venv/lib64/python3.8/site-packages/pandas/core/arraylike.py:364: RuntimeWarning: invalid value encountered in log\n",
      "  result = getattr(ufunc, method)(*inputs, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "bin_number_to_check = 0  # bin that we want to check in more details in the next cell\n",
    "bin_results_to_check = None\n",
    "\n",
    "log_likelihood.errordef = Minuit.LIKELIHOOD\n",
    "decimal_places = 3\n",
    "starting_point = [-0.1,0.0]\n",
    "fls, fl_errs = [], []\n",
    "afbs, afb_errs = [], []\n",
    "for i in range(len(bins)):\n",
    "    m = Minuit(log_likelihood, fl=starting_point[0], afb=starting_point[1], _bin=i)\n",
    "    m.fixed['_bin'] = True  # fixing the bin number as we don't want to optimize it\n",
    "    m.limits=((-1.0, 1.0), (-1.0, 1.0), None)\n",
    "    m.migrad()\n",
    "    m.hesse()\n",
    "    if i == bin_number_to_check:\n",
    "        bin_results_to_check = m\n",
    "    fls.append(m.values[0])\n",
    "    afbs.append(m.values[1])\n",
    "    fl_errs.append(m.errors[0])\n",
    "    afb_errs.append(m.errors[1])\n",
    "    print(f\"Bin {i}: {np.round(fls[i], decimal_places)} pm {np.round(fl_errs[i], decimal_places)},\", f\"{np.round(afbs[i], decimal_places)} pm {np.round(afb_errs[i], decimal_places)}. Function minimum considered valid: {m.fmin.is_valid}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "It is a good idea to check the output of the minimisation. We notice again the invalid value warning so at some point in the minimisation iminuit explored an unphysical set of fit parameters. The fits all show valid function minima though so ultimately iminuit ended up somewhere reasonable. It is important to check that the fitted uncertainties are unreasonable. Very large or very small values would imply that the fit has not converged properly.\n",
    "\n",
    "Let us also look at the profiles for $A_{FB}$ and $F_L$ for a bin (e.g. bin number 0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAADJCAYAAAAuJEhjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABJMUlEQVR4nO3dd3gVZdrH8e+dnpAQWqhp9N4jTRAVRCyoqKio2HCxrMu+q2vZdde1r6uru3YXC6go9oao2EAQ6UiHUAMJLdSEVFKe94+Z6CGm55zMKffnus6Vc2bmzPxOu/PMM02MMSillFJK+ZMgpwMopZRSSrmbNnCUUkop5Xe0gaOUUkopv6MNHKWUUkr5HW3gKKWUUsrvaANHKaWUUn5HGzh1ICIPi8ghEdkvIskiYkQkxOlcSin/JiJdRWS1iBwXkakiMkNEHnY6l1LeSBs4tSQiicAdQA9jTGuHs1wpIrtEJFdEPhGRZlVM209EVopInv23n8u4JiLyuohk2rf7K3juQhHJEpEMEfl7ufGXicgmu+huFJGLvOR13Ski6+1cO0XkznLPfUhE1olIcQWvuY2IfCYie+0GbHJ9XlMV+ZNFZJ6df7OIjK5i2nAReU1Esu3G9e0u44aIyDcickREDorI+yLSxmX8GfZyskQkrYJ5V/peKK9yFzDPGBNjjHnGUwsRkT/Z37Fs+zsXXsl01X3vqqst8+znZYvIGhG5sJ65R9m/ozx73klVTFvpb09ErhCRVPv3kmm/hsblnn+FXfdyRWS7iIxwGefWmlhJ/krrQSXTV/qZVvc5iEiciLxtvx9HReQtl3HNRORdETks1or/W+XfK8cYY/RWixswHMhweZwMGCCkgXP0BI4DpwHRwNvAO5VMGwbsAv4EhANT7cdh9vjpwPtAlP16tgPXuzx/I/AIEAx0BPYBF9jj2gEngHMAAc4D8oCWXvC67gIGACFAV3vcFS7Pv9bO/Slwf7l5twJuBYban29yLV7D/eXnV8W0i4GngEjgEuAYEFfJtP8EFgJNge7AfmCsPe4cYALQ2P4cXwO+cnnuIGASMAVIq2Delb4XevOeG/AtcKPL4xnAw25extnAAfu32BSYDzxWybTVfe+qqy19sGsnMNj+7bepY+4WQJadJwJ4AlhSxfSV/vaABKCFfT8aeAt4xuW5Z9n1ZAhWR0E7oJ09rs41ETgdmF/D11tpPajtZ1rd52Av5ykgFggF+ruMewH42v4OxNrf0aec/q0YY7SBU8mX4R77h3gc65/7eHv4aCAfKAVy7OKSjPUPcAqwF+uf/58bIOOjwNsujzvaP6qYCqYdA+wBxGXYbn7953gIOMVl3F+BhS6P87B6rMoevw/8xb4/GMgst7yDwFCnX1cF0z8DPFvB8JlU8k8dq3HkkQYO0AUodH1tdiG5uZLp9wJjXB4/ROWNvwHA8QqGj6aCBk5N3gu9OXsDvgdKgAK7/nTBMw2ct4FHXR6PAvbX8Lknfe+qqy3lnjvIfm2D6ph7CvCTy+NGdr3uVsG0Nf7tYTVw3gC+cBn2EzC5khx1ronUroFTm3pQ48+0/Odg19k0ILiS6b8EbnV5/Htgrju/k3W96Saqim0HRmC1Rh8AZopIG2PMt1it8r3GmGhjzHUuzzkD6Iz1Zbi7sk0NInKPiByr7FaLjD2BNWUPjDHbsRoCXSqZdq2xv322tfbwX6KVu9/L5fF/gWtEJFREumL1anxrj1sBbBKRC0Qk2O6KLbTnXxfufl3WCxIRrM90Qx1zeUJPYIcx5rjLsDVUnL8p0AaX96ayaW2n4V2vVdWTMeZMrH/Ct9n1Z0tV04vI8KpqjYgMr+SpJ/0G7futRKR5DWJW9L2rqrYgIp+LSAGwFKtnYUUNllOR8rUjF6uWV/Qbqfa3Z79/WVgrupdg1UFEJBhIAeJEZJtYm+2fE5FI+6nurom/UYd6UO1nWsXnMARIBV63N0MtF5GRLvN6HjhfRJrauS7BavQ4Ths4FTDGvG+M2WuMKTXGvAtsxWrVVuUBY0yuMWYdVrfsxErm/Zgxpkllt1rEjMbqjnWVBcTUYdqvgHtEJEZEOgE3YHUpl/kcuBRrbWgz8KoxZrn9ekqw1m7exvoRvw3cZBeXunDn63J1P9b3fXodc3lCbV9r2fgqpxWRPsB9wJ3lx6nAYYz5sapaY4z5sZKnlv9elt2v6Hv5i0q+d9XVFowx59vzPhf42hhTWuMXWXXusux1qh32+xcLxGNt7kqzR7XC2kxzKdZKUz+gP/A3+3nurokVqXE9cJm+ys+0is8hHmvFfR7QGngS+FREWtjjV2HtLnDYvpVgbbZynDZwKiAi14h1pEJZr0ovrO27VUl3ub8LaOvGPCNEJMe+la0d5WBt83TVGGtto7zqpp2K1XjZirUPxiwgw152M6wi9SDWdu0E4GwRudUePxp4HKtrNQwYCbwiLjv7uryORJfXkVPJy3Xn6ypb7m3ANcB5xpjCSpZbL/baT9n35R6sol62pvx5JU+r7WstG1/ptPY/kS+BPxpjFtbyZSgFv/1elt2v6HsJVPm9q7S2uDLGFBljvgTGiMgFlSwjx+WWWIPcZdnrXDvsbHuwauA79qB8+++zxph9xphDWPunnGvnrHFNtKe/x6V2fA6c1PNW0XOoYT0oN321n2kln0M+1mbtV+3x72D9vzvVHv8esAWrcdQYq9dsZiU5GpQ2cMoRa6/7l4HbgOZ2r8p6Tu5mrUiCy/1ErO2jFc3/r+V+qDnV/dM3xiy0u6SjjTFlXZAbgL4u8+2AtaNtRd3WG4A+9maaMn3s4RhjjhhjrjLGtLbnHwQss6frAJQYY94wxhQbYzKwfujn2uP7AQuMMSvsHq/lWF2cv9lEZ4zZ7fI6osuPd/frsp9/A1aDY5Sd3SOMMee79MI9hrUDX9ma8vmVPG0D0EFEXNe6+lLBpiVjzFGs/bv6Vjat/d39FnjIGPNmvV6Q8nnlVowquo2o5Kkn/Qbt+weMMYcrWU6l37tqaktFQrD2u/sN19phjNldXW4RaWTPq6JNtTX+7ZXPZf8WM7D2zfslnsv9ftSwJtrze8yldpwPnNTzVslzqq0H5dTqM+Xkz2FtudcHv329/7O3YOQAL/Hr/wdn1XanHX+/AT2wdrDqinXU0PVAMfaRC1it8oqOonoLq+u1J5CJy85fHsrZE8jG6iJthNViru5ooz9iNRZu4+SjjToCze3Xew7WjoE97XGNsY4uuBKrOLXGOvrgUXv8SHv6fvbj/ljdlHV6/W5+XVdhHVnQvZLnh2L1Sr0NPGzfD3YZH2FnMPb3IaKGr+F+an4U1RLg3/ayxlP1UVSPAT9gHQXRDavAle0o3g5rzanCHdztzy7C/nx32ffDavpe6M07blj7Rnj6KKqx9u+mB9AEa+fmyo6iqu57V1Vt6WYPi7S/f1dj7W83oI6547A2vVxif3//RdVHUVX627NrR6J9P8n+3X3k8twHgeVAS/v3uBCrgQf1qInUbifjSutBbT7T6j4HoBlwFOtIy2CsTXNH+PUos3nAs/bzI7E2T/1Uk9fg8d+L0wG88YZ1SPQR+0v6lP0lqq6BU3YU1X7grgbKeSXWUUO5WN2/zVzGfQn81eVxf2AlVnfjKk4+zO8yO3sesBo4u9xyzrR/zFn263sZiHIZfxuwDau7cwdwh5e8rp1AEVb3bNntJZfxM+zPzvV2ncv48uNMDfPfT80bOMlY/7TysXbkG+0y7ipgg8vjcKzDcLOxDvm83WXcP+yMrq81x2X86RW8nvk1fS/05h03GqCBY8/3dvs7lo2131q4y7gNwFX2/eq+d5XWFqxDm5fadeOYXWPG1zP3aKz9BPPt9yrZZdxL5X7/Vf32HsHqpcm1/07D6tEvGx+K9Y/8GFZNfAaXFSDqWBOpXQOnqnqQaH8WidV9pjX5HLBWONfZ81wBjHAZ1x6YjdWIO4K1Oa+z078VY4x1eK1SSimllD/RfXCUUkop5Xe0gaOUUkopv6MNHKWUUkr5HW3gKKWUUsrvhDgdwFPGjh1rDh065HQMn5OXl+d0BDIzMwFo2bKlw0kgKiqq+olUhVauXDnXGDPW6RyepHWm7rTWnExrTd1VVmv8toEDsGJFXS9pErhWrlzpdASmTJkCwLRp0xxOAgMHDnQ6gs86+fyL/kvrTN1orTmZ1pq6q6zW+O0mKl2rUspx1V3exOdpnVHKK1RYazzWwBGR10QkU0TWuwybICIbRKRURFJchl9lX/up7FZadt0OERkoIuvsq7Y+I4GyWqiUUkqpOvNkD84MrNNDu1oPXAwscB1ojHnLGNPPGNMPmATsNMastke/CPwO6Gzf/HqbvlJKKaXqz2P74BhjFohIcrlhm6DabfMTsa/aKiJtgMbGmCX24zeAi7BO1+8WxSWlhAT77ZY6n9SqVSunIyilAoDWGu9UWmowQHBQ/TbYeONOxpcDF9r322FdB6RMhj2sQiIyBeuaUCQmJla5EGMM105fTnzTSB4d37tegZV7PfTQQ05HUD5ERF7DugpzpjGmlz1sAtY1wboDg4wxK+zhVwF3ujy9D9ZFBVeLyECsnudI4Avgj0avZePXtNZ4p+83Z3L/7A28ccMgOsRF13k+XtV1ISKDgTxjzPpqJ66AMWaaMSbFGJMSFxdX3bJoFRPOJz/vIbugqC6LUyrgHPfO38oMdHO4Un7jraW7OFFcSkKz+h0671UNHOAKYJbL4z1AvMvjeHuYW0wamkTeiRI+XuW2WSo3ePLJJ3nyySedjqHKKSwu4Yx//8DT3251OspJjDELsK5i7DpskzEmtZqnVrg53O61KdscrvyY1hrvk34kj/lbDnLFoERC67n7iNc0cEQkCLgMu+AAGGP2AdkiMsQ+euoa4FN3LbNPfBP6xsfy5pJdaE+090hNTSU1tbr/TaqhfbFuH4dyChmQ1MTpKO5yOb+uUNV4c7iITBGRFSKy4uDBgx6OqDxJa433eXvZbgSYOCih3vPy5GHis4DFQFcRyRCRySIyXkQygKHAHBGZ6/KU04B0Y8yOcrO6FXgF2AZsx407GANcPSSJbZk5LN15pPqJlQpgM37aRYe4Rpza0fdPb1OfzeG12RSulKq5E8WlvLc8nVHdW9EmNrLe8/PkUVQTKxn1cSXTzweGVDB8BdDLfclONq5vWx6es4k3l+xiSIfmnlqMUj5tdfox1qQf4/5xPQiq55ENXqJBN4crpar31Yb9HM49wdVDktwyP6/ZROWUiNBgJgyMZ+76/WRmFzgdRymv9MbiNBqFBXPJwPjqJ/ZyTmwOV0pVb+aSXSQ2i2JEJ/f0Egd8AwfgqiFJFJca3lme7nQUBSQlJZGU5J4WvKq/QzmFfL5mH5cMjCcmItTpOL/h7ZvDN+3L5tIXf2LX4Vx3zE65kdYa77H1wHGW7TzClYMT3dZL7I3nwWlw7Vs0YkTnFsxatptbT++oJ/5z2L333ut0BOXi3eXpnCgp5Zqh3vmPwNs3hzdrFMaajGNMX5TG/Rf0dPfsVT1orfEeby3dTVhwEBPc2Eus/8ltk4YksS+rgO82ZzodRSmvUVxSyswluxjeqQWdWsY4HccntWocwbg+bXlvRTpZ+V55HiGlHJV3opgPV2VwTu/WNI8Od9t8tYFjO7NbS9rERjBzyS6nowS8Rx55hEceecTpGAr4ZuMB9mUVeG3vja+4YXh78k6U8O7y3U5HUS601niH2Wv2cryg2G07F5fRBo4tJDiIKwclsnDrIXYczHE6TkDbtWsXu3ZpQ9MbvL44jXZNIhnVXa/ZUx+92sUypEMzZixKo7ik1Ok4yqa1xjvMXLKbLq2iSUlq6tb5agPHxeWDEggJEt5aqmtZSqXuP86SHUe4ekhSvS96p+DG4R3Ym1XAl+v3Ox1FKa+xNuMY6/ZkcfWQpOouxF1r2sBx0TImgrG9WvPBygzyT5Q4HUcpR72xOI3wkCCuOKX+ZxRV1mbw9i0a8cqPO/XM6UrZZi7ZRWRoMBf1r/Q62nWmDZxyJg1JIiu/iNlr9zodRSnHZOUX8dGqPVzQty1NG4U5HccvBAUJN5yazJr0Y6zafdTpOEo5LiuviM/W7OWi/m1p7IFTUGgDp5xB7ZvRpVW07mzsoK5du9K1a1enYwS0D1ZmkF9UwrXDkp2O4lcuGRhPbGQoryzc6XQUhdYap330cwYFRaVcNdgzBzHoeXDKERGuHpLEfZ9uYE36MfomNHE6UsC54447nI4Q0EpLDW8uTmNgUlN6tYt1Oo5fiQoL4crBifzvh+2kH8kjoVmU05ECmtYa5xhjeGvpbvomNPFYndEenAqM79+OqLBg3tReHBWAfth6kLTDeXpouIdcOzSZIBGmL0pzOopSjlmy4wjbMnO4enCix5ahDZwKxESEcsmAeD5bvZeDxwudjhNw/v73v/P3v//d6RgB642f0mgRHc45vdo4HcUvtY6N4Pw+bXhvRTrZBXriPydprXHOjJ92EhsZyvl92npsGdrAqcT1pyZTVFqqvTgOOHDgAAcOHHA6RkDacTCH+VsOcuXgRMJCtDx4yuThHcgpLOY9vf6do7TWOGP34Ty+3niAqwYnEhkW7LHlaAWrRIe4aEZ1a8XMJbsoKNJDxlVgeHnhTkKDg5jk5jOKqpP1jo9lUPtmTNcT/6kANP2nnQSLcM3QZI8uRxs4VbhxRHuO5J7g45/3OB1FKY87lFPIh6syuGRAO+Ji3Hc9GFWxycPbs+dYPl9v1B4EFTiyC4p4b3k65/dpQ+vYCI8uSxs4VRjcvhm92jXm1R93UlqqJ+ZS/u2Nxbs4UVzKjSM6OB0lIIzu3oqk5lG8snCH01GUajDvLU8n90QJk4d7vs5oA6cKIsKNwzuwLTOHH7YcdDpOwOjTpw99+vRxOkZAyT9RwpuL0xjdvRUd46KdjhMQgoOE64cls2q3nvjPKVprGlZxSSnTF6UxKLkZveM9fwoKbeBU49zebWjdOIJXftS1rIZy2223cdtttzkdI6B8sDKdo3lF3DRSe28a0oSUBGIjQ3lx/nanowQkrTUN6+uNB9hzLJ8bhrdvkOV5rIEjIq+JSKaIrHcZNkFENohIqYiklJu+j4gstsevE5EIe/jlIrLWHv4vT+WtTFhIENcOS2bRtsNs3Jvd0ItXyuNKSg2v/LiTfglN3H41X1W1RuEhXH9qMt9sPEDq/uNOx1HKo15ZuIPEZlGc1aNVgyzPkz04M4Cx5YatBy4GFrgOFJEQYCZwszGmJ3A6UCQizYEngFH28NYiMsqDmSt05aBEosKCefVHPb16Q7jzzju58847nY4RML7esJ9dh/O46bQObr+ar6redcOSiQoL5sX525yOEnC01jScn3cfZdXuY1x/ajLBQQ1TZzzWwDHGLACOlBu2yRiTWsHkY4C1xpg19nSHjTElQAdgqzGmbAeYb4FLPJW5MrFRoVyWksBna/aQmV3Q0IsPOFlZWWRlZTkdIyAYY/jfgh0kNY9iTM/WTscJSE2iwrh6SBKfrdnL7sN5TscJKFprGs6rP+4kJjyECSkJDbZMb9kHpwtgRGSuiKwSkbvs4duAriKSbPfyXARU+u6IyBQRWSEiKw4edO9OwdefmkxxqeGNxXriP+U/Vuw6yur0Y9w4vH2DrVW5mz9sDr9xeHtCgoJ4aYHui6P8z55j+Xy5fj9XDEogOrzhLoHpLQ2cEGA4cJX9d7yIjDLGHAVuAd4FFgJpQKVn3TPGTDPGpBhjUuLi4twaMKl5I8b0aMXMpbvIP6En/lP+4X8/7KBpVCiXDmy4tSoPmIGPbw5v2TiCCSnxfLAigwPaS6z8zBs/pWGM4dphyQ26XG9p4GQAC4wxh4wxecAXwAAAY8xsY8xgY8xQIBXY4lTIG0d04FheER+uynAqglJus/1gDt9uOsCkockePV26p/nL5vCbTutIiTF6XhzlV3ILi3l72W7O6dWG+KZRDbpsb2ngzAV6i0iUvYY1EtgIICIt7b9NgVuBV5wKmZLUlL7xsbymJ/7zqEGDBjFo0CCnY/i9VxbuIDwkKNCuGl7vzeGe2hSe2DyKC/q25a2luzmae8Jt81WV01rjeR+szOB4QXGDHRruypOHic8CFmMVjQwRmSwi40UkAxgKzBGRuQD2pqingOXAamCVMWaOPaunRWQjsAh4zBjjWA+OiDB5RAd2HMplXmqmUzH83o033siNN97odAy/dvB4IR+u2sMlA+NpER1Ql2Wo9+ZwT24Kv+X0juSdKGHGT2luna+qmNYazyopNby2aCf9E5sw0IFTUHhsbx9jzMRKRn1cyfQzsbaN13Q+jjinV2vaxkYwbcEORnVvmGP5lXK3NxanUVRSyo0OrFU57JfN4QAiUrY5/DtjzGxgtj18ClXs7+cpXVrFMKZHK2b8lMbvTuvQoDtkKuVu3206wK7Dedx5dldHlu8tm6h8RmhwEJNHdGDpziOsSDtS/RNUrU2dOpWpU6c6HcNv5Z0o5s0luzireys6BN5lGbx+c/itZ3QiK7+It5boEZueprXGc4wxvPjDdto1iWSsQ6eg0AZOHUwclEDzRmE8N09PzOUJBQUFFBTokSSeMnPJLo7lFXHTyI5OR3ELf9sc3i+hCcM7teDlhTspKNIjNj1Ja43nLN5+mJ93H+PmkR0ICXamqaH9n3UQFRbCDcPb88TcVNZlZDXIRcOUcof8EyVMW7CDEZ1bOLJN3BP8cXP4rWd05MqXl/L+ygwmDQmoncCVn3hu3jbiYsIb9MR+5WkPTh1dMzSJxhEhPK+9OMqHvL1sN4dyTjB1VGeno6gqDO3QnP6JTfjfD9spKil1Oo5StbJq91F+2n6YKSM6EBHq3CkotIFTRzERoVx3anu+2rCfLQf0InnK+xUUlfDSD9sZ2qE5pyQ3czqOqoKI8PvTO5FxNJ9PV+91Oo5StfL899toEhXKlYMTHc2hDZx6uN6+SJ724rjXiBEjGDFihNMx/M47y3Zz8Hih9t74iFHdW9KzbWOe+W6r9uJ4iNYa99u4N5vvNmdyw6ntaeTwUYDawKmHpo3CmDQkidlr9pJ2KNfpOH5j0qRJTJo0yekYfqWgqIQXf9jOoPbNGNqxudNxVA2ICHeM6cLuI3l8sFLPnu4JWmvc7/n524gOD+HaoclOR9EGTn1NHtGe0OAgXpyvF8lT3uv9lRkcyC7kj9p741PO6NqS/olNeOa7rXpElfJ62w/m8MW6fUwamkRsVKjTcbSBU18tYyKYOCiRD1dlsOdYvtNx/MKUKVOYMmWK0zH8xoniUl6ct42BSU0Zpr03PkVE+POYruzLKuCdZbudjuN3tNa414vztxMeEsRkLzmBqDZw3GDKaR0QgWk/aC+O8j4frMxgb1YBU0d1RkScjqNqaVjH5gzp0Izn5m0n/4T24ijvlHE0j09+3sMVpyR6zeVftIHjBm2bRHLJgHhmLU8n87ieNEp5j6KSUp6ft42+CU04rXMLp+OoOrD2xenKoZxC3lic5nQcpSr0vx92IAI3jezgdJRfaAPHTW45vSPFJaW8snCn01GU+sXHq/aw51g+fxzVSXtvfNgpyc0Y2SWOl37YzvGCIqfjKHWSzOwC3l2RziUD4mkTG+l0nF9oA8dNkpo34sJ+7Zi5ZBdHc084HUcpiktKeW7eNnq3i+WMri2djqPq6Y4xXTiaV8T0RWlOR1HqJK/8uJPiklJu9rLLv2gDx41uPb0jeSdKmL5Ie3Hq46yzzuKss85yOobP+2T1XnYfydN9b/xEn/gmjOnRipcX7OBYnq5EuYPWmvo7mnuCmUt2Ma5vW5JbNHI6zkm0geNGnVvFcE6v1kxflKa9OPUwYcIEJkyY4HQMn1Zs73vTvU1jRnfX3ht/cfuYLuScKOblhTucjuIXtNbU32uLdpJ3ooRbT+/kdJTf0AaOm/3pLKsAvaRHVNWZXuG3/j76eQ87D+Xqvjd+plvrxpzfpy3TF6VxKKfQ6Tg+T2tN/RzOKeS1H3dybu/WdG0d43Sc39AGjpt1aRXD+P7tmPFTGvuz9IdTF1OnTmXq1KlOx/BZBUUl/OebLfRNaMLZPVs7HUe52f+N7mxdV0xPLlpvWmvq5/l528kvKuH2s7o6HaVC2sDxgD+N7kKpMTzz/Vano6gA9MbiNPZlFXD32K7ae+OHOsZFc/GAeN5csktXopRj9hzLZ+aSXVw6MJ5OLaOdjlMhbeB4QEKzKK4clMi7y9PZqdeoUg0oK7+I5+dtZ2SXOIZ11PPe+Ks/jupMqTE8/d0Wp6OoAPXMt9YK/B9Hd3E4SeU81sARkddEJFNE1rsMmyAiG0SkVERSyk3fR0QW2+PXiUiEPXyi/XitiHwlIj5RtW87szNhwUE89Y0WINVwXvphO9kFRdw9tpvTUZQHJTSLYtKQZN5dns7m/dlOx1EBZvvBHN5fmc5VQxJp18R7zntTnid7cGYAY8sNWw9cDCxwHSgiIcBM4GZjTE/gdKDIHv40cIYxpg+wFrjNg5ndJi4mnBuGJzN7zV427M1yOo4KAPuzCnjtx51c2LctPdo2djqO8rCpozoRExHKI3M2YYxxOo4KIE99vYWI0GB+f4b3HTnlymMNHGPMAuBIuWGbjDGpFUw+BlhrjFljT3fYGFMCiH1rJNbOBI2BvZ7K7G5TTutIbGQo/55b0UtWlRk3bhzjxo1zOobPefq7LZQawx1jvHOHP08J1N7iJlFhTB3VmYVbDzF/y0Gn4/gkrTW1t35PFnPW7ePG4e295ppTlfGWfXC6AEZE5orIKhG5C8AYUwTcAqzDatj0AF51LmbtxEaGcsvpHZmXepBlO49U/wQFaNGpi22ZOby7PJ2rBieR0CzK6TgNbQYB2ls8aUgSyc2jeHTOJopLSp2O43O01tTeE3NTaRIVyo2nec81pypTZQNHRBKrurkxRwgwHLjK/jteREaJSChWA6c/0Bar6PylirxTRGSFiKw4eNA71miuHZpMy5hwHv9qs3Yj19CxY8c4duyY0zF8yr/nphIVFsIfzvTuLuPK1KfWBHJvcVhIEPec052tmTm8szzd6Tg+R2tN7SzdcZgfthzklpEdaRwR6nScaoVUM34OYLB++GUMEAe0BILdlCMDWGCMOQQgIl8AA4BsAGPMdnv4e8A9lc3EGDMNmAaQkpLiFa2JyLBgpo7qzN8+Wc+81EzO7NbK6Uhe76677gJg2rRpDifxDat2H+WrDfu5/awuNPfyLuMqNFSt+aW32J73O8aYx40xRSJS1lucC2wFfl/RDERkCjAFIDHRnet5dXN2z1YMat+M/3yzhQv7tSXGB/7xeAutNTVnjOHxuam0ahzOtcOSnY5TI1X24Bhjehtj+th/ewPjgEVADvB/bswxF+gtIlF2V/FIYCOwB+ghInH2dGcBm9y43AZx+SkJJDWP4om5Wygt9Yp2l/ITxhge+3IzLaLDmTy8vdNx6qwBa029e4uNMdOMMSnGmJS4uLiKJmlQIsLfz+vB4dwTvKAn/1MeMi81k5W7jjJ1VGciQt21vuFZNdoHR0Q6i8gM4EtgJdDDGPNsNc+ZBSwGuopIhohMFpHxIpIBDAXm2GtRGGOOAk8By4HVwCpjzBxjzF7gAWCBiKwF+gGP1v5lOis0OIjbz+rCpn3ZzF7r9b3eyofMt/fv+uOoTjQKr65D1vvVpdbU0i+9xcaYPKCst7gfWL3FxtqW/B4wzI3L9aje8bFc3L8dr/64k/QjeU7HUX6mtNTwxNwtJDWP4rKUBKfj1Fh1++D0shsqHwLfAr2MMa/YO/9WyRgz0RjTxhgTaoyJN8a8aoz52L4fboxpZYw522X6mcaYnsaYXsaYu1yGv2SM6W6v3Y0zxhyux+t1zLg+benWOoanvtlCYXGJ03GUHygpNfzrq80kNY/iikHObyqpj/rUmlry297iP5/dlSCBx/WoTeVms9fuZdO+bG4/qwuhwd5ybFL1qlvlWwOkY20fHwQMcj31uzFGL+JRQ0FBwl/O7c61ry1j+qI0bh7Z0elIysd9uDKDzfuP88zE/j5VdCpR51pjN4xOB1rYPcT/wNrp+Fms/WzmiMhqY8zZxpijIlLWW2yAL4wxc+z5lPUWFwG7gOvc/SI9qW2TSKaM6MAz32/j+lOTGZDY1OlIyg8UFJXw+FepdGsdw7g+bZ2OUyvVNXAmYxUB5QYju8QxuntLnv1uKxf3b0fLxhFOR/JKl156qdMRvF5WfhH/+mozKUlNGdenjdNx3KHOtcYYM7GSUR9XMv1MrEPFyw9/CXipLhm8xU0jOzJreToPf76RD28Zptciq4bWmupNW7CDPcfymfW7IQQF+db3qboGzjtAjDHmpGOu7W7c4x5L5cfuPa8HY/7zA4/PTeXfE/o6HccrjRkzxukIXu8/32zhaN4J3rhwkL/8E9Na4waNwkP485gu3P3hOj5fu49xfX1rjbuhaa2p2t5j+bwwfxvn9m7N0I7NnY5Ta9X1az8DjKhg+HDgP+6P4//at2jEDcPb88HKDFanH3M6jlfav38/+/fvdzqG10rdf5w3l+ziysGJ9Gwb63Qcd9Fa4yaXDkygZ9vGPDxnIzmFxU7H8Wpaa6r2zy83Ywz89dzuTkepk+oaOAONMR+VH2iM+Rg4zTOR/N8fzuxMXEw493+2QQ8br8B9993Hfffd53QMr2SM4R+frScmIoQ7zvKrSzJorXGT4CDh4Yt6kXm8kKe+1ov9VkVrTeWW7TzC7DV7uWlkR+Kb+ubZ0atr4FT1qnx+r0anRIeHcPfYbqxOP8bHP+9xOo7yIXPW7WPJjiP8eUxXmjYKczqOO2mtcaP+iU25clAiM37ayfo9erFfVTslpYYHZm+gTWwEN4/0/ksyVKa6wpEpIoPKDxSRUwDvuBaCj7q4fzv6JjThsa82azeyqpG8E8U8MmcTPds2ZqKPHxZeAa01bnbX2d1o1iiMez9ZT4n2FKtaeG9FOhv2ZvOXc7sTFea759eqroFzJ/CeiNwvIuPs2wNYJ8G60/Px/FdQkHD/uB4cPF7I8/O2OR1H+YDn521jX1YBD1zQk2AfO5qhBrTWuFlsVCh/O68Ha9KPMWvZbqfjKB+RlV/EE3NTGZTczOeP0KzuUg3LsM5JIVjnhLjOvj/YGLPU0+H8Xf/Eplw8oB2vLtxJ2qFcp+MoL5Z2KJeXF+xkfP92pCQ3czqO22mt8YwL+7VlWMfm/OurzRw8Xuh0HOUDnv52K0fzTnDfuB4+f4RmlX1PIpJojNmNdeIs5QF3j+3GV+v38/CcTbxybYrTcbzC1Vdf7XQEr/PQ5xsJDRb+ck43p6N4hNYazxARHrqoF+f8dyGPzNnIf6/o73Qkr6K15mTbMo/zxuI0rjglkV7tfP8Izeo2UX1SdkdEPvRslMDUqnEEt53ZiW83HWDBFt3VAOC0007jtNP0wJky328+wHebM/nj6M7+fHLIT8ruaK1xr45x0dx8ekc+Wb2XRdsOOR3Hq2it+ZUxhgdmbyQyLJg/j+nidBy3qK6B49o/5bu7Unu5ycPbk9Q8igdmb+BEcanTcRyXlpZGWlqa0zG8QmFxCQ/O3kiHuEZcN8x3rxZeA1prPOjW0zuS1DyKv32ynoIivRZeGa01v/puUyYLtx7iT6O70Dw63Ok4blFdA8dUcl+5UXhIMP8Y14PtB3N5Yb7ucPzoo4/y6KM+d9F4j3hh3nbSDudx/7iehIX49dHSWms8KCI0mIcu7MXOQ7n874cdTsfxGlprLLmFxdw/ewOdWkYzaWiS03HcprqK2VdEskXkONDHvp8tIsdFJLshAgaKM7u14oK+bXl+3ja2HNAz0yvYtC+b5+dtY3z/dpzWJa76J/g2rTUedlqXOM7v04bn52/TgxrUSf79dSp7juXz2MW9/eHCvb+o7iiqYGNMY2NMjDEmxL5f9rhxQ4UMFP8Y14Po8BDu+mCtnrciwBWXlHLnB2toEhXKfef3cDqOx2mtaRj3nd+D8OAg7vlorZ5FXQGwctdRZvyUxjVDkvzuCE3/aar5gebR4fxjXE9Wpx/j9Z/SnI6jHDRt4Q7W78nmwQt7+dsZi5WDWjaO4G/nd2fJjiPM0BoT8AqLS7j7w7W0jY3kzrH+d4SmNnC8zIX92nJ61ziemJtK+pE8p+MoB2zLPM5/v93KOb1ac25v3z7RlvI+l6UkcGa3lvzrq81sy8xxOo5y0PPfb2NbZg6PjO9FdLjvnrG4MtrA8TIiwiPjexMk8NeP12FM4HUjT548mcmTJzsdwxElpYa7PlhLVFgwD1zY0+k4yg+JCI9d3JvIsGDueH8NxSWBe+RmINeaTfuyeWH+di4e0I7Tu7Z0Oo5HaAPHC7VrEsk953Rj4dZDfLAyw+k4DW7w4MEMHjzY6RiOmPFTGqt2H+Mf43rQMsZvz3mjHNaycQQPX9SLNenHeHH+dqfjOCZQa01xSSl3f7iW2MhQ/n6e/+7j57EGjoi8JiKZIrLeZdgEEdkgIqUiklJu+j4istgev05EIkQkRkRWu9wOich/PZXZm1w1OIlTkpvy8JxNZB4vcDpOg0pNTSU1NdXpGA1u1+Fcnpi7mTO7teSifu2cjqP83Pl92jKub1ue/m5rwF5xPFBrzWuLdrI2I4sHLuzp1/v4ebIHZwYwttyw9cDFwALXgSISAswEbjbG9AROB4qMMceNMf3KbsAu4CMPZvYaQUHCY5f0Ib+ohPs/2+B0nAb15JNP8uSTTzodo0GVlhru+XAdoUFBPDK+l89fA6Yh6cpU3T10YU+aNQrjjvfWUFgceCcADMRak3Yolye/3sJZPVpxnp/v4+exBo4xZgFwpNywTcaYiprLY4C1xpg19nSHjTEn/dpEpAvQEljoochep2NcNH8c1Zkv1u3nq/X7nY6jPGjW8t0s3nGYv57XnTaxkU7H8TUz0JWpOmkSFca/LulD6oHjPPXNFqfjKA8rLTXc/eFawkKCePgi/1+R8pZ9cLoARkTmisgqEbmrgmmuAN41AbbX7ZTTOtC9TWPu+3Q9R3NPOB1HecCeY/n884vNnNqpOVeckuB0HJ+jK1P1c0a3lkwclMC0BTtYkXak+iconzVr+W6W7jzCved2p5X/XtfuF97SwAkBhgNX2X/Hi8ioctNcAcyqaiYiMkVEVojIioMH/ePClaHBQfx7Qh+O5RVx5wdrA/KoKn9WXFLKH2f9jDGGf47v4/drVF5AV6YqcO95PYhvGskd768ht7DY6TjKA7YfzOGROZsY1rE5lwfIipS3NHAygAXGmEPGmDzgC2BA2UgR6QuEGGNWVjUTY8w0Y0yKMSYlLs5/Tm3fs20sd5/TjW83HeDNJbucjqPc6OnvtrJi11EeGd+bxOZRTscJBPVemfLHFano8BD+fWlfdh/J4+E5G52Oo9yssLiEP7z9M+EhQTx5Wd+AWZHyljP7zAXuEpEo4AQwEviPy/iJVNN74+9uODWZH7ce5OE5mzgluRnd2/jv2et///vfOx2hQSzadojn5m3j0oHxXNRfj5pqIL+sTAGISNnK1Hf242pXpowx04BpACkpKX7TyzO4Q3NuHtmRF+dv55TkZlw8IN7pSB4XKLXmn19sZuO+bF65JiWg9vHz5GHis4DFQFcRyRCRySIyXkQygKHAHBGZC2CMOQo8BSwHVgOrjDFzXGZ3GQHewBERnpjQl9jIUP4w62fyTvhvN3Lfvn3p27ev0zE86uDxQv7v3dV0aNGIB/WEfg1pLtBbRKLsHY5HAq5dFgG9MnXHWV0Y3L4Zf/14HZv3+/81TgOh1nyz8QAzfkrjumHJjO7Ryuk4DcqTR1FNNMa0McaEGmPijTGvGmM+tu+HG2NaGWPOdpl+pjGmpzGmlzHmrnLz6mCM2eyprL6iRXQ4/728H9sP5vDQ5/7bjbxmzRrWrFnjdAyPKS013PH+GrLyi3juygFEhXlLR6pv0pUp9wkJDuLZK/sTExHKrTNXcbygyOlIHuXvtWZfVj53frCGnm0b85dz/e9aU9Xxln1wVA2d2qkFN4/syKxl6cxZu8/pOB7x/PPP8/zzzzsdw2OmLdzBgi0H+fv5Pfx6U2ND0ZUp92oZE8GzE/uz60ge93zo35eL8edaU1Jq+L93VnOiuJRnJ/YnPCTY6UgNThs4Puj2s7rQL6EJ93y0Vi/I6WNW7T7Kv+emck6v1lw9ONHpOEpVaEiH5tx5dlfmrNvH9EVpTsdRdfDc99tYuvMID17Yiw5x0U7HcYQ2cHxQaHAQz07sDwb++M7PAX2xPF+SlV/EH97+mdaxETx2iR4SrrzbTad1YHT3Vjz6xSZW7jrqdBxVC0t3HObp77Ywvn87LhkQuAcwaAPHRyU0i+KRi3uzavcxnv5uq9NxVDWMMdzz4VoOZBfwzMT+xEaGOh1JqSqJCE9e1pe2TSK57e1VHM4pdDqSqoGjuSf4v3dXk9gsiocC4GzFVdEGjg+7oG9bLkuJ57l52/h+8wGn46gqTF+Uxpfr9/Pns7syILGp03GUqpHYyFBeuGoAh+1/miWl/rs/jj8oLTXc+cEaDuUU8uzEAUSHB/YBDIH96v3AAxf0YuO+bKbOWs1Htw6jS6sYpyPV2x133OF0BLeal5rJw3M2claPVkwZ0cHpOErVSq92sTx4QU/u+Wgd//lmC38+u6vTkdzG32rN43NT+XZTJveP60Hv+Fin4zhOe3B8XGRYMC9fk0JUWDCTX1/uF93IXbt2pWtX/yiiqfuP84e3f6Zb68b89/J+BAUFbnex8l2Xn5LA5SkJPDdvG++vSHc6jtv4U615f0U6L/2wnSsHJ3LtsGSn43gFbeD4gTaxkUy7JoXM7EJumbmKE8W+vdPx0qVLWbp0qdMx6u1QTiGTX19OZFgwr16XQqMA7y5WvktEeOiiXgzv1IK/fLSOhVv94xIV/lJrlu08wl8/XsepnZrzwAU9A3q/G1fawPET/RKa8PilfViWdoS/feLb56549dVXefXVV52OUS+FxSXc/OZKDh4vDLjToyv/FBYSxItXD6BTy2humbmKjXt9/0zH/lBrdh/O46Y3V5DQNIoXrhxIaLD+Wy+j74QfubBfO/5wZifeW5HBqz/udDpOwDLG8JcP17Fi11GeuqwffROaOB1JKbeIiQhl+vWnEB0ewg0zlrMvK9/pSAEtu6CIG15fTqmBV687hdgoPTrTlTZw/MyfRnfhnF6tefSLTXpklUNemL+dj37ewx1ndeG8Pm2cjqOUW7WJjWT69aeQU1jM9dOXk+3nl3PwVsUlpdz29s+kHcrlxasG0L5FI6cjeR1t4PiZoCDr3BXd2zRm6qzVbDlw3OlIAeXLdft4Ym4qF/Vry21ndnI6jlIe0b1NY168egDbMnO41Q/2+/NFD8/ZxIItB3nwwl4M69TC6TheSRs4figqLIRXrk0h0j6yKvN4gdORAsKa9GP86b3VDEhsomcqVn5vROc4HrukDz9uO8RfPvLt/f58zZtLdjHjpzRuOLU9V+olXyqlh3X4qTaxkbx8TQpXvryEq19ZyqzfDaF5dLjTsWrkr3/9q9MRam39niyueW0ZcTHhTLsmhYjQwLuwnQo8lw6MZ8/RfP7z7RbaxEb43DlyfLHWfLp6D/d9up4zusZx73ndnY7j1bQHx4/1S2jCK9emsOtwHle/uoxjeSecjlQjycnJJCcnOx2jxjbvz2bSq0uJDg/h7RuH0MJHGpJKucPUUZ244hTrHDlPfZ3qUz05vlZrZq/Zy5/eXc3g9s144aqBBOt5taqkDRw/N6xjC16+JoXtmTlMenUZWfnev0PgggULWLBggdMxamTLgeNc9fJSwkOCeft3g0loFuV0JKUalIjw6PjeXJ6SwDPfb+Pxub7TyPGlWvPlun3837urGZjUlFevPYXIMO0lro42cALAaV3ieGnSADbvz+a66cvIKSx2OlKVZs6cycyZM52OUa1tmTlc+fJSgoOEWVOGkNRcj2JQgSkoSPjnxb25anAiL87fzqNfbPKJRo6v1JqvN+znD7N+pm98LNOvH6QnDa0hbeAEiDO7teLZiQNYm5HF9dOXkXfCuxs53m7noVyufHkJAG//bogeoqkCXlCQ8PBFvbh2aBIvL9zJg59v9IlGjrf7btMBfv/2Knq1i+X1GwYF/AU0a0MbOAFkbK/WPH1FP1buOsqNr6+goKjE6Ug+adfhXCZOW0JJqWHW7wbTqWW005GU8goiwv0X9GTy8PZMX5TGfZ9uoFSvQF5n81MzuWXmKrq3aczrNwwiJkJP5Fcb2hQMMOf3aUtRSSm3v7eG372xgpf1iJ9aST+Sx5UvL6WwuIRZU4bQ2Q+u3q6UO4kIfzuvOyHBwv9+2EFxaSmPXNRbLzRbSz9uPcSUN1fSuVU0b9wwiNhIbdzUlsd6cETkNRHJFJH1LsMmiMgGESkVkZRy0/cRkcX2+HUiEmEPDxORaSKyRUQ2i8glnsocKMb3j+dfF/dh4dZDTHx5CYf84ArkDeHn3UcZ/8IicgqLmXnjYLq1bux0JIXWGm8kItwzthu3ndGJWcvS+fP7aygs1h7jmvp09R5ueH05HVo0YubkwTSJCnM6kk/yZA/ODOA54A2XYeuBi4H/uU4oIiHATGCSMWaNiDQHyg73uRfINMZ0EZEgoJkHMweMy05JICYihP97dzXjX1jE9OtOoVNL7+iNePDBB52O8Bufr93LHe+toVXjCF67LsVr3isFaK3xSiLCHWO6EB4SxJPfbGH3kTxevHogcTHecxoFb6s1xhj+++1Wnv5uK4PbN+OlqwfStJE2burKYw0cY8wCEUkuN2wTUNEZXscAa40xa+zpDruMuwHoZg8vBQ55KHLAOad3G9o0ieTG11cw/oWfeOnqgUQ4HQpo3bq10xF+YYzhhfnbeWJuKilJTfnfpIE+c8LEQKG1xnuJCH8Y1Zn2cY348/truPC5H5l2TQq92sU6HQ3wrlpTUFTCnR+sZfaavVw6MJ5Hx/cmLER3k60Pb3n3ugBGROaKyCoRuQtARJrY4x+yh78vIq0qm4mITBGRFSKy4uDBgw0Q2/f1S2jCJ78fRtvYSK59bRnf7sxzOhJff/01X3/9tdMxKCox/Pn9tb9cW2rmjYO1ceP76l1rtM7U3vl92vLBzcMAuPSln5i9Zq/DiSzeUmuOFZQw8eUlzF6zl7vHduOJS/to48YNvOUdDAGGA1fZf8eLyCh7eDzwkzFmALAY+HdlMzHGTDPGpBhjUuLi4hogtn+IbxrF+7cMZWjH5ry4IpuZ645T6uDhnR988AEffPCBY8sHOF5YyoMLjvDhqgz+NLoL/7m8n+6M7R/qXWu0ztRNr3axfHrbcHq1jeUPs37m7fXO1hnwjlqzO6uIe747zKZ92bx09QBuOb2jXsfOTbylgZMBLDDGHDLG5AFfAAOAw0Ae8JE93fv2cOVmjSNCmX7dKYzpEMnHm3N5askxCgL0CsEZ2cXc8/1hth4p4ukr+vHH0Z214PgPrTUOiosJ563fDebylAQ+3JTL4z8dI78oMOsMwMp9hfz1+yOUlML7Nw1jbK82TkfyK97SwJkL9BaRKHsnwJHARmOdJWo2cLo93ShgozMR/V9IcBBTBjTm2r4xLMko5I5vDrP5kG9cv8odSo3hi6253PntIfJOlHL/yGZc2K+d07GUe2mtcVh4SDCPXdKbyf1iWLmvkLu/O8yWw4FTZwAKiw2vrc7m0R+P0jo6mMdGN6d3vHfsl+RPPLaTsYjMwioWLUQkA/gHcAR4FogD5ojIamPM2caYoyLyFLAcMMAXxpg59qzuBt4Ukf8CB4HrPZVZWTsFXtClER2bhvLcsiz+Pu8IF3RtxBU9owkN9t9ejMzcEl5YkcW6zBP0bx3GrSmxNIvUTVK+QGuN7xERzu3ciITYEJ5dlsW93x/h/C5RXNEzhvAQ/60zAKmHT/Dssiz25ZQwtmMUk/pEE6H723iEJ4+imljJqI8rmX4m1uGb5YfvAk5zYzRVAz3jwnhqTHNmrDnOJ6m5rNpXyNTBsbRv4l8nmzLG8H1aPtNXH8cAtwxszKj2kbpJyodorfFdvVuG89+zW/Dm2uN8tiWPZXsLuTUllp5x/ndo9IkSw7sbcvgsNZdmUUH847Sm9GmlBy14kp7JWFUqMjSIW1JiGdQunBdWZHP3t4e5rEc047s1ItiDZyV9/PHHPTZvV0fzS3hxZTYr9xXSMy6U358SS6tG+pNQqiFFhQZx08BYhiVE8OKKbO6bf4SxHaO4unc0kaGe7dloqFqz7UgRzy7PIiO7mNHtI7m2bwxRHn5tShs4qgYGtongv2PCePnnbGZtyGHFvkJu7N+YTs0805vTpEkTj8y3TEmp4Yfd+by+5jgnig3X94vh3E5RBGmvjVKO6d0ynKfGNGfW+hzmbM1j5b4CbhoYS//Wnuvl8HStyS8q5ePUXD7enEuTiCD+NqKpR1+POpk2cFSNxIQHcfuQJgxul88rq7K5+7vDnNI2nCt6RpPs5s1Ws2fPBmDcuHFunW+JMfyUXsB7G3PYe7yELs2sXpv4xvozUMobRIQEcX2/xgyNj+D5FVk8vPAoPVqEMqFHNL1bhrl907Gnak1+cSlfbcvj09Rcjp8wnJ4UwQ39GtMoTHttGpJWdlUrpyZE0r91OHO25vHZllzu+OYwQ+MjuLxnNAluaii4u+iUGsPSPYW8uyGH9OxiEmNDuGtYEwa1Ddd9bZTyQt1ahPHkWS34dkceH2/O5YEFR+na3Gro9GvlvoaOu2tNQXEpc7fn88nmHLJPGPq3DuPyntF0buZ/+xT5Am3gqFqLCg1iQo9ozukUxewtuXy+NY8lGQWMSIxgQo9o2sZ4x9fKGMOKfYW8sz6HtKxi2sUEc/uQWIbGR+jmKKW8XFiwdaTV6A5RfL8zn4825/DwwqN0bmY1dAa0dn+PTl0VFhvm7sjjk825ZBWW0reV1bDp2lwbNk7yjv9EyidFhwUxsVcM53VuxKepuXyxLZcf0wvo3zqcYQkRDGob7siOdHuPF7MovYAf0wvIyC6mdXQwUwfFMjwxgmAvKYhKqZoJCxbGdopiVPtI5u3K56NNOTz641GSY0MYkRTJ0PhwRw4OKDWGzYeKWJxRwKL0ArIKS+ndMowrekbTrYU2bLyBNnBUvTUOD2JSnxjO7xLF51vy+DE9n5X7CgkNggFtwjk1IYKBbcI9eq6HzFyrUbMovYCdx4oB6N4ilN+f0piRiZEePepLKeV5ocHCmA5RnJkcyfxd+czdlseba4/z5trjdGwawpD4CIbGR9Am2nP/1kqNYcvhIhalF7Ako4AjBaWEBUH/NuGc37kRPfzw8HZfpg0c5TZNI4KZ1CeGq3pHs/VIET/uLmBxRgFL9xQSHiwMbBNOz7gwEmJDSGwcQkx43Ro8pcaQmVvC7qxidmUVs3JfIVuPFAHQuVko1/WNYVh8BM2j9ER9SvmbkCBhdPsoRreP4kBuMUsyClmcUcBb63J4a10OybEhDGoXTnKTUBIah9CqUXCdV3AKiw0Zx4tJzypm+9Eiluwp4Eh+6S8rb8PiIxjYNpxIPVGfVxLj8MXOPCUlJcWsWLHC6Rg+Z+XKlW6dX4ndjbsoPZ/FGYVkF/563ZkmEUEkNA6xbrEhRNubs04UFgAQFh4BgMFwKK+U3VlFpGcXk55dzImSX5fRvkkIwxIiGBYfQWs3r70NHDjQrfMLJCKy0hiT4nQOT9I6U3furjWZuSUs3WOtVKUeLvpleEgQtI2x6kx84xDaRAcTLILB/KbWnCiBPXaDJj27mMzcEsr+Q4YGQb/Wv/ZIu3vzu9aauqus1mgPjvKoYBF6xoXRMy6M3/U3HM4vJT27mN12AdmdVcR3O/MpLKmooV1w0qOmEUEkxoYwpkMUibEhJMZaBUvXnpRSLRsFM65LI8Z1aUReUanVUMkuJiO7hIzsYrYesTYtVazwl3vBAm1iQujQNJTTkyJJiLUaR22i694TpJyhDRzVYESEFlHBtIgKPulkV6XGcDCvhMJiq5Hz5acfAnDOhZf8Mk3TiOA6b9JSSgWWqNAgOjcL+83h2QXFpRzMK6Vsy8VXn5XVmksRIEggrlEwodqQ8QvawFGOCxI56SiIlT9+D8BN11zhVCSllB+KCAkiofGvK0orFlq1ZsokrTX+SFeJlVJKKeV3tIGjlFJKKb+jDRyllFJK+R1t4CillFLK7/jteXBE5CCwq5ZPawEc8kAcX8oQ6Mv3hgxOL99dGZKMMXHuCOOt6lhn6sobvhd14au5wXezB1ruCmuN3zZw6kJEVjh9YjKnMwT68r0hg9PL95YM6mS++pn4am7w3eya26KbqJRSSinld7SBo5RSSim/ow2ck01zOgDOZwj05YPzGZxePnhHBnUyX/1MfDU3+G52zY3ug6OUUkopP6Q9OEoppZTyOwHXwBGRZiLyjYhstf82rWS6r0TkmIh8Xm64iMgjIrJFRDaJyNSGzuAy/hkRyWno5YvIWyKSKiLrReQ1EQlt4OW3F5GlIrJNRN4VkbCKnu+mDNfa02wVkWtdhk8UkXUistbO2aKBlx8mItPs7+FmEbmkoud7avku4z8TkfW1WbY6mYiMtX9P20TkngrG3y4iG+3v2ncikuQyrsrPx9Pqmb1ERFbbt8+8LPfN9u97tYj8KCI9XMb9xX5eqoic3ZC57eXXKbuIJItIvst7/pI35XaZ7hIRMSKS4jKsbu+5MSagbsDjwD32/XuAf1Uy3ShgHPB5ueHXA28AQfbjlg2dwR6XArwJ5DjwHpwLiH2bBdzSwMt/D7jCvv9SbZdf0wxAM2CH/bepfb8p1kVqM4EWLvO6v6GWb497AHjYvh9UlqWhlm+Pvxh4G1hf2/dfb7+8h8HAdqADEAasAXqUm+YMIMq+fwvwbk0+H2/Obj+ude1qwNyNXe5fAHxl3+9hTx8OtLfnE+wj2ZOd+q3WJLc9XQywAFgCpNT3PQ+4HhzgQuB1+/7rwEUVTWSM+Q44XsGoW4AHjTGl9nSZDZ1BRIKBJ4C76rDsei/fGPOFsQHLgPiGWr6ICHAm8EF1z3dDhrOBb4wxR4wxR4FvgLH82rhrZOdpDOxtwOUD3AD8E8AYU2qMqe3Jseq1fBGJBm4HHq7lctXJBgHbjDE7jDEngHewPptfGGPmGWPy7IdL+PX3VtX3oyHUJ7uTapI72+VhI6BsZ9ULgXeMMYXGmJ3ANnt+DaU+2Z1UbW7bQ8C/gAKXYXV+zwOxgdPKGLPPvr8faFXL53cELheRFSLypYh0diDDbcBnLvNo6OUDYG+amgR81YDLbw4cM8YU248zgHa1XH5NM7QD0l0eZwDtjDFFWA3ddVgNmx7Aqw21fBFpYj9+SERWicj7IlLbz7DOyy9bNvAkkFf+SapWqnqPKzIZ+LKOz3W3+mQHiLDr6BIRucgD+SpTo9wi8nsR2Y7V2zm1Ns/1oPpkB2gvIj+LyA8iMsKzUU9SbW4RGQAkGGPm1Pa5lQmpfU7vJyLfAq0rGHWv6wNjjBGR2rZuw4ECY0yKiFwMvAb85oviqQwi0haYAJxezXSefA/KvAAsMMYsdGj5VfLgZxCK1cDpj7VZ4FngL5TrzfDgexCCtSb8kzHmdhG5Hfg3VmPT48sXkX5AR2PMn0QkuRa5VT2IyNVYm6ZHOp2ltirJnmSM2SMiHYDvRWSdMWa7Mwl/yxjzPPC8iFwJ/A1o8H2c6qqS7PuARGPMYREZCHwiIj3L9fg4QkSCgKeA69w5X79s4BhjRlc2TkQOiEgbY8w+EWmDtS9FbWQAH9n3PwamN3CG/kAnYJu1dYQoEdlmjOnUQMsvm8c/gDjgporGe3D5h4EmIhJi9+LEA3s8lGEPJzck44H5QD97/tvteb2HtR9LQy3/MFbPSdn38H2steOGWv5QIEVE0rBqSEsRmW+MOR1VW3uABJfHFX6fRWQ0VsN0pDGm0OW5p5d77nyPpKxYfbJjjNlj/90hIvOxaltDNHBqlNvFO8CLdXyuu9U5u/3eF9r3V9o9PF2AFZ6JepLqcscAvYD59v+11sBnInJBDZ5bqUDcRPUZv7bErwU+reXzP8HacQ6stZEtDZnBGDPHGNPaGJNsjEkG8so3bjy5fAARuRFr+//Esn2RGmr59n4/84BL6/L8WmaYC4wRkaZiHWU0xh62B+ghImUXdzsL2NRQy7ffg9n8+s9tFLCxAZf/ojGmrf39Gw5s0cZNnS0HOot1ZGAYcAXWZ/MLEekP/A+4oNw+f5V9PxtKnbPbmcPt+y2AU6n9d9iTuV13PTgP2Grf/wy4QkTCRaQ90BlrP8SGUufsIhJn77+J3WvWGasHuiFUmdsYk2WMaeHyf20J1ndmBfV5z40De1Q7ecPah+M7rA/9W6CZPTwFeMVluoXAQSAfq9fmbHt4E2AO1v4Xi4G+DZ2h3LzqchRVfd+DYqw1rdX27b4GXn4H+wu+Dav3ItyD78EN9nK2Ade7DL8Zq1GzFqux0byBl5+EdbTBWns+iQ25fJfxyehRVPW6YR2VuMX+Td1rD3sQq8Bjfz4HXH5vn9X08/HW7MAwrBq6xv472ctyPw1ssDPPA3q6PPde+3mpwDle+J5XmB24xGX4KmCcN+UuN+187KOo6vOe65mMlVJKKeV3AnETlVJKKaX8nDZwlFJKKeV3tIGjlFJKKb+jDRyllFJK+R1t4CillFLK72gDR3kdEZkg1pXa54nIdSLynNOZlFL+RUSm2nVmj9YY/6QNHOWNJgO/M8acUe2USilVN7dinaTz3uomVL5JGzjKUSLyiYisFJENIjJFRO7DOjvuqyLyhD1ZgojMF5Gt9iUilFKqzkTkJawThn4JNHU4jvIQv7wWlfIpNxhjjohIJNbpvEcCZwJ/NsasEJHrgEFY1ynJA5aLyBxjncJbKaVqzRhzs4iMxbrszvlO51GeoT04ymlTRWQN1rVHErCuM1LeN8aYw8aYfKwLTA5vyIBKKaV8j/bgKMeIyOnAaGCoMSbPvqJwRAWTlr+eiF5fRCmlVJW0B0c5KRY4ajduugFDKpnuLBFpZm/GughY1FABlVJK+SZt4CgnfQWEiMgm4DGszVQVWQZ8iHXl7A91/xullFLV0auJK6WUUsrvaA+OUkoppfyONnCUUkop5Xe0gaOUUkopv6MNHKWUUkr5HW3gKKWUUsrvaANHKaWUUn5HGzhKKaWU8jvawFFKKaWU3/l/p3SIe7Tqb1YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 5))\n",
    "plt.subplot(221)\n",
    "bin_results_to_check.draw_mnprofile('afb', bound=3)\n",
    "plt.subplot(222)\n",
    "bin_results_to_check.draw_mnprofile('fl', bound=3)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "We can also plot the resulting fitted PDF over the measured $cos(\\theta_l)$ distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEKCAYAAAAVaT4rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2GUlEQVR4nO3dd3wUdf7H8dcnm4QWuhgpKigoPylGiF3PYEU9ARUURQ4UxYZ6igV7VyzYz/PwRFDBqPRiQyDYTwFRUFSKeIBIBwkt2eTz+2Mm3Aopu8nOzm7283w89pHd2SnvTOCzs9/5zndEVTHGGJM8UvwOYIwxJras8BtjTJKxwm+MMUnGCr8xxiQZK/zGGJNkrPAbY0ySSfU7QDj22WcfbdmyZaWW3bZtG3Xq1IluoCiwXJGxXJGxXJGprrnmzp27XlWb7PWGqsb9o3PnzlpZs2bNqvSyXrJckbFckbFckamuuYA5WkpNtaYeY4xJMlb4jTEmyVjhN8aYJJMQJ3dLU1hYyMqVK9m5c2e589WvX59FixbFKFX4EilXzZo1adGiBWlpaT6lMsZEU8IW/pUrV1K3bl1atmyJiJQ539atW6lbt24Mk4UnUXKpKhs2bGDlypW0atXKx2TGmGhJ2KaenTt30rhx43KLvqk6EaFx48YVfrMyxiSOhC38gBX9GLH9bEz1ktCF32+BQICsrKzdj+XLl3PccccBsHz5csaMGVPmsnfddRft2rXjlltuqdS2R44cyaBBg0p9ryRDuO6++246duxIVlYW3bt357fffqtUJmNMYkjYNv54UKtWLebPn/+naZ9//jnwv8J/8cUXl7rsyJEj2bhxI4FAIKxtBYNBUlPD+3OVZAjXLbfcwoMPPgjA448/zgMPPMBLL70U0TqMaTlkWoXzDO4QpL873/KhZ3sdyZTBjvijLCMjA4AhQ4bwySefkJWVxdNPP/2nebp160Z+fj6dO3fmrbfeYvny5Zx88sl07NiRU045hf/+978A9O/fn6uuuoqjjz6aW2+9da9trVixgpycHNq0acP999+/V4a8vDxycnLo2bMnbdu2pU+fPmgpd1yrV6/e7ufbt2+3ph1jqrnqccT/3hD4fUGpb9UqCkKgEr/mfh3gzKHlzrJjxw6ysrIAaNWqFRMmTNj93tChQ3nyySeZOnXqXstNnjyZjIyM3d8WzjnnHPr160e/fv0YMWIE119/PRMnTgSc3kuff/55qd8MvvrqKxYuXEjt2rU58sgjOfvss8nOzv7TPN988w3ff/89zZo14/jjj+ezzz7jhBNO2Gtdd955J6+99hp169Zl9uzZ5f7expjEVj0Kv09Ka+qpjC+++ILx48cD0Ldv3z8d3ffq1avM5qDTTjuNxo0bA3Deeefx6aef7lX4jzrqKFq0aAGw+zxEaYX/4Ycf5uGHH+a+++7jhRde+NM3CFN2M0Zo00Uoa8Yw8czTwi8iNwKXAwosAC4FmgK5QGNgLtBXVQuqtKFyjsx3xGl/+XCVNzLfnk0ypTXR1KhRY/fzQCBAMBgsd3sXXHABF1xwgRV+Y6oxzwq/iDQHrgcOU9UdIvI20Bs4C3haVXNF5CVgAPBPr3L4pW7dumzdujWseY877jhyc3Pp27cvo0eP5sQTTwxruenTp7Nx40Zq1arFxIkTGTFiRKWyLl68mDZt2gAwbdo02rZtW6n1mPgVzonXPdm3lurL65O7qUAtEUkFagOrgZOBse77o4AeHmfwRceOHQkEAhx++OF7ndzd0/PPP8+rr75Kx44def3113n22WfD2sZRRx3F+eefT8eOHTn//PP3auYJ15AhQ2jfvj0dO3Zk5syZYW/fGJOYPDviV9VVIvIk8F9gB/AhTtPOZlUtaW9YCTT3KoPX8vPzy5yWlpbGzJkzy1x29erVu58feOCBpc47cuTIMpfv378//fv3LzdXTk4OOTk5u6e/8MILpc4/bty43c/jdSgJY0z0SGnd+6KyYpGGwDjgQmAz8A7Okf59qtranWd/4D1VbV/K8gOBgQCZmZmdc3Nz//R+/fr1ad26dYU5ioqKwu4rH0uJlmvJkiVs2bLFh0SO/Pz83d1U/bBgVem/e2YtWLNj7+kdmtf3OFH59txfZeUvT6S/QzjbCN1ffu+jUH7/+ypLVXN16dJlrqru1RTg5cndU4FfVHUdgIiMB44HGohIqnvU3wJYVdrCqjocGA6QnZ2toUeuAIsWLQrryDRej2ATLVfNmjU54ogjfEjkKLkmwS+l9dwBp1fPsAV7/zda3ifH40Tl23N/lZW/PJH+DuFsI3R/+b2PQvn976ssXuXyso3/v8AxIlJbnO4mpwA/ALOAnu48/YBJHmYwxhizB88Kv6r+B6dpZx5OV84UnCP424CbRGQJTpfOV6qwjSgkNRWx/WxM9eJpP35VvRe4d4/Jy4CjqrrumjVrsmHDBhua2WMl4/HXrFnT7yjVQ8mHqP2bNT5K2Ct3W7RowcqVK1m3bl258+3cuTMui1Yi5Sq5A5f5s3QKaVDwO51kB/VlGw3Id35KPrw7G3Zsgp2bnZ+7H5shJQC1GkHtRs7PWg3+97x2I6jV8M/v124MdfaxDwsTNQlb+NPS0sK6I1ReXp6vJyXLYrkSi1BMS1nD4bKUw1OWkpWylMNkOTV+DdK3RikLfFffKeA1Gzg/6+/vFPiaDUCLYPvG/30YbFgKK792phUXlh4gIxOaZ0OLztDiSGh2BNSIv84BJjEkbOE3xlP5a2HV3N2P+TW+pL5sB2Cb1mCBHsSrRV05qHlzRq9owGatw2Yy2KJ1+IM6LLvvnMi3qQoF22DHRveDwf1wyF8Lv813Phx+cnvOSAo0aQstst0PhGzndUr8dRE28ccKvzEAm36FRVOc4rpqHmxxhsZGApB5GNOKjmG+Hsz84tYs0eYUu/0iBtcLMrs4Sv+NRKBGhvNocEDp82zf6ORb+TWsmgM/TIZ5rznvpWc43wRaZNNoSwYEj4PU9OhkM9WKFX6TvLZtgB8mwHfvwIovnWkNDnCOno++Epp3hqYdIb0Od1SiH7wnajeCNqc6D3C+JWxY6nwIrJzjfCB8/jwdi4Ow5AVody50vBBnnEQ7R2AcVvhN3CltQLGyhj8uEfaAYgXb4Kf34Lu3YekMKA5Ck/+DU+6B9j2h4YGVje0PEdintfM4vLczrXAH303+Bx11Ecx/E+aM4OP0JkwsPp5JRcezVBN2lBQTJVb4TfVXVAjL8pxi/+M0KNwG9ZrDsddChwsgs1316jGTVouNjbMh52bYtRUWTWX5uH9wbWAS16dO5LviVkwsOoEpRceyjgZ+pzU+sMJvqidVp9nju7fh+wmwfb3To6ZjL+jQCw44DlKS4M6jNepC1kX8LbceTdhEt8AX9Ah8yj1pr3Nn6ht8WtyBCUUn8GFxNtuJv+7FoSIdWtqGlS6bFX5TvRQFYeE4+ORJWP8zpNaEQ7pCxwug9amQWlrfy+Swjoa8UnQWrxSdRWtZSY/AZ/QIfMYzgRfZrjWYWHQcLxb1YKU28Tuq8ZgVflM9FAVhwTvw8ROwcSlktofuL8L/nQM161W8fJJZoi14Mnghw4K9yJafOS/wCecHPqFX4GPeKfoLLxb18Dui8ZAVfuO5ytz9KVwBiuiR8hn84y7YuAwyO8CFb8ChZydHU04VKSl8rW35OtiWZ4PncVXqFC4KzKRX4GOYPA/+cnPZXUtNwrLCbxJSgCLODXzKoMBEWqasgfQOcOFoOPQsK/iV9DuNuS/Yn5eC53B16mT6ffsmzB8DR/SBEwfbB0A1YoXfJJRUgrsL/oEpa1lY3JIrCm7i5SvviahnjpffQhLd7zTm3uCl9LvlGfj0aecCsW9GwxGXuB8A+/sd0VSRFX6TEFI0yAWBTxgUmMgBKetYUNySAQWDmVHcCZDq1R0zXtRvAWcPgxNuhE+ecj8A3oBOfeGEm+wDIIFZ4TdxTjkn5QsuWZ5L/bT1fFt8EPcV9GNm8RHYlagxUr8F/PUp5wOg5BvAvNedD4Cc2yFjX78TmghZ4TdxqxnreShtBCcH5rM2cCCXbr+UWcVZWMGvuKmqoiudK6XB/iEfAE85xX/heDjjEci6OLrbMp7y7CyYiBwqIvNDHn+IyN9FpJGITBeRxe7Phl5lMIkphWL6BT5geo1bOCZlEQ8U9uXt/e9nlh3lx4cG+8Nfn4arP3NGBJ10Dbx+Li1krd/JTJg8O+JX1Z+ALAARCeDcVH0CMASYoapDRWSI+/o2r3KYxHKIrOCxtJc5ImUJs4s6cmdwACu1CYMlWO5ydrI2+sLZp8Ig+gT+jyFL3+TD9M8ZFuzFq0Vdd49eauJTrP46pwBLVfVXoDswyp0+CugRowwmjqVTyI2p7zA1/Q4OlN+5oeAa+hXeZleRxjklhTeKTuO0XU/wRfFh3J32BuPT7+VQ+a/f0Uw5YlX4ewNvus8zVXW1+/x3IDNGGUycypYfeTf9dm5IncDU4mM5ddeTTCo+AWvWSRyracyAwpu5vmAQ+8tapqbfyU2pb5NOGXcUM74SLbn5s1cbEEkHfgPaqeoaEdmsqg1C3t+kqnu184vIQGAgQGZmZufc3NxKbT8/P5+MjIxKLeulZMq1YNWWUqenF23n2A1v03HLTP5I3YdZ+/bnv3U6ljpvZi1YsyOqsaKirFwdmtf3dLtl7dMSfu6vmkVbOXHdGNpu/YyNaU2ZmTmA1bUO2StXpPuoot95T5Gsv7r+f+zSpctcVc3ec3osCn934FpVPd19/ROQo6qrRaQpkKeqh5a3juzsbJ0zZ06ltp+Xl0dOTk6llvVSMuUqra34tJQ5PJA2kkw28WpRV4YFe5U7OuTgDkGGLYi/Tmhl5fJ6ZMhwevX4vb9OSvmWh9NeoRkbeL3oVB4P9uaqDmm7c0W6j7wcnbO6/n8UkVILfyyaei7if808AJOBfu7zfsCkGGQwcaIBW/lH2jO8nP4Um7UO5xXcz4PBvnE/JLCJ3Oziwzl91+OMKjqdvoGP+LDGrRy47Vu/Yxk8LvwiUgc4DRgfMnkocJqILAZOdV+bJNBOfmFK+l2cmjKPJwovoFvBw8zX1n7HMh7aTk3uD/ajZ8G9bNOadPttGENSxxCgyO9oSc3T74Kqug1ovMe0DTi9fEwS6RXI46HUV1lPPS4ouIdvreAnlXl6CH8teISxmaO4astUOsgvsO04qLOP39GSUvw1mprqJbiLR1L/zcWpM/m0qB3XF17HRmx8/D0lw3UIBaQxe9/+vLa+DQ+njoB/nQQXvubc1N7ElF1lYbyzZSWM6MrFqTN5MdiNfoVDrOgbxhadxPkF94KkwIiuMHdUxQuZqLIjfuONZXkw9jIIFnBlwY18UHyk34lMHFmoB5G15k6eS3uBv0y5njETJnBfsD8FpPkdLSnYEb+JLlVnBMfXz4U6TWBgnhV9U6rN1KV/4W38I9iNi1Nn8Xb6/TRlg9+xkoIVfhM9O/+Aty6Bj+6Dw3rA5TNgHzuJa8pWTApPBHtzZcGNHCyrmVrjDo5N+d7vWNWeFX4THWsXwctd4Kf34IxHoecIqBF/V0Ka+PRB8ZF0L3iQjVqPN9Ie4YrAVMDbi0uTWYWFX0TqiEiK+/wQEekmItYQZ/5n4Th4+RTniL/fFDj2GrsjlonYMm1Gj4IHeL/4SO5MG8MLac9Rhzgcp6MaCOfk7sfAie64+R8CXwMXAn28DGZiJ7QrYTg38Nh9KXxxMUy/G754AfY/BnqNhHpNPUxqqrtt1OLawhv4tngqt6Xmcmj6Sq4svJFl2szvaNVKOE09oqrbgfOAF1W1F9DO21gm7hUVwvgrnKJ/5BXOkb4VfRMVwvCic7ik8A4ayR+MS7+PLFnid6hqJazCLyLH4hzhlxwKBryLZOJewXbIvRgWjoVT74Ozn4TUdL9TmWrmi+J29Ch4gD+0NqPTH+a4lIV+R6o2win8fwduByao6vcichAwy9NUJm7VYxu8cR4sng7nPOvcf9UYj6zQTHoW3MsK3ZdX0x7njJSv/Y5ULVRY+FV1tqp2A553Xy9T1es9T2biThM281b6g7ByjtOe37m/35FMElhHQy4suJvvtSUvpj1Dr0Ce35ESXji9eo4VkR+AH93Xh4vIi54nM3GlhazlnfT7OUDWQJ+3oV0PvyOZJLKFDPoU3MFnxe15Im04AwLVf2wjL4XT1PMMcAY4l9Sp6rfAXzzMZOLMIbKCcen30UDyuaTgDjj4ZL8jmSS0g5pcXngz04qO4u600QxOfRvr6185YY3Vo6or5M/9sm0w7STRSX5mRPoT7CSdXgX3slhb+B3JJLEC0riu8Hr+0Fe4LnUiDcjnnmB/1K5FjUg4hX+FiBwHqHvh1g3AIm9jmXhwYsp3/CvtadZoA/oW3s5K3dfvSMZQTAq3By9nCxlclTqFerKdwYVXEdyjnEUy1PXgDkFyopwznoXzMXkVcC3QHFgFZAHXhLNyEWkgImNF5EcRWeSeL2gkItNFZLH7c68brRv/nZ3yJa+kPcFy3Y9eBfdZ0TdxRhgavIihhb3pHvic4WlPUZNdfodKGOEU/kNVtY+qZqrqvqp6CfB/Ya7/WeB9VW0LHI7zTWEIMENV2wAz3NcmjlwUmMHzac8zX1vTu+Au1lPf70jGlOqlom7cXjiAnJRveS19qNPd2FQonML/fJjT/kRE6uOcBH4FQFULVHUz0B0oufPCKKBHOEFNLCidN07h0bRXyCs+nL8VDOEP6vgdyphyvVl0CtcXDiJLlvBm+kPswxa/I8U9US39rLh7te5xOBdwPR3yVj3gXFU9vNwVi2QBw4EfcI725+KcH1ilqg3ceQTYVPJ6j+UHAgMBMjMzO+fm5ob/W4XIz88nIyP+RomMp1wLVjn/UY5Z/w5HbprCT3WP5aPMKyiW+LlPT2YtWBOH43WVlatD88i+JZX8DaIl0fZXNByw7TvOWv0c+akNmdj8NvLTwr+fb2Yt2LdR/H2zrWqd6NKly1xVzd5zenmF/yQgB6eN/6WQt7YCU1R1cXkbFJFs4EvgeFX9j4g8C/wBXBda6EVkk6qW286fnZ2tc+bMKW+WMuXl5ZGTk1OpZb0UT7laDpnG1YHJ3JaWy8J6OZyz9vK46yUxuEOQYQvi54OoRFm5dg9kF6Zo33M30fZXtHSSnxmZ/jjrtR69Cu5lQ5jNlIM7BLmuT3fPclVWVeuEiJRa+Mv8C6jqbGC2iIxU1V8rsc2VwEpV/Y/7eixOe/4aEWmqqqtFpCmwthLrNlHUJ/ARt6XlMqnoOFbs2x9dG19F35hwzdNDuLTgFt5If5TX0odyUcFd1lxZinD+h28XkSdE5F0RmVnyqGghVf0dpyvooe6kU3CafSYD/dxp/YBJlQluomTBWB5MfZUZRUcwuPAqVKzom8Q2Vw/lysIbaSMr+Xf6k9bbpxTh/C8fjTNcQyvgfmA5zpj84bgOGC0i3+F0A30EGAqcJiKLgVPd18YPP70PE67kK23LNYU37NUP2phE9XHx4dxYeC3Z8jMvpT1DGkG/I8WVcP6nN1bVV0TkhpDmn7AKv6rOB/ZqX8I5+jd+Wv4pvNMPMttz+S+D2IUNq2yql2nFx1A3uJ2haf/maV7k+sJBFMfZuSu/hLMXCt2fq0XkbBE5AmjkYSbjtVXzYExvaHAgXDKefGr7ncgYT+QWncwjhRfx18CXPJQ6AhvbxxHOEf9Dbp/8wTj99+sBNgh7olr7I7xxPtRuCH+bCHUa+53IGE8NLzqH+rKNa1Mns4U6PBa8yO9Ivquw8KvqVPfpFqCLt3GMpzYth9d7QCAN+k6EenYfU69Eu3umqZonghdSj+1cnTqFLVqHl4q6+R3JV2UWfhF5nnK+F9nNWBLM1t/htR5QuAMufRcaH+x3ImNiSLgn2J96sp0habn8QR3GFCXvqcby2vjn4FxtWxPoBCx2H1lgZwITyvaN8Pq5kL8WLhkHme38TmRMzCkpDC68iplFWTyUOoK/pnzhdyTflHcB1ygAEbkaOEFVg+7rl4BPYhPPVNmufBjdCzYsgT7vQIvSOlkZkxyCpHJN4Q2MSn+Mp9NeJL+wFnnFWX7HirlwevU0xDmhWyLDnWbiXeFOyL0YfvsGer4KB+X4ncgY3+2kBpcX3MyPuj//THuGI+VHvyPFXDiFfyjwjYiMFJFRwDycC7FMPCsKwrgB8Mts6PEi/N9f/U5kTNzYSm36FQzhN23MK+lPsM/O5X5HiqkKC7+qvgocDUwAxgPHljQDmTj2we3w41Q48wk4vLffaYyJOxupxyUFd7CV2nT/7Qmn11uSKLPwi0hb92cnoBmwwn00c6eZeDVnBHw1HI4dBEcP9DuNMXFrNY3pW3A7KVoEb14Mu7b6HSkmyuvHPxi4AhhWynsKnOxJIlM1v3wC794CrU+D0x7wO40xcW+ZNuP9/QbRY/UwGH8lXPgGpFTvoR3K69VzhfvTLtpKFBt/gbf/Bo0Ogp6vQErA70TGJIQVddrDGY/A+7fBrIfhlLv9juSp8i7gOq+8BVV1fPTjmErb+Qe8eRFoMSetuopf7/vU70TGJJajr4Q1C+GTJyHzMGh/vt+JPFNeU8857s99cW7BWDIGfxfgc5wTvSYeFBfB+IGw/mfoO55fh9sNp42JmAic/ZRzzcvEa5xvzs2O8DuVJ8psyFLVS1X1UiANOExVz1fV84F27jQTL2Y+CD+/B2c+Zn31jamK1HS44HWo08Q52bv1d78TeSKcMxj7q+rqkNdrgAM8ymMi9e1b8OnT0PlSOPJyv9MYk/gymkDvMbBzM+T2cS6ErGbCKfwzROQDEekvIv2BacBH4axcRJaLyAIRmS8ic9xpjURkuogsdn/aVcCVtXIOTL4OWp4IZz3hfFU1xlRd045w7r9g1RyY+nfQ6jWOfzgXcA0C/gUc7j6Gq+p1EWyji6pmhdzpfQgwQ1XbADPc1yZSW1Y5wzHU3Q96jXKGWjbGRM9h3SDnDvj2Tfj8eb/TRFVYN1l1e/BE62RudyDHfT4KyANui9K6k0PBdqfoF2yDv02ym6kY45WTboW1P8D0e6BJWzjkdL8TRYVoBV9h3G6dj+H07hH3oapar9wFnWV/ATbhXPD1L1UdLiKbVbWB+74Am0pe77HsQGAgQGZmZufc3NwIfq3/yc/PJyMjo1LLeqnSuVQ57IcnabLuMxa2v5MN+xy51ywLVm2pdK7MWrBmR6UX94zliozlikx5uVKLd9Fz5YPUK1zHO/vfw6b05gB0aF7f81xVrV9dunSZG9Lasls4hX8JcI6qLop0oyLSXFVXici+wHTgOmByaKEXkU2qWm47f3Z2ts6ZMyfSzQOQl5dHTk5OpZb1UqVzzX4CZj0Ep94PJ/y91FmqcvenwR2CDFsQ1hfBmLJckbFckakoVzPWM6nGXeRrLXoUPMgWMlg+9GzPc1W1folIqYU/nJO7aypT9AFUdZX7cy3OIG9HAWtEpKkbqimwtjLrTkqLpjhFv2NvOP4Gv9MYkzR+Yx+uLLiJZrKBF9KeI0CR35GqJJyP3jki8hYwEdhVMrGiK3dFpA6Qoqpb3eenAw8Ak4F+OMM99wMmVS564mo5ZBqDOwTpH+aR+fKhZ8PvC5yLtJpnwznPWg8eY2Jsnh7CXcHLeCJtOHfpG0Di3rc3nMJfD9iOU7hLKBWf7M0EJjjN+KQCY1T1fRH5GnhbRAYAvwIXRJw62Wzf6FxMUrMB9B4NaTX9TmRMUnqnKIdDZQWXp74Hc0dC5/5+R6qUCgu/e/VuxFR1GU73zz2nbwCS9y7HEVOnr/7W1XDZB073TWOMbx4NXkwbWcVJ794CzTo5ff4TTIVt/CJSU0SuFZEXRWREySMW4Qz0C3zo3FDltPuhRWe/4xiT9IoIcGPhNVC7MYy91LmvdYIJ5+Tu68B+wBnAbKAFkBx3K/BZO/mFO1JHwyFd4Zhr/I5jjHFtpB6c9zJsXObc/yLBhFP4W6vq3cA295aLZ+PcitF4qA47eCHtOecfWPcX7WSuMfGm1Ynwl1vh2zHwbeWuM/JLOIW/0P25WUTaA/VxLuYynlEeTnuFA2Qt1xcMsitzjYlXf7kFDjwept4E6xf7nSZs4fTqGe4OpHYXTlfMDOAeT1MluV6B2fQIfM6wwp58rW2rdEGWMcZDgVSnyeelE5z2/gEfJUSvu3AGafu3qm5S1Y9V9SBV3VdVX4pFuGTUWlbyQOpIPitqxz+KevgdxxhTkfrNocc/nWttpifGLRvD6dXziIg0CHndUEQe8jRVkqpBAS+kPc82avL3wmsoDqslzhjju0O7wjHXwlfDYdFUv9NUKJzKcqaqbi55oaqbgLM8S5TE7k19jbYpK7ip8GrWYbcpMCahnHovNM2CSdfC5hV+pylXOIU/ICI1Sl6ISC2gRjnzm0o4O+VLLk6dyT+D5/Bx8V7XvRlj4l1qDeg5wrkH9rgBUBT0O1GZwin8o3HuwjXAHWZhOs44+iZKDpA1PJr2MnOL2zAs2MvvOMaYymp8MJzzDKz4D+Q94neaMoUzZMNjIvItcKo76UFV/cDbWMkjjSDPpz2PIlxfMIhgePfGMcbEqw49YVkefPKUc1vUg7v4nWgv4d6B633gfY+zJKXbUt/k8JRlXFlwI6to4nccY0w0nPk4rPjKGVH36s8gI74ufbLDSx+dnDKPy1PfY2TwdD4o3vtOWsaY+FXR9TWHyGVMTr+Lrx47j36Ft6GkxOTmLeGw/oI+qVO4kWFpL/F98YE8GrzY7zjGmCj7Wffn/uDf+EtgAVcG4quLZ5mFX0RmuD8fi12c5BCgiDN+/yfpFDKo8Hp2ke53JGOMB94sOpmpRUdzc+rbdJKf/Y6zW3lH/E1F5Digm4gcISKdQh+xClgdXZc6geY7f+LOwgH8ok39jmOM8Yxwe+EVrNbGPJf+AuzY5HcgoPzCfw9wN84wzE8Bw0IeT4a7AREJiMg3IjLVfd1KRP4jIktE5C0RSa7D3d++YVBgIj/WPZ6JxSf4ncYY47Gt1GZQ4XVksiluhnAus/Cr6lhVPRN4XFW77PE4OYJt3ACE3qz9MeBpVW0NbAIGVCp5IgrugonXsJ76fNzkEr/TGGNi5FttzQvBHrDgnbgY0iGcQdoeFJFuIvKk+/hruCsXkRY44/f/230twMnAWHeWUUCPiFMnqtmPw9ofGFJ4ObsCdfxOY4yJoX8UdYf9OsDUG537aPtIVLX8GUQeBY7CuYIX4CLga1W9o8KVi4wFHgXqAjcD/YEv3aN9RGR/4D1VbV/KsgOBgQCZmZmdc3Mrd6OD/Px8MjIyKrVsNNX9YzGd5t3K7/vlMLZufzJrwZodfqfam+WKjOWKTLLnOrbeBjrNu5l1TY5n0WGDK5y/qvWrS5cuc1U1e8/p4fTjPxvIUtViABEZBXwDlFv43W8Ga1V1rojkRBpYVYcDwwGys7M1JyfiVQCQl5dHZZeNmuAu+NcQqLsfTfu9yrD7P2NwhyDDFsTfZRSWKzKWKzLJnmv50Mug7loy8x4ls8uV8H/lN6B4Vb/C7cffIOR5/TCXOR6nR9ByIBeniedZoIGIlOzhFsCqMNeXuPKGwrpF0O05qNXA7zTGGD+dONj3Jp9wCv+jwDciMtI92p8LPFzRQqp6u6q2UNWWQG9gpqr2AWYBPd3Z+gGTKpU8UayaC589A1mXQJvT/E5jjPFbIM25j/aOjfDebb5ECOfk7pvAMcB4YBxwrKq+VYVt3gbcJCJLgMbAK1VYV3wr3AkTr4GM/eCMCj8rjTHJomlH5369C96GH2N/a9VwB2lbjXO/3UpR1Twgz32+DOdkcfU3eyis+xH6jLUmHmPMn504GH6cClP+DgccC7UbxWzTNlaPV1bOhc+ehSOsiccYUwofm3ys8HuhcCdMvBrqNoUz4vdmDMYYn/nU5FNu4XeHW/gxVmGqjbxHYf1PTi+emuF2gjLGJKWSXj5T/h6zXj7lFn5VLQJ+EpEDYpKmOljxNXz+HBzRF1qfWvH8xpjk5kOTTzgndxsC34vIV8C2komq2s2zVAkk9GYMNShgWvod1JKGdP3iJLZ+Efuz9caYBFTS5JP3KLTrAW29vWFLOIX/bk8TVCM3po6ldcpv9C0YwlZq+x3HGJNIYtjLJ5x+/LOB5UCa+/xrYJ5niRLUEbKYKwLTGBPswifFHf2OY4xJNDFs8qmw8IvIFTijaf7LndQcmOhhpoRTgwKeTHuJ1TTmkWAfv+MYYxJVjHr5hNOd81qccXf+AFDVxUB83TLeZzelvsPBKasZUngF+dbEY4ypipBePqmFWz3ZRDiFf5eqFpS8cAdYK38s5ySSJUu4PPAuY4In82lxB7/jGGMSXUiTT5vFL3uyiXAK/2wRuQOoJSKnAe8AUzxJk2iKgjyS9gpracgjwYv9TmOMqS7cJp/MtbNh5Zyorz6cXj1DcG6PuAC4EngX945aSe/rlzks5VeuLrjBmniMMdF14mAWrA/QocVe91GpsgoLv6oWu8Mx/weniecnrei2Xcngj9Uw82FmF3XkveLkGHPOGFM1odf9hGNk1yM9yVFh4ReRs4GXgKWAAK1E5EpVfc+TRIniwzuhqIB7gv1xdosxxiSGcJp6hgFdVHUJgIgcDEwDkrfwL50FC8dBzu38+v5+fqcxxpiIhHNyd2tJ0XctA7zpY5QIgrvg3ZuhYSs4/u9+pzHGmIiVecQvIue5T+eIyLvA2zht/L1wrt4tl4jUBD4GarjbGauq94pIK5x78DbGuY1j39DuonHvs+dgwxK4ZByk1fQ7jTHGRKy8I/5z3EdNYA1wEpADrANqhbHuXcDJqno4kAV0FZFjgMeAp1W1NbAJp8dQYtj4C3zyJBzWw0beNMYkrDKP+FX10qqs2O35k+++THMfCpwMlHR6HwXcB/yzKtuKCVV471ZISYWuj/qdxhhjKk0q6pnpNs1cB7Qk5IMinGGZRSSA05zTGvgH8ATwpXu0j4jsD7ynqu1LWXYgMBAgMzOzc25ubni/0R7y8/PJyMgIe/4Fq7aUOv2g/LmcvfpZPtnnIuY3PLNSWUJl1oI1O6q8mqizXJGxXJGxXJFpVT8QUf3aU5cuXeaq6l4XAoTTq2ci8ArO1brFkWzUvZFLlog0ACYAbSNYdjgwHCA7O1tzcnIi2fRueXl5RLJs/1L62dZiJx/VeINFuj+XrjyT4Mqw7lFfrsEdggxbUPX1RJvliozliozliszIrnUiql/hCuc33amqz1VlI6q6WURmAccCDUQkVVWDQAtgVVXWHQs3pE6guWzg+oJBBMPaZcYYE7/C6c75rIjcKyLHikinkkdFC4lIE/dIHxGpBZwGLAJmAT3d2foBkyoXPTbayEoGBN7lrWAOc/VQv+MYY0yVhXP42gHoi3NStqSpp+QkbXmaAqPcdv4U4G1VnSoiPwC5IvIQ8A1OM1KcUh5KG0E+tRga7O13GGOMiYpwCn8v4KBI+9qr6nfAEaVMXwYkxOA256Z8ytEpPzKk8HI2Uc/vOMYYExXhNPUsBBp4nCPu1COfO9JGM6+4NW8V5fgdxxhjoiacI/4GwI8i8jXORVlAeN05E9nNqe/QiK30KxyChvX5aIwxiSGcwn+v5yniTEdZyiWBjxhZdAY/aEu/4xhjTFSFMx7/7FgEiRcpFPNQ2gjWUZ+ngj0rXsAYYxJMOOPxb+V/99hNxxl6YZuqVsuznRcHZtAx5ReuKxhkd9UyxlRL4Rzx1y15LiICdAeO8TKUb7at59bUt/ikqD1Tio/1O40xxngiorOW6pgInOFNHJ/lDaU2O7kv2A+7q5YxproKp6nnvJCXKUA2sNOzRH5Z9xPMGcGYolNYqs39TmOMMZ4Jp1fPOSHPg8BynOae6mX6PZBeh2d2nO93EmOM8VQ4bfxVGpc/ISydBT+/D6c9wMYp1fKctTHG7FberRfvKWc5VdUHPcgTe8VF8OFd0OAAOOpKmDLD70TGGOOp8o74t5UyrQ7OrRIbA9Wj8M8fA2sWQs9X7R66xpikUN6tF4eVPBeRusANwKU4N0ofVtZyCWVXPsx8EFocBe3O9TuNMcbERLlt/CLSCLgJ6INzf9xOqropFsFi4rNnIX8NXDgaxLpvGmOSQ3lt/E8A5+Hc/rCDquaXNW9C2rIKPn8e2p8P+x/pdxpjjImZ8i7gGgw0A+4CfhORP9zHVhH5IzbxPDTzQdBiOCXpxqAzxiS5Mgu/qqaoai1Vrauq9UIedcMZp0dE9heRWSLyg4h8LyI3uNMbich0EVns/mwYzV8oLKvmwbdvwjFXQ8MDY755Y4zxk5d3Dg8Cg1V1nntyeK6ITAf6AzNUdaiIDAGGALd5FWLBqi30HzItZIryVvqDHCz16DKjPVtnTCtzWWOMqY48u8OIqq5W1Xnu8604N1pvjnPV7yh3tlFAD68ylOb0lDkcnfIjTwd7stVG3zTGJCFR1YrnqupGRFoCHwPtgf+qagN3ugCbSl7vscxAYCBAZmZm59zc3Epte+3GLazZ4TxP0SB9fr2dIknlzQMeQiVQqXVGQ2YtdueKJ5YrMpYrMpYrMq3qB8jIyKj08l26dJmrqtl7TveyqQcAEckAxgF/V9U/JKTbpKqqiJT6yaOqw3F6FJGdna05OTmV2v7zoycxbIHza14WmE6DtDX0L7iVvIU1KrW+aBncIbg7VzyxXJGxXJGxXJEZ2bUOla195fH0ZrIikoZT9Eer6nh38hoRaeq+3xRY62WGEvXJ5/rU8Xxc1IG84sNjsUljjIlLnhV+txnnFWCRqj4V8tZkoJ/7vB8wyasMoa5PnUBdtvNwsA821r4xJpl5+d3meKAvsEBE5rvT7gCGAm+LyADgV+ACDzMA0FJW87fAh7xVlMNPeoDXmzPGmLjmWeFX1U8p+9D6FK+2W5ohqbnsIo2ng71iuVljjIlLnrbxx4Nm23+ka+Br/hnsxjoa+B3HGGN8V70Lf3ExJ64fwyptzL+LzvI7jTHGxIXqXfgXvM2+u5bzeOGF7CLd7zTGGBMXqnfhnz+GNTVaMbn4OL+TGGNM3Kjehf+ScUxregNazX9NY4yJRPWuiIE0tqU18juFMcbElepd+I0xxuzFCr8xxiQZK/zGGJNkrPAbY0ySscJvjDFJxgq/McYkGSv8xhiTZKzwG2NMkrHCb4wxScYKvzHGJBkvb704QkTWisjCkGmNRGS6iCx2fzb0avvGGGNK5+UR/0ig6x7ThgAzVLUNMMN9bYwxJoY8K/yq+jGwcY/J3YFR7vNRQA+vtm+MMaZ0oqrerVykJTBVVdu7rzeragP3uQCbSl6XsuxAYCBAZmZm59zc3EplWLtxC2t2VGpRT2XWwnJFwHJFxnJFJl5ztaofICMjo9LLd+nSZa6qZu853bObrVdEVVVEyvzUUdXhwHCA7OxszcnJqdR2nh89iWELfPs1yzS4Q9ByRcByRcZyRSZec43sWofK1r7yxLpXzxoRaQrg/lwb4+0bY0zSi3Xhnwz0c5/3AybFePvGGJP0vOzO+SbwBXCoiKwUkQHAUOA0EVkMnOq+NsYYE0OeNWqp6kVlvHWKV9s0xhhTMbty1xhjkowVfmOMSTJW+I0xJslY4TfGmCRjhd8YY5KMFX5jjEkyVviNMSbJWOE3xpgkY4XfGGOSjBV+Y4xJMlb4jTEmyVjhN8aYJGOF3xhjkowVfmOMSTJW+I0xJslY4TfGmCTjS+EXka4i8pOILBGRIX5kMMaYZBXzwi8iAeAfwJnAYcBFInJYrHMYY0yy8uOI/yhgiaouU9UCIBfo7kMOY4xJSqKqsd2gSE+gq6pe7r7uCxytqoP2mG8gMNB9eSjwUyU3uQ+wvpLLeslyRcZyRcZyRaa65jpQVZvsOdGzm61XlaoOB4ZXdT0iMkdVs6MQKaosV2QsV2QsV2SSLZcfTT2rgP1DXrdwpxljjIkBPwr/10AbEWklIulAb2CyDzmMMSYpxbypR1WDIjII+AAIACNU9XsPN1nl5iKPWK7IWK7IWK7IJFWumJ/cNcYY4y+7ctcYY5KMFX5jjEky1aLwi0gvEfleRIpFpMyuT2UNFeGeaP6PO/0t96RzNHI1EpHpIrLY/dmwlHm6iMj8kMdOEenhvjdSRH4JeS8rVrnc+YpCtj05ZLqf+ytLRL5w/97ficiFIe9FdX9VNLSIiNRwf/8l7v5oGfLe7e70n0TkjKrkqESum0TkB3f/zBCRA0PeK/VvGqNc/UVkXcj2Lw95r5/7d18sIv1inOvpkEw/i8jmkPc82V8iMkJE1orIwjLeFxF5zs38nYh0Cnmv6vtKVRP+AfwfzkVeeUB2GfMEgKXAQUA68C1wmPve20Bv9/lLwNVRyvU4MMR9PgR4rIL5GwEbgdru65FATw/2V1i5gPwypvu2v4BDgDbu82bAaqBBtPdXef9eQua5BnjJfd4beMt9fpg7fw2glbueQAxzdQn5N3R1Sa7y/qYxytUfeKGUZRsBy9yfDd3nDWOVa4/5r8PpcOL1/voL0AlYWMb7ZwHvAQIcA/wnmvuqWhzxq+oiVa3oyt5Sh4oQEQFOBsa6840CekQpWnd3feGutyfwnqpuj9L2yxJprt383l+q+rOqLnaf/wasBfa6MjEKwhlaJDTvWOAUd/90B3JVdZeq/gIscdcXk1yqOivk39CXONfKeK0qQ7GcAUxX1Y2qugmYDnT1KddFwJtR2naZVPVjnIO8snQHXlPHl0ADEWlKlPZVtSj8YWoOrAh5vdKd1hjYrKrBPaZHQ6aqrnaf/w5kVjB/b/b+R/ew+1XvaRGpEeNcNUVkjoh8WdL8RBztLxE5CucobmnI5Gjtr7L+vZQ6j7s/tuDsn3CW9TJXqAE4R44lSvubxjLX+e7fZ6yIlFzIGRf7y20SawXMDJns1f6qSFm5o7Kv4nbIhj2JyEfAfqW8daeqTop1nhLl5Qp9oaoqImX2nXU/zTvgXN9Q4nacApiO05/3NuCBGOY6UFVXichBwEwRWYBT3CotyvvrdaCfqha7kyu9v6ojEbkEyAZOCpm8199UVZeWvoaomwK8qaq7RORKnG9LJ8do2+HoDYxV1aKQaX7uL88kTOFX1VOruIqyhorYgPM1KtU9aotoCInyconIGhFpqqqr3UK1tpxVXQBMUNXCkHWXHP3uEpFXgZtjmUtVV7k/l4lIHnAEMA6f95eI1AOm4Xzofxmy7krvr1KEM7RIyTwrRSQVqI/z78nLYUnCWreInIrzYXqSqu4qmV7G3zQahazCXKq6IeTlv3HO6ZQsm7PHsnlRyBRWrhC9gWtDJ3i4vypSVu6o7KtkauopdagIdc6YzMJpXwfoB0TrG8Rkd33hrHevtkW3+JW0q/cASu0B4EUuEWlY0lQiIvsAxwM/+L2/3L/dBJz2z7F7vBfN/RXO0CKheXsCM939MxnoLU6vn1ZAG+CrKmSJKJeIHAH8C+imqmtDppf6N41hrqYhL7sBi9znHwCnu/kaAqfz52++nuZys7XFOVn6Rcg0L/dXRSYDf3N79xwDbHEPbKKzr7w4Yx3rB3AuTlvXLmAN8IE7vRnwbsh8ZwE/43xi3xky/SCc/5hLgHeAGlHK1RiYASwGPgIaudOzgX+HzNcS55M8ZY/lZwILcArYG0BGrHIBx7nb/tb9OSAe9hdwCVAIzA95ZHmxv0r794LTdNTNfV7T/f2XuPvjoJBl73SX+wk4M8r/3ivK9ZH7/6Bk/0yu6G8ao1yPAt+7258FtA1Z9jJ3Py4BLo1lLvf1fcDQPZbzbH/hHOStdv8tr8Q5F3MVcJX7vuDcsGqpu+3skGWrvK9syAZjjEkyydTUY4wxBiv8xhiTdKzwG2NMkrHCb4wxScYKvzHGJBkr/MYYk2Ss8BtTRSJSS0Rmi0jAfR0QkWfFGTp6gXu5f+j86SLysXu1rzExZ4XfmKq7DBiv/xvj5XZgmaq2A57DGb55N3VGiZwBXIgxPrDCb5KeiDQTkXEi8o2I/CgiR4lIWxGZKc4NOD5yL9kvuQnGXHeEyU/dVfTBHV5CROoA56rqs+57vwCtS9nsRHc5Y2LOrtw1Sc1tbpmLcyn/VBGpDaQBnwN9VHW+iNwG1AOGAv/BGSKiQEQaANuB/6rqfu76ugMvA7+5m2gEfKSql+2x3QDwu6p6cS8BY8plR/wm2fUAFqnqVAB1bmDSFfhUVee78/wA7AsUAbWAYSKSraqbgX2AzSHrywLuUdUsVc0CPsQZLwcRGVkyk9ssVCAidT35rYwphxV+k+yycO5SFeownIGxSnTAGZl0O9Ae+AwYLiLXADtwBmsr0RDnW0DJt4nTgSnuN4lte2ynBrAzOr+GMeGzwm+S3e9Au5IXItIEZ6TUw9zXBwF9gddEpI2qblPVXGAqUFOd298FRKSk+P+Mc49UgBuBaercfrETMC9kO42B9Rpy/wVjYsUKv0l2I4FMt+vlfOBYnDt7NXPvOJYLXKbOTUTuFJGfRGQezi36XnTX8SFwgvv8TaCTiCwBOgI3udOPxBkbvkQXnJvJGBNzdnLXmCoSkU7Ajarat5x5/g1cWdLlU0TGA0NU9ecYxTRmNzviN6aKVHUeMKvkAq4y5rk8pOinAxOt6Bu/2BG/McYkGTviN8aYJGOF3xhjkowVfmOMSTJW+I0xJslY4TfGmCRjhd8YY5KMFX5jjEkyVviNMSbJ/D/3sArwCTqv7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bin_to_plot = 3\n",
    "number_of_bins_in_hist = 25\n",
    "cos_theta_l_bin = bins[bin_to_plot]['ctl']\n",
    "hist, _bins, _ = plt.hist(cos_theta_l_bin, bins=number_of_bins_in_hist)\n",
    "x = np.linspace(-1, 1, number_of_bins_in_hist)\n",
    "pdf_multiplier = np.sum(hist) * (np.max(cos_theta_l_bin) - np.min(cos_theta_l_bin)) / number_of_bins_in_hist\n",
    "y = d2gamma_p_d2q2_dcostheta(fl=fls[bin_to_plot], afb=afbs[bin_to_plot], cos_theta_l=x) * pdf_multiplier\n",
    "plt.plot(x, y, label=f'Fit for bin {bin_to_plot}')\n",
    "plt.xlabel(r'$cos(\\theta_l)$')\n",
    "plt.ylabel(r'Number of candidates')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Finally, let us visualize the results in each bin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAalElEQVR4nO3df5QddX3/8efL0NQ2IFGTrmySkpy6369irJHvfkOp1mwh2FhtwqmiCYeWVDBQiYXy1Z5YPJyIf5TiKdJjt5qU0uLPlKLFLUbTFrn5+u0BTgIsQpKG3W8Es1ksmLLYRQsmvPvHzMLNenfvvXvvzr0z9/U4Z8/emfkw8/7k7of3zOfzmRlFBGZmZpYvL2t1AGZmZlY/J3AzM7MccgI3MzPLISdwMzOzHHICNzMzyyEncDMzsxzKNIFLWiPpoKRhSVsqbP9FSXdLelDSdyT9ZpbxmZmZ5YWyug9c0hzgUeA8YATYA2yIiP1lZbYDD0bEZySdAeyMiKXT7XfBggWxdOm0RXj22WeZN29eYxVoc51QR3A9y91///0/iIiFGYXUctXauv82isX1fMlUbf2kWYvqp60EhiPiEICkHcA6YH9ZmQBekX4+FRitttOlS5eyd+/eacuUSiX6+vpmEHJ+dEIdwfUsJ+nxbKJpD9Xauv82isX1fMlUbT3LBL4IOFy2PAKcNanMVuCfJH0ImAesrrQjSZuATQBdXV2USqVpDzw+Pl61TN51Qh3B9TQzm5BlAq/FBuBvI+LPJJ0NfF7S8oh4obxQRGwHtgP09vZGtbOXTjiT64Q6gutpZjYhy0lsR4AlZcuL03XlLgFuA4iIe4CXAwsyic7MzCxHskzge4AeScskzQXWAwOTynwPOBdA0utJEvhTGcZoZk1Q7Y6TsnLvlhSSerOMz6wIMkvgEXEM2AzsAg4At0XEPknXSVqbFvs/wAckPQR8GdgYfl2aWa6kd5z0A+8AzgA2pHeVTC53CnAlcF+2EZoVQ6Zj4BGxE9g5ad21ZZ/3A2/JMiYza7pa7jgB+ATwp8BHsg3PrBjabRKbzURfHyvGxmBwsNWRmEENd5xIOhNYEhFflzRlAq/njpNOmbnves6OFVddBcDgTTdldkxorJ5O4GaWKUkvA24ENlYrW88dJ50yc9/1nCXz5wNk/m/bSD39LHQza7Zqd5ycAiwHSpIeA34FGPBENrP6OIGbWbNNe8dJRDwTEQsiYmn6qOR7gbURMf0jFc1m09AQ7N4NV1zR6khq5gRuZk1V4x0nZu1lNH1y97ZtrY2jDh4DN7Omq3bHyaT1fVnEZDat7u4kiV92WasjqZkTuM3MxKSLDpgNa2YdoKcn+envb3UkNXMXupmZWQ45gZuZmeWQE3gRDA1x6kMP5Wr2pJlZWymVcjck6AReBKOjCHI1e9LMzBrjSWxF0N1NjI6iHM2eNDPreA0+BtsJvAh6enhm4ULm52j2pJmZNcZd6DYzOXxqkZlZkTiB28zk8KlFZmZF4gRuM9Pdnfz2uLuZWUtkmsAlrZF0UNKwpC0Vtn9K0mD686iksSzjszr09MCqVbl6apGZWVtp8BbgzBK4pDlAP/AO4Axgg6QzystExB9GxIqIWAF8GvhqVvGZmZllqsFbgLO8Al8JDEfEoYh4HtgBrJum/Abgy5lEZmZmlrXubgJmPBSZ5W1ki4DDZcsjwFmVCko6HVgGfCuDuMzMzLLX4C3A7Xof+Hrg9og4XmmjpE3AJoCuri5KVR5/Nz4+XrVMrm3dyvj4OCdnWMcVY2MADGb871r47zLVKfU0s5nLMoEfAZaULS9O11WyHphyVD8itgPbAXp7e6Nv4tWWUyiVSlQrk3eZ13H+fIDM/1074buE/NdT0hrgz4E5wM0Rcf2k7ZeTtPHjwDiwKSL2Zx6oWY5lmcD3AD2SlpEk7vXAhZMLSXod8Ergngxja55OeU920etnM1Y2YfU8kqGyPZIGJiXoL0XEZ9Pya4EbgTWZB2vWSqUSg6USfTP8zzObxBYRx4DNwC7gAHBbROyTdF3agCesB3ZERGQVm5k1VdUJqxHxw7LFeYDbu1mdMh0Dj4idwM5J666dtLw1y5jMrOlqmrAq6QrgamAucE6lHdUz36VT5g24nsXSSD3bdRKb2U9r8M091l4ioh/ol3Qh8DHg4gplap7vkvd5A7VyPYulkXr6Uapm1mz1TFiFpIv9/NkMyKyInMDNrNlenLAqaS7JvJaB8gKSesoW3wkMZRifWSG4C93MmioijkmamLA6B7hlYsIqsDciBoDNklYDPwGepkL3uZlNzwnczJqu2oTViLgy86DMCsZd6GZmZjnkBN5sQ0Owe/eMXw9n02jw1XtmZkXiBN5so6PJ7xm+Hs6m0eCr98zMisQJvNm6u5PfM3w9nE2jwVfvmZkViSexNVtPT/Izw9fD2TQafPWemVmR+ArczMwsh5zAzczMcsgJ3MzMLIecwM3MzHLIk9gsP0olBksl+lodh5lZG/AVuJmZWQ75CrzZOuAF9GZms6qvjxVjYzA42OpI2lqmV+CS1kg6KGlY0pYpyrxX0n5J+yR9Kcv4zMzM8iKzK3BJc4B+4DxgBNgjaSAi9peV6QE+CrwlIp6W9AtZxWdmZpYnWV6BrwSGI+JQRDwP7ADWTSrzAaA/Ip4GiIgnM4zPzMzagV9cVJMsx8AXAYfLlkeAsyaV+R8Akv4VmANsjYhvTt6RpE3AJoCuri5KVcadx8fHq5bJu06oI7ieeSBpDfDnJG345oi4ftL2q4FLgWPAU8D7I+LxzAO19lX+4iI/OnlK7TaJ7SSgB+gDFgP/V9IbI2KsvFBEbAe2A/T29kZfX9+0Oy2VSlQrk3edUEdwPdtdLUNlwINAb0T8SNLvAzcA78s+Wmtb3d3E6Cjyi4umlWUX+hFgSdny4nRduRFgICJ+EhHfBR4lSehmlg9Vh8oi4u6I+FG6eC/J/wvMXtLTwzNvepOvvqvI8gp8D9AjaRlJ4l4PXDipzB3ABuBvJC0g6VI/lGGMZtaYWobKyl0CfGOqjfUMl+V52KEenVDPFWNjHD9+vPD1hMa+z8wSeEQck7QZ2EUyNnZLROyTdB2wNyIG0m1vl7QfOA58JCKOZhWjmWVH0kVAL7BqqjL1DJflddihXh1Rz8HBzqgnjX2fmY6BR8ROYOekddeWfQ7g6vTHzPKnlqEyJK0GrgFWRcRzGcVmVih+lKqZNdOLQ2WS5pIMlQ2UF5D0ZmAbsNa3iprNnBO4mTVNRBwDJobKDgC3TQyVSVqbFvskcDLw95IGJQ1MsTszm0a73UZmZjlXw1DZ6syDMisgX4GbmZnlkBO4mZlZDjmBm5mZ5ZATuJmZWQ45gZuZmeWQE7iZmVkOOYGbmZnlkBO4mZlZDjmBm5mZ5ZATuJmZWQ45gZuZmeWQE7iZmVkOOYGbmZnlUKYJXNIaSQclDUvaUmH7RklPpa8YHJR0aZbxmXUKSa+XdF76+edaHY+Z1S+z14lKmgP0A+cBI8AeSQMRsX9S0b+LiM1ZxWXWoT4N3Cnpg8AxSQfKX/lpZu0vyyvwlcBwRByKiOeBHcC6DI9vZi/ZHxE3AU9ExAXAq5q14xp62t4m6QFJxyS9p1nHNes0mV2BA4uAw2XLI8BZFcq9W9LbgEeBP4yIw5MLSNoEbALo6uqiVCpNe+Dx8fGqZfKuE+oIrmcTnS3pL4DXSnojoGbstMaetu8BG4EPN+OYZp2q4QQu6dci4tvNCAb4R+DLEfGcpMuAW4FzJheKiO3AdoDe3t7o6+ubdqelUolqZfKuE+oIrmcTvRVYCPwv4ALg9Cbt98WeNgBJEz1tLybwiHgs3fZCk45p1pGacQV+AVBLAj8CLClbXpyue1FEHC1bvBm4oeHozKyST0XEByX9UpPHvmvtaatJPb1t7p0pFtezuroTuKQB4LvAA8D9dexjD9AjaRlJ4l4PXDhp36dFxBPp4lrgQL3xmVld3gfsbnUQU6mnt829M8XielZXdRKbpM+VL0fEWuBG4IckSbimrreIOAZsBnaRJObbImKfpOskrU2L/YGkfZIeAv6AZJzMzJpvpaS/BF4v6U2SmjWhtWpPm5k1Ry1Xz2+c+CDpnyLi7RHxOPA48A/1HCwidgI7J627tuzzR4GP1rNPM6tfRPRKWkwyBv5u4GMkw2GNqtrTZjkzcXXYAd3ZeVNLAo+yzwtnKxAzy1ZEjAAjknZHxFiT9nlM0kRP2xzglomeNmBvRAxI+t8kJ/+vBH5L0scj4g3NOL5ZJ6klgb9G0kbgIZp0q4mZZU/SPOANwPKy38uBecD8Zh2nhp62PSRd62bWgFoS+FaSbrbfAxZLehjYl/7sj4ivzF54ZtYMkh4Dfobkdq5/I5mHsgFYERFPtjA0a3dDQzA6CldcAf39rY7GylSduBIR2yPiQxGxKiIWAO8guT/7eeD8WY7PzJrjH4H/AP4qbc9/CTzn5G1VjY4mv7dta20c9lPqvo1sYtwM+EbzwzGz2RARH5J0OvBxSR8BruXE+S1mlXV3J0n8sstaHYlNkuWjVM2shdK7RzZKegPwCZL5Lb8eEXe3ODRrZz09yY+7z9uO3wdu1mEiYl9E/Dbw68A1ktr2QS5mNjVfgZt1qIi4D1gt6dxWx2Jm9XMCN+twEXFXq2OwNuYHuLQtd6GbmZnlkBO4mZlZDjmBm5mZ5ZATuJmZWQ45gZuZmeWQE7iZmVkOOYGbmZnlUKYJXNIaSQclDUvaMk25d0sKSb1ZxmdmZpYXmSVwSXOAfpK3mZ0BbJB0RoVypwBXAvdlFZuZNVe1k3VJPyvp79Lt90la2oIwzXItyyvwlcBwRByKiOeBHcC6CuU+Afwp8F8ZxmZmTVLjyfolwNMR8VrgUyRt3szqkOWjVBcBh8uWR4CzygtIOhNYEhFfT195WJGkTcAmgK6uLkpVHvU3Pj5etUzedUIdwfXMiRdP1gEkTZys7y8rsw7Ymn6+HfgLSYoIv+LUrEZt8yx0SS8DbgQ2VisbEduB7QC9vb3R19c3bflSqUS1MnnXCXUE1zMnqp6sl5eJiGOSngFeDfygvFA9J+s5P+mpmetZLI3UM8sEfgRYUra8OF034RRgOVCSBPAaYEDS2ojYm1mUZtY26jlZz/lJT81cz2JppJ5ZjoHvAXokLZM0F1gPDExsjIhnImJBRCyNiKXAvYCTt1n+VDtZP6GMpJOAU4GjmURnVhCZJfCIOAZsBnYBB4DbImKfpOskrc0qDjObddOerKcGgIvTz+8BvuXxb7P6ZDoGHhE7gZ2T1l07Rdm+LGIys+ZKx7QnTtbnALdMnKwDeyNiAPhr4POShoH/IEnyZlaHtpnEZmbFUe1kPSL+C7gg67jMisSPUjUzM8shJ3AzM7MccgI3MzPLISdwMzOzHHICNzMzyyEncDMzsxxyAjczM8shJ3AzM7MccgI3MzPLISdwMzOzHHICN5tOX1/yY2bWZor/LPS+PlaMjcHgYKsjMTMza5riJ3Azs6LwBYmVcQI3s/xzYrMO5DFwMzOzHMo0gUtaI+mgpGFJWypsv1zSw5IGJf0/SWdkGZ+ZNUbSqyT9s6Sh9Pcrpyj3TUljku7MOkazosgsgUuaA/QD7wDOADZUSNBfiog3RsQK4AbgxqziM7Om2ALcFRE9wF3pciWfBH4ns6iKYmiIUx96CK64otWRWBvI8gp8JTAcEYci4nlgB7CuvEBE/LBscR4QGcZnZo1bB9yafr4VOL9SoYi4C/jPph21UxLb6CgC2Lat1ZFYG8hyEtsi4HDZ8ghw1uRCkq4ArgbmAudkE5rZFIaGYHQ0SQz9/a2OJg+6IuKJ9PP3ga5GdiZpE7AJoKuri1KpVLHcqjSxxWc/y+4LLmjkkG3t7Fe/mrlHj3LkXe9ieIp/i6IYHx+f8vsukkbq2Xaz0COiH+iXdCHwMeDiyWVqbdQAZz/yCKcePcrI+eczfNVVsxJzO/Af++xoVWJo5+9T0r8Ar6mw6ZryhYgISQ31okXEdmA7QG9vb/RN9VCd7m5idBRdfjlTlimC5csZGxtj8R13sLjVscyyUqlU7O8y1Ug9s0zgR4AlZcuL03VT2QF8ptKGmhs1wNGjycHuvJPFd9xRR7j54j/2WdLdnXRbZpwY2vn7jIjVU22T9O+STouIJySdBjyZSVA9PTyzcCHz3UtiHSTLMfA9QI+kZZLmAuuBgfICknrKFt8JDDV81O7uZCD9sssa3pV1oJ4eWLXK3ee1G+ClXrOLga+1MBazQsvsCjwijknaDOwC5gC3RMQ+SdcBeyNiANgsaTXwE+BpKnSf181n5mZZuh64TdIlwOPAewEk9QKXR8Sl6fK3gdcBJ0saAS6JiF0titkslzIdA4+IncDOSeuuLft8ZZbxmLWlHD9VLCKOAudWWL8XuLRs+deaeuBSicFSib6m7rQNdUo9rSZ+EpuZmVkOOYGbmZnlkBO4mZlZDrXdfeBmbaVN78U2M/MVuJmZWQ45gZuZmeVQ8bvQfduFmZkVkK/AzczMcsgJ3MzMLIecwM3aTae829rMGuIEbtZu0leYsm1bqyOx6fT1JT9mLeIEbtZu/AY9M6tB8Wehm+WN36BnZjXwFbiZmVkOOYGbmc3E0BDs3u3JhtYyTuBm1jSSXiXpnyUNpb9fWaHMCkn3SNon6TuS3teKWBs2Opr89mRDaxEncDNrpi3AXRHRA9yVLk/2I+B3I+INwBrgJknzswuxSbq7k9+ebGgtkmkCl7RG0kFJw5J+qmFLulrS/vSs/C5Jp2cZn5k1bB1wa/r5VuD8yQUi4tGIGEo/jwJPAguzCrBpenpg1SrwZENrkcxmoUuaA/QD5wEjwB5JAxGxv6zYg0BvRPxI0u8DNwD57F4z60xdEfFE+vn7QNd0hSWtBOYC/3+K7ZuATQBdXV2Upnm96/j4+LTbm23F2BgAgxm/cjbreraK61ldlreRrQSGI+IQgKQdJGfrLybwiLi7rPy9wEUZxmdmNZD0L8BrKmy6pnwhIkJSTLOf04DPAxdHxAuVykTEdmA7QG9vb/RN8+CUUqnEdNubbv58gGyPSQvq2SKuZ3VZJvBFwOGy5RHgrGnKXwJ8o9KGes7KoTPO5DqhjtAh9dy6lfHxcU5u03pGxOqptkn6d0mnRcQTaYJ+copyrwC+DlwTEffOUqizq02/H+scbfkgF0kXAb3Aqkrb6zkrh844k+uEOoLrmQMDwMXA9envr00uIGku8A/A5yLi9mzDMyuOLCexHQGWlC0vTtedQNJqkq64tRHxXEaxmVlzXA+cJ2kIWJ0uI6lX0s1pmfcCbwM2ShpMf1a0JFqzHMvyCnwP0CNpGUniXg9cWF5A0puBbcCaiKjY9WZm7SsijgLnVli/F7g0/fwF4AsZh2ZWOJldgUfEMWAzsAs4ANwWEfskXSdpbVrsk8DJwN+nZ+UDWcVnZmaWJ5mOgUfETmDnpHXXln2ecnKMmZmZvcRPYjMzM8shJ3AzM7McUsSUz1nIBUlPAY9XKbYA+EEG4bRSJ9QRXM9yp0dE/h5BOkM1tHX/bRSL6/mSim099wm8FpL2RkRvq+OYTZ1QR3A9bWqd8m/mehZLI/V0F7qZmVkOOYGbmZnlUKck8O2tDiADnVBHcD1tap3yb+Z6FsuM69kRY+BmZmZF0ylX4GZmZoXiBG5mZpZDhU7gktZIOihpWNKWVsczGyQtkXS3pP2S9km6stUxzSZJcyQ9KOnOVscyWyTNl3S7pH+TdEDS2a2Oqd25rReL23mN+yjqGLikOcCjwHnACMnb0DZExP6WBtZkkk4DTouIBySdAtwPnF+0ek6QdDXJu+JfERHvanU8s0HSrcC3I+Lm9N3ZPx8RYy0Oq225rRernuB2XqsiX4GvBIYj4lBEPA/sANa1OKami4gnIuKB9PN/krzpbVFro5odkhYD7wRurlY2rySdSvKu7L8GiIjnnbyrclsvELfz2hU5gS8CDpctj1DAP/ZykpYCbwbua3Eos+Um4I+AF1ocx2xaBjwF/E3ahXizpHmtDqrNua0Xy024ndekyAm8o0g6GfgKcFVE/LDV8TSbpHcBT0bE/a2OZZadBJwJfCYi3gw8CxRyTNdmpsht3e28PkVO4EeAJWXLi9N1hSPpZ0ga9Bcj4qutjmeWvAVYK+kxki7ScyR9obUhzYoRYCQiJq6sbidp6DY1t/XicDuvQ5ET+B6gR9KydILAemCgxTE1nSSRjKMciIgbWx3PbImIj0bE4ohYSvJdfisiLmpxWE0XEd8HDkv6n+mqc4HCTVJqMrf1gnA7r89JTY2qjUTEMUmbgV3AHOCWiNjX4rBmw1uA3wEeljSYrvvjiNjZupCsQR8Cvpgmo0PA77U4nrbmtu62nlMNt/PC3kZmZmZWZEXuQjczMyssJ3AzM7MccgI3MzPLISdwMzOzHHICNzMzyyEn8AKRdFzSoKSHJD0g6VfT9d2Sbm91fACStkr6cKvjMMszt3WDAt8H3qF+HBErACT9BvAnwKqIGAXe08rAmiF9kIUiosjPSDarhdu6+Qq8wF4BPA3Jiw8kPZJ+3ijpq5K+KWlI0g2V/mNJj0n6eHp2/7Ck16XrTzirlvRIuv+l6Xtt/1bSo5K+KGm1pH9Nj7OybPdvknRPuv4DZfv6iKQ9kr4j6eNlsR+U9DngEU58ZKaZua13LF+BF8vPpU9oejlwGnDOFOVWkLzJ6DngoKRPR8ThCuV+EBFnSvog8GHg0irHfy1wAfB+ksdbXgi8FVgL/DFwflrul4FfAeYBD0r6OrAc6CF5NaSAAUlvA76Xrr84Iu6tcnyzTuG2bk7gBVPerXY28DlJyyuUuysinknL7QdO58TXMU6YeFnC/cBv13D870bEw+l+96XHCUkPA0vLyn0tIn4M/FjS3SQN+a3A24EH0zInkzTm7wGPu0GbncBt3ZzAiyoi7pG0AFhYYfNzZZ+PM/XfwXMVyhzjxKGXl0+x3xfKll+YdIzJz+8NkjPxP4mIbeUblLz3+Nkp4jPreG7rnctj4AWVjmPNAY42edePkb72TtKZJC+mr9c6SS+X9Gqgj6QLbhfwfiXvOkbSIkm/0JSIzQrMbb1z+Qq8WCbGxSA5y704Io4nEzqb5ivA76bdZvcBj85gH98B7gYWAJ9IZ86OSno9cE8a7zhwEckVgZmdyG3d/DYyMzOzPHIXupmZWQ45gZuZmeWQE7iZmVkOOYGbmZnlkBO4mZlZDjmBm5mZ5ZATuJmZWQ79N7z0MUykAtQOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 3))\n",
    "ax1.errorbar(np.linspace(0, len(bins) - 1, len(bins)), fls, yerr=fl_errs, fmt='o', markersize=2, label=r'$F_L$', color='red')\n",
    "ax2.errorbar(np.linspace(0, len(bins) - 1, len(bins)), afbs, yerr=afb_errs, fmt='o', markersize=2, label=r'$A_{FB}$', color='red')\n",
    "ax1.grid()\n",
    "ax2.grid()\n",
    "ax1.set_ylabel(r'$F_L$')\n",
    "ax2.set_ylabel(r'$A_{FB}$')\n",
    "ax1.set_xlabel(r'Bin number')\n",
    "ax2.set_xlabel(r'Bin number')\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}