{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic pyerrors example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import pyerrors, as well as autograd wrapped numpy and matplotlib. The sys statement is not necessary if pyerrors was installed via pip."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pyerrors as pe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use numpy to generate some fake data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_sample1 = np.random.normal(2.0, 0.5, 1000)\n",
    "test_sample2 = np.random.normal(1.0, 0.1, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this we can construct `Obs`, which are the basic object of `pyerrors`. For each sample we give to the obs, we also have to specify an ensemble/replica name. In this example we assume that both datasets originate from the same gauge field ensemble labeled 'ens1'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs1 = pe.Obs([test_sample1], ['ens1'])\n",
    "obs2 = pe.Obs([test_sample2], ['ens1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now combine these two observables into a third one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs3 = np.log(obs1 ** 2 / obs2 ** 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pyerrors` overloads all basic math operations, the user can work with these `Obs` as if they were real numbers. The proper resampling is performed in the background via automatic differentiation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are now interested in the error of obs3, we can use the `gamma_method` to compute it and then print the object to the notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Obs[1.415(20)]\n"
     ]
    }
   ],
   "source": [
    "obs3.gamma_method()\n",
    "print(obs3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With print level 1 we can take a look at the integrated autocorrelation time estimated by the automatic windowing procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 1.41522010e+00 +/- 2.03946273e-02 +/- 1.01973136e-03 (1.441%)\n",
      "  t_int\t 5.07378446e-01 +/- 4.51400871e-02 S = 2.00\n"
     ]
    }
   ],
   "source": [
    "obs3.print(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected the random data from numpy exhibits no autocorrelation ($\\tau_\\text{int}\\,\\approx0.5$). It can still be interesting to have a look at the window size dependence of the integrated autocorrelation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcpUlEQVR4nO3de5hcdZ3n8fcnfe8k3R1IAyGJJDOGRQS5bA+KrOMFcYIXMo+ODlFHdFF2LqiDDjvR2QUXd2YYncHLLqNmGRRvoKKO2SUO4wKKC4RNIxgICERASEhIQ0ink753vvPHOd2pdKpPupOuOp2qz+t56qlz+fU53+qkz6d+v3PqlCICMzOziczKuwAzM5vZHBRmZpbJQWFmZpkcFGZmlslBYWZmmRwUZmaWyUFhNgFJ75H0r3nXYZY3B4VVBEm7Cx57JfUVzL/nULYZEd+KiDdNcv/vl/T/DmU/Bdu4WNKvJPVIek7SWklzp7iNyyU9lG7jSUmXH6T9uek+eyXdIemEw3kNVpkcFFYRImLO6AN4GnhbwbJv5V3fwUh6LfA3wMqImAu8DPjOoWwKeB8wD1gOXCrpwgn2OR/4AfBfgaOAzkPcp1U4B4VVNElnSbpH0k5JWyX9T0n16bolkkJSbUH7n0r6YDq9Xy8hbfvHkh5Pt3etEi8DvgycnfZgdh5Cqb8D3BMR9wNExI6IuCEieqaykYj4TET8IiKGI+JR4EfAORM0fzuwMSK+FxH9wKeA0ySddAj1WwVzUFilGwEuA+YDZwPnAn96GNt7K8lB/RXAu4Dfi4hHgD8mOdDPiYi2Q9juvcDvSfpvks6R1FC4UtKqNJyKPoptUJKA1wAbJ9jny4Ffjs5ExB7g1+lyszEOCqtoEXFfRKxL32E/BXwFeO1hbPLqiNgZEU8DdwCnT0OZRMTPSd7hnwncArwg6RpJNen6qyOibaLHBJv9FMnf+FcnWD8H6B63rBuY0nkRq3y1B29iduSSdCJwDdABNJP8n7/vMDa5rWC6l+RgOy0i4sfAjyXNAl4PfA94lCTcpkTSpSTnKl4TEQMTNNsNtIxb1gJMabjLKp97FFbpvgT8ClgWES3AJ0lO+ALsSZ+bC9ofd4j7mbbbMEfE3oi4DbgdOAVA0ifHXdm136Pw5yX9R2AVcG5EbM7Y1UbgtIKfmw38NhMPVVmVclBYpZsL7AJ2pydp/2R0RUR0AVuA90qqSQ+wv32I+3kOWDR6orwYSZ+S9NMJ1q2QdKGkeekJ8rNIhsjWpbX+TeGVXeMfBdt5D8nVU+dFxBMHqfmHwCmS3iGpEbgC2BARv5rKC7fK56CwSvcXwLtJhlP+Fwde/vkh4HLgBZKTuHcf4n5uJ3knvk3S8xO0WQzcNcG6F9NaHicJtm8Cnz2ES3v/O3A0sL6gx/Hl0ZWSNo5+riQNyncAf53u/5VA0UtprbrJX1xkVh6SHiAZDnoh71rMpsJBYWZmmTz0ZGZmmRwUZmaWyUFhZmaZKu4Dd/Pnz48lS5bkXcaRo7c3eW5uzm5nZhXtvvvuez4i2outq7igWLJkCZ2dnXmXceR44IHk+fTT86zCzHIm6TcTrfPQk5mZZXJQmJlZJgeFmZllclCYmVkmB4WZmWVyUJiZWaZcg0LS9ZK2S3pogvXvkbRB0oOS7pZ0WrF2ZmZWOnn3KL4GLM9Y/yTw2og4Ffg0sLocRZmZ2T65fuAuIu6UtCRjfeF3A6wDFpW8KDMz20/ePYqpuBj4cbEVki6R1Cmps6urq8xlmZlVtiMiKCS9niQo/rLY+ohYHREdEdHR3l70ViVmZnaIZvy9niS9ArgOON/fDGZmVn4zukch6SXAD4A/iojH8q7HzKwa5dqjkHQj8DpgvqTNwJVAHUBEfBm4guSL4v9REsBwRHTkU62ZWXXK+6qnlQdZ/0Hgg2Uqx8zMipjRQ09mZpY/B4WZmWVyUJiZWSYHhZmZZXJQmJlZJgeFmZllclCYmVkmB4WZmWVyUJiZWSYHhZmZZXJQmJlZJgeFmZllclCYmVkmB4WZmWVyUJiZWSYHhZmZZXJQmJlZJgeFmZllclCYmVmmXINC0vWStkt6aIL1kvRFSZskbZB0ZrlrNDOrdnn3KL4GLM9Yfz6wLH1cAnypDDWZmVmBXIMiIu4EdmQ0WQF8PRLrgDZJC8pTnZmZQf49ioNZCDxTML85XbYfSZdI6pTU2dXVVbbizMyqwUwPikmJiNUR0RERHe3t7XmXY2ZWUWZ6UGwBFhfML0qXmZlZmcz0oFgDvC+9+ulVQHdEbM27KDOzalKb584l3Qi8DpgvaTNwJVAHEBFfBtYCbwY2Ab3AB/Kp1MyseuUaFBGx8iDrA/izMpVjZmZFzPShJzMzy5mDwszMMjkozMwsk4PCzMwyOSjMzCyTg8LMzDI5KMzMLJODwszMMjkozMwsk4PCzMwyOSjMzCyTg8LMzDI5KMzMLJODwszMMjkozMwsk4PCzMwyOSjMzCyTg8LMzDI5KMzMLFOuQSFpuaRHJW2StKrI+pdIukPS/ZI2SHpzHnWamVWz3IJCUg1wLXA+cDKwUtLJ45r9F+C7EXEGcCHwj+Wt0szM8uxRnAVsiognImIQuAlYMa5NAC3pdCvwbBnrMzMz8g2KhcAzBfOb02WFPgW8V9JmYC3w4WIbknSJpE5JnV1dXaWo1cysas30k9krga9FxCLgzcA3JB1Qc0SsjoiOiOhob28ve5FmZpUsz6DYAiwumF+ULit0MfBdgIi4B2gE5pelOjMzA/INivXAMklLJdWTnKxeM67N08C5AJJeRhIUHlsyMyuj3IIiIoaBS4FbgUdIrm7aKOkqSRekzT4OfEjSL4EbgfdHRORTsZlZdarNc+cRsZbkJHXhsisKph8Gzil3XWZmts9MP5ltZmY5c1CYmVkmB4WZmWVyUJiZWSYHhZmZZXJQmJlZJgeFmZllclCYmVmmXD9wZ2Zmh+9zP3mML9z2+Nj8R89dxmXnnTht23ePwszsCHfZeSfy1NVvAeCpq98yrSEBDgozMzsIDz2ZmU2TUg8B5cU9CjOzaVLqIaC8OCjMzCyTg8LMzDI5KMzMLJNPZpuVQaWe5JyJ/Luefg4KszK47LwTuey8E1my6paxk53lUI0Hzbx+15XMQWFWwfI8aFZjSFWqXM9RSFou6VFJmyStmqDNuyQ9LGmjpG+Xu0YzOzSVeqloNcqtRyGpBrgWOA/YDKyXtCYiHi5oswz4BHBORLwo6Zh8qrVK4He4Zocmzx7FWcCmiHgiIgaBm4AV49p8CLg2Il4EiIjtZa7RKojf4ZodmjyDYiHwTMH85nRZoROBEyXdJWmdpOXFNiTpEkmdkjq7urpKVK6ZWXWa6Seza4FlwOuARcCdkk6NiJ2FjSJiNbAaoKOjI8pco5lZWQ2N7KWnf5ie/iF6+ofZ1T/E7v5hAB7c3M0xLQ0cPbue2prp6QscNCgknRMRdx1s2SHYAiwumF+ULiu0Gbg3IoaAJyU9RhIc6w9z35YjnyuwatY/NLLfQX5semB43PKJ1w+NBHMba5NHQx1zGmtpaUwO53/5/Q107R5gZ+8grU31tM9t4Jj0MTbd0lgw30hTfU1mzZPpUfwP4MxJLJuq9cAySUtJAuJC4N3j2vwzsBL4qqT5JENRTxzmfi1nvs7dKsXA8Ajbdw2wtbufbbv62dbdx7buAQAu+XpncmAfSA74u9ODPrDvIN9Yx9zGWuY07JtuaazluJZGlh0zZ9/6dPnofFNdDZIOqGfJqltY+9HXADA8spcdewbZ3jNAV88A23v62b5rgF937WHdEzuS+Z4BtvcM0HCQnseEQSHpbODVQLukjxWsagGy42cSImJY0qXAren2ro+IjZKuAjojYk267k2SHgZGgMsj4oXD3bdZqezdG7zYOzj2B7h9V//YH2pXT3IAWfX9DekBoW6/g8W+A0Xd2IFkuoYObOp6B4fZ1t3Ptu7+giAYne5jW3c/3X1DHDO3keNak8eClkYWtDYC8PYzFxYNgsa6wz58Tkptzayk59DSmNkuItjVN0zbVRnbyvj5emBO2mZuwfJdwB9MutrsAtcCa8ctu6JgOoCPpQ+z3AyP7OWFPYNs35W8M3sufU7CYICudPr53QPMbqhNu/pp976lgcVHNXPmCfO45cGtnLqodWwIoev5gbHpXf37Dz3sHhimvmbWAe8+C4NkopBx2EwsIugZGN530O/uS5/3D4O+oREWtDZyXHrwP661iZceM4f/sGz+2LKj5zRQM+vAd/Z/vfYRlp+yIIdXN3WSaG2uy2wzYVBExM+An0n6WkT8ZrqLM5sJRocOknf9+w782/ebTsZ725rr0/HdhrEgOOm4ufzusvn7jfk21E78jvEjN97Pe155wqRqiwh6B0fYVTBmPT5MevqHeD4jbHr6h8bqedPnfkZrUx2tTXW0pM9tTfW0NtXS2lyXrqsfa9PaVEd97cwPmb17g/7hEfYMjNA7OMyegRH6hpIhnm+u+81+vYCt3f08191PAAtaG1nQ2pT0BFobOXVRK29qPZZjW5Ll85rrig7vVKPJnKNokLQaWFLYPiLeUKqirDw+95PH+Mm3fw7Aw8duOaJPKo/sDXr6h+juSx67+obHpseWpeshOWhu7xlgz8Aw8+fsO7k3GgBnvGTevl7BNF9BMlmSmN1Qy+yGWha0Hto2IoI9gyOccuWtfHHlGezsLfwdJc+/7tpzwO9qdH197awi4bIvSPYFzL5HW3M9LUV6MhHBwPBe9gwM0zs4Qu/gCHsGh+lND/D7z6cH/YL5PaNtBobpG9oXDH1DIzTUzmJ2fS3NDTU01yXPABuf7ea4liZeufQojm0d7Rk0Mreh1iEwBZMJiu8BXwauIzlPYNMozyuALjvvRC5r7+XNX/j5jDipPDSyd78D2PgD2q7+Ybp7ix/8ewdHmF1fM3YwG/9oaapj4bwmWpvquGXDVr5w4RkcM7eBec31zCoydFApJDGnIfkzP+m4lin97GjIdPcNsbN3cL9/i9HA2batv2jI9PQP01RXM3YlzqlX3sqewWFqa2Yxu76G5vpaZjfU0FRfOzbfXF/D7IZ0XX0Nbc31LJzXNDbfVF/D7Iak3eiy5obkxG6x4Z8lq27hb9/+isP/JdqkgmI4Ir5U8kqqVKVfARQRvLBnkK07+9mys49nRx/dfQAs//ydYweXgeG9tDTWHvAOdvT5qNn1LJ0/e9/6xn1BMKextujBopiP3Hg/L1swtYNmNRoNmTkNtSxsa5rSz+7dG+weTIL9NZ+5g7s+8Qaa62p8vuQINZmg+N+S/hT4ITAwujAidpSsKjti9A2O8Gx3QQDs7B8LgtHp5voaFrQ2cXxbEwvbGjm+rYnTFrex9sFtXPOu08eGL2bXF7/kz448s2aJlsYkzIGxZzsyTSYoLkqfLy9YFsBvTX85lW1kbxwwrDL6APjO+qdpqk+60s1pV/vA6cm/c56Oert6Bg4IgtGewdbufnYPDHN8a3LwP76tieNbG/n3J8zjbW3Hp8saaa4v/t/swzfez8nH+5292Ux30KCIiKXlKORIMTyyNxkrn+DkX3ff0AHj6KPreodGmNNQW3T8HKDzqRfpGxqhLz3RNzY9NEzfYDLdNzRCbc2s/QKkOQ2RpvpamuvSUKmvmWA6adNcX0NjfQ31OwbpqW/im+t+w9bu/YNg+66BZGy/bV8QLD6qmVf91lFjPYSjZ1f2GL+ZZX/g7g0RcbuktxdbHxE/KF1Z+fjn+7ew+cXeIiEwPBYCvYPDzG0sfrK0tamOec31LDl69rirQJL1cxtqJzyo3vj/n+az7zztoDWOXjmSBMgIfYPD9A3uTa4aGRqhPw2ZwunuviG2dfcX/MzI2PRgzx52NLfy0JZujm9r4pyXzuf4tkYWtiWXDWZd6mlm1SGrR/Fa4HbgbUXWBVBRQfGHX7mHe5/cd9pl6dHNXL78pAPCIOtgXw6SaKyrobGuhnnTscEHHuDNX3iUq9/xzunYmplVoKwP3F2ZPn+gfOXk5zv/6ey8SzAzm5Emc/fYK4otj4iMO4OYWaHCz8ssWXXLEf3hRqs+k7nqaU/BdCPwVuCR0pRTfXwAKZ88f9ejn5cxOxJN5qqnfyicl/T3JHd1tWlQjQeQvA7Y1fi7tupQ6r+pQ/mGu2aSLxkyOyQ+YJePe6zlValvgiZzjuJBkqucIPneiHbA5ycqQOFNAX0QqUwO5fKq1N+3kq98yGggFd4TeRh4LiKGS1rVYejo6IjOzs68yzhyPPBA8nz66XlWYRUorxte+qt2D42k+yKio+i6gwVFwUaOITmZDUBEPD095U0vB8UUOSjMjOygOOitHCVdIOlx4EngZ8BTwI+ntUIzM5uxJnPP308DrwIeS+/7dC6wrqRVmZnZjDGZoBiKiBeAWZJmRcQdQNHuyVRJWi7pUUmbJK3KaPcOSSFpWvZrZmaTN5nLY3dKmgPcCXxL0nZg9+HuWFINcC1wHrAZWC9pTUQ8PK7dXOCjwL2Hu08zM5u6yfQofgn0ApcB/wL8GvjVNOz7LGBTRDwREYPATcCKIu0+Dfwd0D8N+zQzsymaTFC8PiL2RsRwRNwQEV8Efmca9r0QeKZgfnO6bIykM4HFEXFL1oYkXSKpU1JnV1fXNJRmZmajJgwKSX+SftjuJEkbCh5PAhtKXZikWcA1wMcP1jYiVkdER0R0tLe3l7o0M7OqknWO4tskl8H+LVB4orlnmr4vewuwuGB+Ubps1FzgFOCn6fcoHweskXRBRPiDEmZmZZL1fRTdQDewskT7Xg8sk7SUJCAuBN49bv/zR+cl/RT4C4eEmVl5TeYcRUmktwG5lOROtI8A342IjZKuknRBXnWZmdn+DuXusdMmItYCa8ctm+iLkl5XjprMzGx/ufUozMzsyOCgMDOzTA4KMzPL5KAwM7NMDgozM8vkoDAzs0wOCjMzy+SgMDOzTA4KMzPL5KAwM7NMDgozM8vkoDAzs0wOCjMzy+SgMDOzTA4KMzPL5KAwM7NMDgozM8vkoDAzs0wOCjMzy5RrUEhaLulRSZskrSqy/mOSHpa0QdJtkk7Io04zs2qWW1BIqgGuBc4HTgZWSjp5XLP7gY6IeAVwM/CZ8lZpZmZ59ijOAjZFxBMRMQjcBKwobBARd0REbzq7DlhU5hrNzKpenkGxEHimYH5zumwiFwM/LrZC0iWSOiV1dnV1TWOJZmZ2RJzMlvReoAP4bLH1EbE6IjoioqO9vb28xZmZVbjaHPe9BVhcML8oXbYfSW8E/gp4bUQMlKk2MzNL5dmjWA8sk7RUUj1wIbCmsIGkM4CvABdExPYcajQzq3q5BUVEDAOXArcCjwDfjYiNkq6SdEHa7LPAHOB7kh6QtGaCzZmZWYnkOfRERKwF1o5bdkXB9BvLXpSZme3niDiZbWZm+XFQmJlZJgeFmZllclCYmVkmB4WZmWVyUJiZWSYHhZmZZXJQmJlZJgeFmZllclCYmVkmB4WZmWVyUJiZWSYHhZmZZXJQmJlZJgeFmZllclCYmVkmB4WZmWVyUJiZWSYHhZmZZco1KCQtl/SopE2SVhVZ3yDpO+n6eyUtyaFMM7OqlltQSKoBrgXOB04GVko6eVyzi4EXI+KlwOeAvytvlWZmlmeP4ixgU0Q8ERGDwE3AinFtVgA3pNM3A+dKUhlrNDOrenkGxULgmYL5zemyom0iYhjoBo4evyFJl0jqlNTZ1dVVonLNzKpTRZzMjojVEdERER3t7e15l2NmVlHyDIotwOKC+UXpsqJtJNUCrcALZanOzMyAfINiPbBM0lJJ9cCFwJpxbdYAF6XTfwDcHhFRxhrNzKpebV47johhSZcCtwI1wPURsVHSVUBnRKwB/gn4hqRNwA6SMDEzszLKLSgAImItsHbcsisKpvuBd5a7LjMz26ciTmabmVnpOCjMzCyTg8LMzDI5KMzMLJODwszMMjkozMwsk4PCzMwyOSjMzCyTg8LMzDI5KMzMLJODwszMMjkozMwsk4PCzMwyOSjMzCyTg8LMzDI5KMzMLJODwszMMjkozMwsk4PCzMwy5RIUko6S9BNJj6fP84q0OV3SPZI2Stog6Q/zqNXMrNrl1aNYBdwWEcuA29L58XqB90XEy4HlwOcltZWvRDMzg/yCYgVwQzp9A/D74xtExGMR8Xg6/SywHWgvV4FmZpbIKyiOjYit6fQ24NisxpLOAuqBX0+w/hJJnZI6u7q6prdSM7MqV1uqDUv6v8BxRVb9VeFMRISkyNjOAuAbwEURsbdYm4hYDawG6OjomHBbZmY2dSULioh440TrJD0naUFEbE2DYPsE7VqAW4C/ioh1JSrVzMwy5DX0tAa4KJ2+CPjR+AaS6oEfAl+PiJvLWJuZmRXIKyiuBs6T9DjwxnQeSR2SrkvbvAv4XeD9kh5IH6fnUq2ZWRUr2dBTloh4ATi3yPJO4IPp9DeBb5a5NDMzG8efzDYzs0wOCjMzy+SgMDOzTA4KMzPL5KAwM7NMiqisDzJL6gJ+cxibmA88P03lHCmq7TVX2+sFv+ZqcTiv+YSIKHo/vYoLisMlqTMiOvKuo5yq7TVX2+sFv+ZqUarX7KEnMzPL5KAwM7NMDooDrc67gBxU22uuttcLfs3VoiSv2ecozMwsk3sUZmaWyUFhZmaZHBQpScslPSppk6RVeddTapIWS7pD0sOSNkr6aN41lYukGkn3S/o/eddSDpLaJN0s6VeSHpF0dt41lZqky9L/1w9JulFSY941TTdJ10vaLumhgmVHSfqJpMfT53nTsS8HBcmBA7gWOB84GVgp6eR8qyq5YeDjEXEy8Crgz6rgNY/6KPBI3kWU0ReAf4mIk4DTqPDXLmkh8BGgIyJOAWqAC/OtqiS+Biwft2wVcFtELANuS+cPm4MicRawKSKeiIhB4CZgRc41lVREbI2IX6TTPSQHj4X5VlV6khYBbwGuO1jbSiCpleQLwP4JICIGI2JnrkWVRy3QJKkWaAaezbmeaRcRdwI7xi1eAdyQTt8A/P507MtBkVgIPFMwv5kqOGiOkrQEOAO4N+dSyuHzwH8G9uZcR7ksBbqAr6bDbddJmp13UaUUEVuAvweeBrYC3RHxr/lWVTbHRsTWdHobcOx0bNRBUeUkzQG+D/x5ROzKu55SkvRWYHtE3Jd3LWVUC5wJfCkizgD2ME3DETNVOi6/giQkjwdmS3pvvlWVXySffZiWzz84KBJbgMUF84vSZRVNUh1JSHwrIn6Qdz1lcA5wgaSnSIYX3yCp0r9udzOwOSJGe4s3kwRHJXsj8GREdEXEEPAD4NU511Quz0laAJA+b5+OjTooEuuBZZKWSqonOfG1JueaSkqSSMatH4mIa/Kupxwi4hMRsSgilpD8G98eERX9TjMitgHPSPp36aJzgYdzLKkcngZeJak5/X9+LhV+Ar/AGuCidPoi4EfTsdHa6djIkS4ihiVdCtxKcoXE9RGxMeeySu0c4I+AByU9kC77ZESsza8kK5EPA99K3wQ9AXwg53pKKiLulXQz8AuSq/vupwJv5yHpRuB1wHxJm4ErgauB70q6mOTrFt41LfvyLTzMzCyLh57MzCyTg8LMzDI5KMzMLJODwszMMjkozMwsk4PCrEQkfU7SnxfM3yrpuoL5f5D0sVyKM5sCB4VZ6dxF+olgSbOA+cDLC9a/Grg7h7rMpsRBYVY6dwOj3/3wcuAhoEfSPEkNwMtIPhRmNqP5k9lmJRIRz0oalvQSkt7DPSR3JT4b6AYeTG9rbzajOSjMSutukpB4NXANSVC8miQo7sqxLrNJ89CTWWmNnqc4lWToaR1Jj8LnJ+yI4aAwK627gbcCOyJiJCJ2AG0kYeGgsCOCg8KstB4kudpp3bhl3RHxfD4lmU2N7x5rZmaZ3KMwM7NMDgozM8vkoDAzs0wOCjMzy+SgMDOzTA4KMzPL5KAwM7NM/wZOuOjl0f0peAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This figure shows the windowsize dependence of the integrated autocorrelation time. The red vertical line signalizes the window chosen by the automatic windowing procedure with $S=2.0$.\n",
    "Choosing a larger windowsize would not significantly alter $\\tau_\\text{int}$, so everything seems to be correct here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlated data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now generate fake data with given covariance matrix and integrated autocorrelation times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov = np.array([[0.5, -0.2], [-0.2, 0.3]]) # Covariance matrix\n",
    "tau = [4, 8] # Autocorrelation times\n",
    "c_obs1, c_obs2 = pe.misc.gen_correlated_data([2.8, 2.1], cov, 'ens1', tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and once again combine the two `Obs` to a new one with arbitrary math operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 3.27194697e-01 +/- 1.96872835e+00 +/- 3.38140198e-01 (601.699%)\n",
      "  t_int\t 5.41336983e+00 +/- 1.59801329e+00 S = 2.00\n"
     ]
    }
   ],
   "source": [
    "c_obs3 = np.sin(c_obs1 / c_obs2 - 1)\n",
    "c_obs3.gamma_method()\n",
    "c_obs3.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time we see a significant autocorrelation so it is worth to have a closer look at the normalized autocorrelation function (rho) and the integrated autocorrelation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAApV0lEQVR4nO3deXhV5bn+8e+TiTlACJNMYQhVcGCIgAotKioqSkdFq1XrcKzV9rT9eURrpbW1x1orta22WLUOdazHAQWllDogKhIKMsuMEIYwhkAIZHh+f+ydGGKyybCTlex9f65rX8kasvazJO477/uu9S5zd0RERKqTEHQBIiLStCkoREQkIgWFiIhEpKAQEZGIFBQiIhKRgkJERCJSUEjcM7MnzOxXQdch0lQpKCTmmdlGMztkZgfMbHs4GNoGVMsdZrYhXMsWM3uhDse438zWmFm+ma0ys+8cY//LzWyTmR00s1fNLK3uZyDxSEEh8eIid28LDAGGArc3dgFmdhVwJTAuXEsWMKcOhzoIXAS0B64CHjSz06t5z8HAtPD7dgUKgIfr8J4SxxQUElfcfTswi1BgVNTRzGaE/0qfb2b9yzaY2elmtsDM8sJfq/xQroFTgVnuvq6sFnd/pA7nMMXdV7l7qbvPB+YCp1Wz+7eB1939PXc/APwM+LqZtavjOUgcUlBIXDGznsD5wNpKmyYBvwA6hrfdE94/DZgB/AHoBDwAzDCzTnV4+4+A75jZrWaWZWaJlWp72Mz2VfNaUs35tCIUQMurec/BwCdlC+GQOgIMrEP9EqcUFBIvXjWzfGAzkAtMqbT9FXf/2N2LgWf4vMVxIbDG3Z9292J3fw5YRajrp1bc/e/ALcB5wLtArpndVmH7Te7eoZrXydUc9i+EgmBWNdvbAnmV1uUBalFIjSkoJF581d3bAWOB44H0Stu3V/i+gNAHLMBxwKZK+24CetSlCHd/xt3HAR2AG4Ffmtl5dTmWmf0WOBG4xKuf3fMAkFppXSqQX5f3lPikoJC44u7vAk8A99fwR7YCfSqt6w3k1LOOInf/B7CE0Ic9ZvaX8NVQVb2O6loys18Q6kI71933R3ir5cApFX6uH9ACWF2f+iW+KCgkHv0eOMfMTjnWjsBMYGD4EtMkM7sUGAS8UdXO4Utvn6hm29VmdqGZtTOzBDM7n9AYwnwAd7/R3dtW8xpc4Ti3A5cTunpq9zHqfwa4yMzGmFkb4G7gZXdXi0JqTEEhccfddwJPAXfVYN/dwATgJ8Bu4H+ACe6+q5of6QXMq2bbfuAO4DNgH3Af8D13f7829QO/JtSqWVuhxXFH2cbw8phw/csJdXE9Q2hsph1wUy3fT+Kc6cFFItFhZimEBpZPdveioOsRiRYFhYiIRKSuJxERiUhBISIiESkoREQkoqSgC4i29PR0z8jICLoMkcgKCkJfW7cOtg6RsIULF+5y985VbYu5oMjIyCA7OzvoMkQiW7w49HXIkCCrEClnZpVnICinricREYlIQSEiIhEpKEREJCIFhYiIRKSgEBGRiAINCjN73MxyzWxZNdvNzP5gZmvNbImZDWvsGkVE4l3QLYongPERtp8PZIZfNwB/boSaRESkgkCDwt3fA/ZE2GUi8JSHfAR0MLPuNTl2wZFiFm/eF4UqRUTiW1O/4a4HoWccl9kSXret4k5mdgOhFge9e/fm0mkfMn/D5/kzsm8aL/zXaQ1frYhIDGrqQVEj7v4I8AhAVlaWl4XCmPv+zd+uPpUBXfQceRGRugp6jOJYcgg9MaxMT2rxrOKRfTvx4fpIPVsiInIsTT0opgPfCV/9NArIc/dtx/qhMqP6deKj9cd6pLCIiEQSaNeTmT0HjAXSzWwLMAVIBnD3vxB6sP0FwFqgALimNscf2TeNe99cibtjZtEsXUQkbgQaFO5+2TG2O/D9uh6/V1prWiQlsm7nQQZ0aVvXw4iIxLWm3vVUb+p+EhGpnzgIijQFhYhIPcRBUHRi/oY9hHqxRESktmI+KHp2bEVKYgLrdx0MuhQRkWYp5oPCzBip7icRkTqL+aAAGNW3Ex/pxjsRkTqJj6Do14n563drnEJEpA7iIih6pbUiKcHYoHEKEZFai4ugCI1TqPtJRKQu4iIoIHQ/xfwNGtAWEamtOAqK0B3aGqcQEamduAmK3mmtSTBj4+6CoEsREWlW4iYozIyRfXU/hYhIbcVNUADkFxZz+8tLyZg8g4zJM5g6e3XQJYmINHkx8SjUmvrZhEEs37qf7fsL2XjvhUGXIyLSLMRVi6JPp9Y4GswWEamNuAoKM2NUv05BlyEi0qzEVVAACgoRkVqKu6AYPSAdgNJSdUGJiNRE3AVFr7TWAKzanh9wJSIizUPcBUWZuWt2Bl2CiEizEFdBMXX2ajImzwDgf99cpfsoRERqIK6C4kfnDGTjvRey5Ofn0iYlke+N7R90SSIiTV6gQWFm483sUzNba2aTq9je28zeNrNFZrbEzC6Ixvumtkzm+O6pfLxB046LiBxLYEFhZonAQ8D5wCDgMjMbVGm3O4EX3X0oMAl4OFrvPyYznffX7orW4UREYlaQLYoRwFp3X+/uR4DngYmV9nEgNfx9e2BrtN58TGZn3lutAW0RkWMJMih6AJsrLG8Jr6vo58AVZrYFmAncUtWBzOwGM8s2s+ydO2v24X9Kz/Zs3XeI3PzCWhcuIhJPmvpg9mXAE+7eE7gAeNrMvlCzuz/i7lnuntW5c+caHTgpMYHT+ndinrqfREQiCjIocoBeFZZ7htdVdC3wIoC7fwi0BNKjVcDozM7MXaOgEBGJJMigWABkmllfM0shNFg9vdI+nwFnA5jZCYSCImoDC1/OTOf9Nbv0eFQRkQgCCwp3LwZuBmYBKwld3bTczO42s4vDu/0EuN7MPgGeA672KH6q9+nUhhbJCazecSBahxQRiTmBPrjI3WcSGqSuuO6uCt+vAM5oyBrGZHbmVzNWHNUF9cOzM/nROQMb8m1FRJqNpj6Y3eDGDEgnwaz8iXcb771QISEiUkHcB8Xp/dPJ3riHwqKSoEsREWmS4j4o2rdOJrNrOxZu2ht0KSIiTVLcBwWErn7SZbIiIlVTUABjBnbW8ylERKqhoACG9OrAZ7sLgi5DRKRJUlAAyYkJjOzXKegyRESaJAVF2PgTuwVdgohIk6SgCDtvcFcA8gqKAq5ERKRpUVCEtWuZDMBby7cFXImISNOioKhk+idRezaSiEhMUFBUsnRLnh5mJCJSgYKiknEndGXGEnU/iYiUUVBUctGQ49T9JCJSQaDTjDcVU2ev5sE5awC45m8LaJWcyOY9Bby0cEv5etD04yISnyzWnu6WlZXl2dnZ9TrGT19ZynEdWvH9MwcAkDF5Rvk05CJRsXhx6OuQIUFWIVLOzBa6e1ZV29T1VIWLTzmO19X9JCICKCiqdGpGGnmHili9Iz/oUkREAqegqEJCgjHh5O5MX6xWhYiIgqIaF5/Sg+mfbCXWxnBERGpLQVGNE3ukkphgfLIlL+hSREQCpaCohplx0SnHqftJROKegiKCi085jjeWKChEJL4FGhRmNt7MPjWztWY2uZp9LjGzFWa23Myebcz6BnRpS5fUFo35liIiTU5gQWFmicBDwPnAIOAyMxtUaZ9M4HbgDHcfDPx3Y9d5+Yg+jf2WIiJNSpAtihHAWndf7+5HgOeBiZX2uR54yN33Arh7biPXyMQhxwGwLe9QY7+1iEiTEGRQ9AA2V1jeEl5X0UBgoJnNM7OPzGx8o1UX1qZFaDqs5z4OlTp19moyJs8of02dvbqxSxIRaVRNfTA7CcgExgKXAX81sw6VdzKzG8ws28yyd+7c2SCFvLDgM4pKSvnROQPL533aeO+FmiRQRGJekEGRA/SqsNwzvK6iLcB0dy9y9w3AakLBcRR3f8Tds9w9q3Pnzg1SbK+OrZmzckeDHFtEpCkLMigWAJlm1tfMUoBJwPRK+7xKqDWBmaUT6opa34g1lrtiVB+emf9ZEG8tIhKowILC3YuBm4FZwErgRXdfbmZ3m9nF4d1mAbvNbAXwNnCru+8Oot7xJ3Zjxdb9bNh1MIi3FxEJTKAPLnL3mcDMSuvuqvC9Az8OvwLVMjmRbw7vyXMff8YdF5wQdDkiIo2mqQ9mNymXjejNSwu3UFhUEnQpIiKNRkERQdmlsBB6yt0ri3IYfFwqby7bFnBlIiKNR8/MjuBH5wz8wuWvby3bzl/nBjKeLiISCLUoamncCV3I2au7tEUkfigoaikpMYFJI3ode0cRkRihoKiDSaf2BuDA4eKAKxERaXgKijro1r4lAK8trnwjuYhI7NFgdj38/aPPuHxEb8yMqbNX8+CcNeXbfnh2puaBEpGYoBZFPRQcKWbR5n0AmixQRGKWgqIeLh/Rm2c+0vxPIhLbFBT18K2sXvxzxXb2FRwJuhQRkQajoKiHtDYpnH18F15auCXoUkREGoyCop7Kph8PzV8oIhJ7FBT1NLxPR1ISE/hgXSCzn4uINDgFRT2ZGVeM6s0z8zcFXYqISINQUETBV4f24P01u8jdXxh0KSIiUaegiIJ2LZO58OTjeGHB5qBLERGJOgVFLVV+RsXU2asB+PbI3jz3se6pEJHYoyk8aqmqZ1QAnNijPV1SW7I17+juJ03tISLNnVoUUXTFqD5fWKepPUSkuVNQRNFFp3QHYMmWfcEWIiISRQqKKGqRlAjAQ2+vDbgSEZHoUVA0gIWb9rF6R37QZYiIRIWCogFcc0YGD6tVISIxokZBYWbtzWyqmWWHX78zs/b1fXMzG29mn5rZWjObHGG/b5iZm1lWfd+zMVx5Wh/eXb2TTbsPBl2KiEi91bRF8TiwH7gk/NoP/K0+b2xmicBDwPnAIOAyMxtUxX7tgB8C8+vzfo0ptWUyV4zqw1/eXRd0KSIi9VbToOjv7lPcfX349QugXz3fewSwNny8I8DzwMQq9vsl8BugWc2Pcc0ZfZm5dDvb8g5Vub3sxr2yV9mNeyIiTU1Ng+KQmY0uWzCzM4CqPwFrrgdQcc6LLeF15cxsGNDL3WdEOpCZ3VDWLbZz5856lhUdaW1S+Nbwnjzy3voqt+v+ChFpLmoaFDcCD5nZRjPbBPwpvK7BmFkC8ADwk2Pt6+6PuHuWu2d17ty5Icuqleu/3I+X/5PD7gOHgy5FRKTOajSFh7t/ApxiZqnh5f1ReO8coFeF5Z7hdWXaAScC75gZQDdgupld7O7ZUXj/Btc1tSUXntydx+dtCLoUEZE6q1FQmFkL4BtABpAU/uDG3e+ux3svADLNrC+hgJgEXF620d3zgPQKNbwD/L/mEhJlvveV/lz8p/eDLkNEpM5q2vX0GqGB5mLgYIVXnbl7MXAzMAtYCbzo7svN7G4zu7g+x25KeqW15szjuwRdhohIndV09tie7j4+2m/u7jOBmZXW3VXNvmOj/f7RVHGW2IzJM46aJfbW877Ey//JYW1uPgO6tAuyTBGRWqtpUHxgZie5+9IGraYZq276cYDu7VsBcPvLS3nhhtNISLBqj6NpyUWkqYnY9WRmS81sGXA28J/wXdRLwuuXNE6JsaO41Hn+GE/B02WzItLUHKtFMYFQmCwFBjR8ObHtf79+Epf/dT7jTuhCl9SWQZcjIlIjEVsU7r7J3TcA/wd0CS+XvxqnxNhxfLdUJp3ai1+8sSLoUkREaqymVz2NBD40s3XqeqqfH5ydybKcPOas3BF0KSIiNVLTwezzGrSKONIyOZFff+0k/uelJYzq1ynockREjqmmd2armymKzhiQzsh+aTygiQBFpBmoaYtCouzOCwdx7tT3arSvLpkVkSDpCXcBSWuTwr1fPwmAXceYNFCXzIpIkBQUARo3qCsAtzy7iOKS0oCrERGpmoKiCUhMMO7/p8YrRKRpUlA0AQ9OGsLrn2xl1vLtQZciIvIFCooGVvbIU6DaR552atuCP10+lDteXsqGXfWalFdEJOp01VMDizRZYEVDe3fkv88ZyI1PL+SV75/eCJWJiNSMgqIJuWJkbxZt2ssdL9dskl5dNisijUFdT02ImXHP105i7c4DNdpfl82KSGNQUDQxrVISefyqUwF4ZdGWgKsREVFQNEllU5DfM2MlH6zdFXA1IhLvFBRN2B8vG8Ytzy3i0+35QZciInFMQdGEnda/E3ddNIjvPrGA7XmFQZcjInFKVz01cROH9GDL3kNc88SCGv+MroYSkWhSi6IZuGlsf4b06gBQozmhdDWUiESTgiIgNblju4yZ8cuJgwH41YyVjVKfiEiZQLuezGw88CCQCDzq7vdW2v5j4DqgGNgJfDdWHqJU0zu2yyQlhjL9vdU7eXb+Z1w+sndDlSYicpTAWhRmlgg8BJwPDAIuM7NBlXZbBGS5+8nAS8B9jVtl0/PoVVk8MPtTPliny2ZFpHEE2fU0Aljr7uvd/QjwPDCx4g7u/ra7F4QXPwJ6NnKNTU6/zm35/aVD+cFzi9ioCQRFpBEEGRQ9gM0VlreE11XnWuDNqjaY2Q1mlm1m2Tt37oxiiU3T6Mx0fjhuINc9lc3+wqIa/1zZuEjZK9K4iIhImWZxeayZXQFkAV+paru7PwI8ApCVleWNWFpgrhzVhzU78rnl2UU1/pmycZGMyTPKr4oSETmWIFsUOUCvCss9w+uOYmbjgJ8CF7t75IdLx5m7JgyipDQuclFEAhRkUCwAMs2sr5mlAJOA6RV3MLOhwDRCIZEbQI1NWlJiAg99exgAT8zbEHA1IhKrAgsKdy8GbgZmASuBF919uZndbWYXh3f7LdAW+IeZLTaz6dUcLmbU5v4KgPatkgF4+J11/GvFjgavT0TiT6BjFO4+E5hZad1dFb4f1+hFBay291eUmXblcK59MpsnU0dwUs/2DVCZiMSrZjGYLcc2tHdHfv21E7n+qWxevul0juvQqkY/p3mhRORYNIVHDBl/YneuHd2X7z6xgPwaXjareaFE5FjUoogx143py6Y9B7npmf8EXUq11IoRaV4UFDHGzPj5RYO5/qlsAEpLnYQEq/VxovFhXt0xdD+HSPOioIhBSYkJPPzt4Zxw11vc+doy7vnqiZjVLiyq+zCv7sM/UijUNBBqe2wRaRwao2gGanvJLECrlEQAVmzdz69mrMQ9OjfmVTemEY2xjtoeW1OSiDQOtSiagbpeMgvw5DUjmPTXj5g6ezU/PvdLUa4sWFW1WNT6EIk+tShiXPvWyTx97QhmLN3Gn99ZF3Q5DU5XcUlJqVNYVEJ+YRF7Dx4hN7+QrfsOsfvAYUo15U2dqEURB9LbtuCZ60ZxybQPaR3ukoo3amnEriPFpSzctJcfv7iYbXmF5esTDNq1TObg4WKKKwRESmICJ/ZIpWtqS7qmtiStTQrtWyXToXUyqa2Sad8qmY6tU+jZsRXJifpbGhQUcaNb+5Y8c91ILp32YdClBEJXWsWW/YVFvLYoh3c+3cnHG/bQr3Mbvjm8J2MyOzO0d4dqP+ALi0rYmX+YHfsLmfx/S3hz2fbybZ3apDC4R3sWbtrDwcMl5es7t23BLWcPYFD3VI7vnkrbFrH1sbljfyFzVkaeSi+2zlgi6pXWmqevG8nZv3uX1z/ZykWnHBd0SYFTS6N52bDrIE/M28Df53921MzJLZMT+UkNxuBaJifSK601vdJa86+fjI24b8GRYlZtz2f51v38cc4acvM/n7y6R4eWTL10aMRQagr2Fxbx4brdzF+/h/zCIopKSikqcY6UlFJUUsrO/MNs2XuIrwzsHPE4Coo4079zWwB+8foKWiUnMm5Q14ArCpZaGk2fu/P+2l38bd5GPtm8j0kjejHvtrPo1r5lg75v65QkhvXuyLDeHblyVB8g1CJZuGkv76/dxS/fWMHGXQc5tW8aXxnYmXMHd6V7+5pNnVMf7s6+giIOHC4mOTGBpEQjOSGB5CQjwYwV2/Yzd/Uu5q7Zycpt+0lKNPIOFZf//IDObTAz1uQeKF+3Y39hVW9VTkHRjFX8azhj8oxa/TX82FVZfPeJBTyYPJTRmekNWaZInS3evI87X11KUbFzzRkZPPztYbRMDm6crWVyImcMSOeMAencNh6+/vA8/r0ql3+vymXK9OW0aZHILWdlMn5wNzLS29T6+AVHitl94Ah7DoZeuw8eYWf+YXL2FZCz9xA5+w6xNvcAFcfkkxKMVimJHDxcfNT647u1Y+HPzqnxf68Xb6x+m4KiGavPZbOn9OrAn68Yzo1/X8i0K4dHubLmr6ouKUDdVI1kf2ER98/6lOc+/oyiktCn3+SXl/LKohxe+K/TAq7ucy/fdEb590Ulpcxfv4e3lm/jW9M+pFObFAZ1TyUlKYGUpARahL+mJCZy4HARuw+EguDzUDiMe2isJK1tCpv3HCLv0Odztg3s0pYHLxtKj46tSG2Z3KjnqaCIYyP6pvHgpCHc+PTCoEtpcqrrklI3VcNyd2Yu3c7dbyznrOO7sOCn4+jQOiXosmokOTGB0ZnpjM5MZ/X2fD7euJdV2/MB6JPWmstH9ubpDzexZd+h8p8Z2LUt933zFKa8toyc8PqteYX0SmvNJ1PODeQ8qqKgiHNjMjtz7zdO5vqnsvnjnDVcN6Zf+V3dIo1p854CfvbaMrbuO8RDlw8jKyMt6JLq7MUbT69y/X99pX+V61+7eXRDllNvCgrhnPCA9qod+Zx5/zv85NyBfH1YTxLrMJlgPNKVU/Xj7px1/zts2F1Qvu63sz5tUl1M8U5BIeUeunwYCzft5dczV/L4vI389IITgi6pWdCVU3VXdj9Dm5ZJ/PNHX2Zg13ZBlyRVUFDIUYb36chLN57GW8u2c+erSwH42avLOLFHKoOPa8/Aru1ISWq61403Z/HUMnF3Xl+yjbtfX863R/bh5rMGNOn7EeKdgiIG1eeyWQg90+L8k7pz9gldGXjnm2Skt2H++j089v4GPttTwIAuoXsxHn5nLf3S2zKgSxt6p7X5QoCUljqHikJ3uLp7rac6jzfx0jLZlneIX72xkk935PP41adycs8OQZckx6CgiEH1uWy2orIP/mtH9y1fd+hICSu37+frD3/AvoIi/pG9mfW7DpKz7xDdwzdAnXrPvyg4XExBUQktk0ID4yfc9RY9O7amZ8dW9OrYml5poRuTikpKY/YvyUjP15j97FwAVnTNiemWQ0Ubdh3kkr98yM4Dn9/hfM+MlRqLaAYUFFIrrVISGda7IwB3VBjDOFxcwuY9hxj3wLu8fvNo2rRIpHVKEokJRsbkGWTfeQ5b9hawec8hNu8pYPPe0MDlqF/P4cKTuzNxSA+G9e4QU62OSJfY/qhzARc8ODemWw5lVmzdz8PvrOWDdbu5YlQfrjk9g45tmsclrxKioJCoaJGUWN4lVdXUCm1bJHF8t1SO75Zavu5v8zby8k2n89rirdz60icUlzgTh4Tmn1JXVfO3M/8wt7+8hCVb8rhuTF/u/cbJMTehXrwI9F/NzMYDDwKJwKPufm+l7S2Ap4DhwG7gUnff2Nh1SsPp06kNPzg7k1vOGsCynP28ujgHgDH3vc24E7oy7oSujOibpgH0ZubiP77Pkpy88uU5K3O54ctV30MgTV9gQWFmicBDwDnAFmCBmU139xUVdrsW2OvuA8xsEvAb4NLGr1YamplxUs/2nNSzPY+9v4FHr8riXyt2cP8/P2X9zgOMCc9uuXlPAb3SWgdcrUTy8YY9bM0r5LffPJlvZfUKuhyJgiBbFCOAte6+HsDMngcmAhWDYiLw8/D3LwF/MjPzaD0AOs7U92qoxlTWTXXzWZnk5hfy9qpcZizZxtce/oAWSQmc1r8Tp/fvxGn9OwVdqlTw+idb+fn05Tw4SZNNxhR3D+QFfJNQd1PZ8pXAnyrtswzoWWF5HZBexbFuALKB7Pbt2ztQ/srOzvbs7Oyj1k2ZMsXd3bt3716+btiwYe7ufv311x+1b05Ojk+fPv2oddOmTfNwWJW/JkyY4O7uEyZMOGq9u/u0adOOWjd9+nTPyck5at3111/v7u7Dhg0rX9e9e3d3d58yZUqDn1PaeTdXeU59bnujynNKO+/mGp9TStf+UTunocOG+Zod+/28m3/t6V+93Xve8owfd/00/8HTH/qUx6Z7QusOEf+dWvU/tcp/pz63vVHlv1OPm578wjn1ue2NKv+d2p9xWc3PKT09pn73SktLvce4q73H9x735M4ZMXFOcfgZkV35s7XsFRNBUfE1fPhwl9rrc9sbDba+IY9dUlLqfW57wx+du96vfWKBnzTlLR/3u3e8z21v+IfrdnlxSWnU3zMq57NokZ9/9R9qdZymamd+od/6j8V+3tR3feu+gqDLkTqKFBRBdj3lABU7MHuG11W1zxYzSwLaExrUFgEgITwf1bWj+3Lt6L6UlDrLt+Zx8Z/mcffrK8jNP8z4E7tywYndGdG3+U4y1xTlFRTxyNx1THt3ffkzqU/7338zsm+a7o2IMUEGxQIg08z6EgqEScDllfaZDlwFfEioBfLvcPKJVCkxwcrv9J35wzFs3HWQN5dt5963VpGzNzSN89IteZzYI1WX39bRgcPF/O39DTw+bwPnDurGO7eOpWdHXWAQywILCncvNrObgVmELo993N2Xm9ndhJpA04HHgKfNbC2wh1CYiNRYRnobvje2P98b25/NewoYc9/bfO+ZhaS2TGbSiF5MPKUH7Vs37kNgmiN3Z/nW/fxzxQ6enb+JMwak8/JNZ9C3Dk9xk+Yn0Pso3H0mMLPSursqfF8IfKux64onzelKqPoqu6z2vVvP5IN1u3l+wWf8dtannH18l4Ara5oKi0r4cP1u5qzcwZyVuew5eITDxaUAvLZ4K9vzCtXFFCd0m2Sci9a8UM1JQoKVP4ls78EjvLwoh1cXb+VrD8/j+jH9OG9wt7h4FkdJqfPmsm1MX7yVQ0UlHC4u5UhxKUUloa/b8wr5Urd2jBvUlaevHUn/zm3UXRendLurxLWObVLKJz28YUw/Hp27nrH3v83f5m3g4OHigKtrGEUlpfwjezODp7zFzc8u4p8rdjB3zS7yDxXxswmDMGBN7gHyDxeTvWkvb6/KZUCXtgqJOKYWhVQpnrqkypx/UnfOP6k7Czft5dG56/n9v0Lnf9kjH9GhdTLtWyWXj2fsLyxq9Afc11dhUQkvZm9m2rvr6dOpNY9fdSqn9e/0hQB44wdjAqpQmioFhVQpHrukygzv05HhfYazM/8wp97zL75/5gDyDhWx79AR8g4VATD2t+9wzekZXH1GBu0aOTAOF5ewYut+Fn22j50HDtMiKYEWSYmhr8kJpCQmcOBwMbn5h9mxv5Cd+YfJ3X+YrfsOMaJvGn+8fGj5DMAiNaGgkFqJp5ZG53YtAL4wFcV9b33KP248jT/OWcPY377Dd0f35arTMxqkhvzCIjbvOcTanQdY/Nk+Fm3ey6pt+ZhBwZGS8v26t2/JeYO7MXPpNnLzP3/eQ9/0Ntx10SB+N+tT8g8XM2dVLnNW5epeB6kVBYXUSjy3NCrq37ktv580lLW5B/jDnDV85b63ASguKSWpjg9i2rK3gBcXbAbg4j+9z+Y9BeQdKqK0wp1DJ3RrR/ad42hTzXTdP794cJXrz/ySruySulNQSFRU1dIAYr71MaBLW/5w2VDW7MjnnKnv8bWHP+A33ziZQcelHvuHw9bm5vPnd9YzZ9UOvjGsJxD6wO/VsTXpbVM0iCyBU1BIVFTX0oi1YKhOZtd2AFw5qg9XPjafSSN6cctZmbRMTqz2Zw4lpXDj0wvJ3rSHq0/P4N1bz6R9q2Qee3+DxhCkSVFQSKOL5XGOS07txdgvdeau15ZzwR/m8ptvnAyEnva2JjeftbkHWL0jn+VrctncoRtX9k3jgUtPoXWK/leUpku/ndLoYn2co0tqS/5y5XDeWraNm5/9DwDnTH2XgV3aMaBrWzK7tGV8y1TueeFTvjt6YsDVihybgkKkgYw/sTtjMjszeMosFv3snKPHGhbnoZEHaS50Z7Y0GVNnryZj8gwg1CU1dfbqgCuqv7KrkzQgLc2ZWhTSZMR6l5RIc6UWhTR5sdjSEGlO1KKQJk8tDZFgqUUhIiIRKSik2VKXlEjjUNeTNFvqkhJpHGpRSExRK0Mk+tSikJiiVoZI9KlFIXFBLQ2RulOLQuKCWhoidacWhcQ1tTREji2QFoWZpQEvABnARuASd99baZ8hwJ+BVKAEuMfdX2jUQiXmqaUhcmxBtSgmA3PcPROYE16urAD4jrsPBsYDvzezDo1XooiIQHBBMRF4Mvz9k8BXK+/g7qvdfU34+61ALtC5sQqU+FZVl5S6qSReBTWY3dXdt4W/3w50jbSzmY0AUoB11Wy/AbgBoHfv3lEsU+JVEI92jeUn/0nz1mBBYWb/ArpVsemnFRfc3c3MIxynO/A0cJW7l1a1j7s/AjwCkJWVVe2xRJoyjZdIU9VgQeHu46rbZmY7zKy7u28LB0FuNfulAjOAn7r7Rw1Uqki9qCUgsS6orqfpwFXAveGvr1XewcxSgFeAp9z9pcYtT6Tm1BKQWBfUYPa9wDlmtgYYF17GzLLM7NHwPpcAXwauNrPF4deQQKoViaKps1dzwYNzAQ2KS/Ng7rHVpZ+VleXZ2dlBlyFyVJcUcHSX1OLFoa9DhjR6XSJVMbOF7p5V1TZN4SHSQNQlJbFCU3iIiEhECgoREYlIQSEiIhEpKEREJCIFhYiIRKSgEBGRiBQUIiISkYJCREQiirk7s81sJ7ApvJgO7AqwnMai84wt8XCe8XCO0LzOs4+7V/nMn5gLiorMLLu6W9Jjic4ztsTDecbDOULsnKe6nkREJCIFhYiIRBTrQfFI0AU0Ep1nbImH84yHc4QYOc+YHqMQEZH6i/UWhYiI1JOCQkREIorZoDCz8Wb2qZmtNbPJQdcTLWb2uJnlmtmyCuvSzGy2ma0Jf+0YZI31ZWa9zOxtM1thZsvN7Ifh9bF2ni3N7GMz+yR8nr8Ir+9rZvPDv7svhJ8f36yZWaKZLTKzN8LLsXiOG81safixzdnhdTHxOxuTQWFmicBDwPnAIOAyMxsUbFVR8wQwvtK6ycAcd88E5oSXm7Ni4CfuPggYBXw//O8Xa+d5GDjL3U8BhgDjzWwU8BtgqrsPAPYC1wZXYtT8EFhZYTkWzxHgTHcfUuHeiZj4nY3JoABGAGvdfb27HwGeByYGXFNUuPt7wJ5KqycCT4a/fxL4amPWFG3uvs3d/xP+Pp/QB0wPYu883d0PhBeTwy8HzgJeCq9v9udpZj2BC4FHw8tGjJ1jBDHxOxurQdED2FxheUt4Xazq6u7bwt9vB7oGWUw0mVkGMBSYTwyeZ7hLZjGQC8wG1gH73L04vEss/O7+HvgfoDS83InYO0cIhfw/zWyhmd0QXhcTv7NJQRcg0eXubmYxcc2zmbUF/g/4b3ffH/pDNCRWztPdS4AhZtYBeAU4PtiKosvMJgC57r7QzMYGXE5DG+3uOWbWBZhtZqsqbmzOv7Ox2qLIAXpVWO4ZXherdphZd4Dw19yA66k3M0smFBLPuPvL4dUxd55l3H0f8DZwGtDBzMr+iGvuv7tnABeb2UZCXcBnAQ8SW+cIgLvnhL/mEgr9EcTI72ysBsUCIDN8ZUUKMAmYHnBNDWk6cFX4+6uA1wKspd7CfdiPASvd/YEKm2LtPDuHWxKYWSvgHELjMW8D3wzv1qzP091vd/ee7p5B6P/Df7v7t4mhcwQwszZm1q7se+BcYBkx8jsbs3dmm9kFhPpGE4HH3f2eYCuKDjN7DhhLaPriHcAU4FXgRaA3oSnWL3H3ygPezYaZjQbmAkv5vF/7DkLjFLF0nicTGuBMJPRH24vufreZ9SP013casAi4wt0PB1dpdIS7nv6fu0+ItXMMn88r4cUk4Fl3v8fMOhEDv7MxGxQiIhIdsdr1JCIiUaKgEBGRiBQUIiISkYJCREQiUlCIiEhECgqRBmJmU83svysszzKzRyss/87MfhxIcSK1oKAQaTjzgNMBzCyB0L0vgytsPx34IIC6RGpFQSHScD4gNCUHhAJiGZBvZh3NrAVwAvCfoIoTqSlNCijSQNx9q5kVm1lvQq2HDwnNknoakAcsDU+DL9KkKShEGtYHhELidOABQkFxOqGgmBdgXSI1pq4nkYZVNk5xEqGup48ItSg0PiHNhoJCpGF9AEwA9rh7SXhCuA6EwkJBIc2CgkKkYS0ldLXTR5XW5bn7rmBKEqkdzR4rIiIRqUUhIiIRKShERCQiBYWIiESkoBARkYgUFCIiEpGCQkREIlJQiIhIRP8ffHvdqMV30JoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAe0ElEQVR4nO3dd3wc9bnv8c+jYkmW5CrJ4CrHhQROYpMrbMAQMD2EhJwUgg+54SYEkuCcEEg5Jg3ICYTk0NJugsMhcHMoIYUWY4IxNXTb2GCKbYqNLWxL7upWee4fMxZrWStL8hbtzvf9eu1rd347nnkGVt/57W9mZ8zdERGR6MhJdwEiIpJaCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYkYBb9EhpmdY2YPpbsOkXRT8MuAZGb1MY8OM2uKmT6nP8t099vc/ZRerv//mNk/+7OemGWcZ2avm1mdmW02swfMrLSPy/iOma0Ml/G2mX1nP/OfGK6z0cweNbMJB7INkp0U/DIguXvJngfwDvDxmLbb0l3f/pjZccBVwBx3LwU+APypP4sCvgAMB04Dvm5mZ8dZZxnwN+CHwAhgST/XKVlOwS8ZxcxmmNkzZrbDzDaa2a/NbFD4XqWZuZnlxcz/mJl9OXy9Vy8+nPerZrYmXN5vLPAB4HfAUeE3jB39KPUI4Bl3fxHA3be5+63uXteXhbj7z919mbu3ufsq4F5gVpzZPwW84u5/dvdm4HJgmpm9vx/1SxZT8EumaQcuBsqAo4ATgQsPYHlnEIT0h4CzgFPd/TXgqwTBXeLuw/qx3OeAU83sCjObZWYFsW+a2bxwZ9Pto7sFmpkBxwKvxFnnYcCKPRPu3gC8GbaLdFLwS0Zx96Xu/mzYA14L3AgcdwCLvNrdd7j7O8CjwPQElIm7P0nQA/8wsADYambXmVlu+P7V7j4s3iPOYi8n+Jv9Q5z3S4CdXdp2An06riDZL2//s4gMHGY2FbgOqAIGE3yGlx7AIjfFvG4kCM+EcPeFwEIzywFmA38GVhHsrPrEzL5OMNZ/rLu3xJmtHhjSpW0I0KfhJcl+6vFLpvkt8Dowxd2HAN8jOAAK0BA+D46Z/6B+ridhl6119w53Xww8AvwLgJl9r8uZS3s9Yv+9mX0JmAec6O4beljVK8C0mH9XDEwi/tCQRJSCXzJNKbALqA8PWn5tzxvuXgtUA583s9wwMCf1cz2bgbF7Dhx3x8wuN7PH4rx3ppmdbWbDwwPGMwiGpJ4Na70q9sylro+Y5ZxDcHbQye7+1n5qvhv4FzP7tJkVAj8CXnL31/uy4ZL9FPySab4N/BvB8MXv2fd0xfOB7wBbCQ5qPt3P9TxC0FPeZGZb4swzDngqznvbw1rWEOyo/gf4r36civoTYCTwQsw3gt/tedPMXtnzu4Zwx/dp4Mpw/TOBbk/9lGgz3YhFpH/MbDnB8MvWdNci0hcKfhGRiNFQj4hIxCj4RUQiRsEvIhIxGfEDrrKyMq+srEx3GSLxNTYGz4MH9zyfSAotXbp0i7uXd23PiOCvrKxkyZIl6S5DJL7ly4Pn6dPTWYXIXsxsXXftGuoREYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRiFPwiIhGj4BcRiRgFv4hIxCj4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRiFPwiIhGTtOA3s5vNrMbMVsa0XW5m1Wa2PHycnqz1i4hI95LZ478FOK2b9uvdfXr4eCCJ6xcRkW4kLfjd/QlgW7KWLyIi/ZOOMf6vm9lL4VDQ8HgzmdkFZrbEzJbU1tamsj4RkayW6uD/LTAJmA5sBK6NN6O7z3f3KnevKi8vT1F5IiLZL6XB7+6b3b3d3TuA3wMzUrl+ERFJcfCb2cExk/8KrIw3r4iIJEdeshZsZncAxwNlZrYBuAw43symAw6sBb6SrPWLiEj3khb87j6nm+b/Ttb6RESkd/TLXRGRiFHwi4hEjIJfRCRiFPwiIhGj4BcRiRgFv4hIxCj4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRiFPwiIhGj4BcRiRgFv4hIxCj4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYmYpAW/md1sZjVmtjKmbYSZLTKzNeHz8GStX0REupfMHv8twGld2uYBi919CrA4nBYRkRRKWvC7+xPAti7NZwK3hq9vBT6ZrPWLiEj3Uj3GP8rdN4avNwGj4s1oZheY2RIzW1JbW5ua6kREIiBtB3fd3QHv4f357l7l7lXl5eUprExEJLulOvg3m9nBAOFzTYrXLyISeakO/vuAc8PX5wL3pnj9IiKRl8zTOe8AngEOMbMNZnYecDVwspmtAU4Kp0VEJIXykrVgd58T560Tk7VOERHZv/32+M1sVm/aREQkM/RmqOdXvWwTEZEMEHeox8yOAo4Gys3skpi3hgC5yS5MRER67/pFq/nF4jWd0xedOCXuvD31+AcBJQQ7h9KYxy7gM4koVEREEuPik6ey9uqPAbD26o9x8clT484bt8fv7o8Dj5vZLe6+LuFViohIn3XXs+8p5LvTm7N6CsxsPlAZO7+7n9CnNYmISK/FC/g9j8p5Czp7+H3Vm+D/M/A74CagvV9rERGRuHoK+QMJ+Hh6E/xt7v7bhK5VREQ6JTPku9Ob4L/fzC4E7gZa9jS6e9dLLouISA8SMT6fCL0J/j3X1vlOTJsD70t8OSIimS+Z4/OJsN/gd/eJqShERCTTDPSAj6enH3Cd4O6PmNmnunvf3f+WvLJERAa+gR7w8fTU4z8OeAT4eDfvOaDgF5FIGChj84nS0w+4Lgufv5i6ckRE0ivVp1amw37H+M3sR921u/uPE1+OiEhqZOr4fCL05qyehpjXhcAZwGvJKUdEJDWiEPDx9Oasnmtjp83sGuAfSatIRCSBsm18PhH6cweuwcDYRBciInIgojx001e9GeN/meAsHgiuw18OaHxfRAYUBXzv9abHf0bM6zZgs7u3JakeEZH90vDNgenNGP86ADOrIDi4O9rMcPd3kl2ciESbhm+SozdDPZ8ArgVGAzXABIKzeg5LbmkigXh//PFuNaeeYPZQwCdHb4Z6/hM4EnjY3Q83s9nA55NblmSzvgR5T727ntq7tvV1nRpKSC39906t3gR/q7tvNbMcM8tx90fN7IZkFyaZpT9h3tsgT4T+7Dy6tiucDpyGbgaG3gT/DjMrAZ4AbjOzGqA+uWXJQJWoXnkmirct1y9azaLbnwTg1VHV2iH0IJs+D5msN8G/AmgELgbOAYYCJcksStJPPbPeu/jkqVxc3sjpv3hSw0sxsm17sklvgn+2u3cAHcCtAGb2UlKrkpRRwCdPMoeX0hGqiTo2I+nX0/X4vwZcCEzqEvSlwFPJLkwSLwpXHcw2fd15JOIMqCgP52WqnY2tvFFbxxs1wSj8F//wPG/Uxh+R76nHfzuwEPgpMC+mvU732x3Y1AOLrkScAaXPycDU0tbOhu1NrNvawNotjQCcPf8Z3qhpoLm1nUnlxUyqCEbh58wYz+SKEib9R/fL6ul6/DuBncCcRG+AJIYCXiT7rNpUx9qtDUHAbw0CftbVj1Bb18LoYYVMGFlM5cjBAMydPZnJFSUcNKQQMwPgb8uqOeWwg3pcR38u0iYppoAXyR51za0ALHhp414Bv25rcAX8ubcvY8KIwUwYWcz7DyoF4PbzZzJ6WBH5uTmdy7n1mXUcO6W8XzUo+DOAAl4kM7g7u5raWL+9keodTQBccf8rVG9vYsP2JjZsb6S1Pbjm5X0rqqkcWcz0ccP55PQxTCgrZtbVj/DwJcfttcwf3fsKE0YWJ7ROBb+ISB9sa9jNhu2NbNjeFAZ6MBxz2g1PsGF7EwaMGV7E2OFFAIweWsSMyhFh22CGD85n4qUPcOP/rkrbNij4Bxid+yySXjsbW3k7HHb53eNvsmF7Y2ePHWD2NY8xNgz2McMGd/bGrz1rGmOHD2ZoUX7nsirnLeD8j7wv9RuxHwr+AWYgDOu4Oy1tHTS0BFffXlm9k4aWNupjHgDXPbSKupa2zvfqmoP2T//2adraO2htd9o6OmgLv9oeedVi2jrC9vYOAA790YPkmmEGuTlGTniAauZVD4ftFrYHtZ35639SkJ9LYX4uRfk5FOXnAnD1wtcpLcxjSGEepYXBH96K9TsYWTKIspICCsP5RAAad7exdksja8OA/9ZdK1i7tYG3tzSwu62DyrLg4Om2ht1MqShl9iEVjBlexGk3PMmKy07ZZ3k//vurHDZ6aEq34UAo+NMk2T37jo4gbNduaWBnUys7mlrZGT4ArlzwKjsa92475mePdIa4YZQUBh+P7/zlJUoL8iguyKWkMJ+SgiBEc3KMMcOKKC3Mo7ggeDy5ZgvzPvp+8nKM/Nwc8nKNvJwcTrrucf524dHk5Rr5OUH7By9/iOe/fxId7nR0OB0O7R3OEVc+zD1zZ9HhhO1Oe4dzwrWPc/knDqO5tYPmtnaad7fT1NrOPcvfpbQwj13NrVTvaKI+3AH94J6VbK1vYUvDbvJzjJElBQB8+dYllJUMYmTJIEYWFzCyZBAA67c1MmpIIYPy3juAJpnL3amtb2H1pnpWb64DYM78Z3l7SwPbG3czYeRgJpYFvfUZE4fzuSPGUVk2mPKSAsyMynkL+N7pH0jnJiRNWoLfzNYCdUA70Obu6RvsSpMD7dnvbutg3dYG1oQ/2Lj8vlfYtLOZTbua2byrmS31LQB84ebnGVqUz7DB+Qwpyu/8GlpWUsDkihKGFuUztGgQc37/LHecf2QY4LkU5AXhXjlvAQsvOnaf9d/x/Hq+eVL3O6ojKkd02z56WNE+bSUF3X8EDx6677wAh48fvk/bJXetYO7syXu13bfiXe7/92OAIADqW9rYWr+b4695jM9WjWVr/W621rfwzrZGlr6zHYCz5z9LTV0zI4oHMXpYUWe9f1m6gYllxbyvrJjhxYO6rUvSa3vDblZtrmNNGPBn3fgMqzfXYcDUUaVMHRWcHXPh7ElMLCvm4KFF5IZfIyvnLeBzR4xPV+lpkc4e/2x335LG9WeMbQ27eTEMp6/9z1LW1NTzzrZGxgwrYnL4g41xIwZzROUIDhpawKghhVSUFjL1Bwt54ruz91ne7c+9w1eOm7RP+7gRg5O7IWliZpQW5ncOAZ3azTnOC15awFPzTqCtvYOauhbe3dFE9Y4mFry0kX+uqeWPz6zlrdoGcnODsLjkruVMKi9hYlkxE8uKqWzrSOk2RdWu5lbWbK5j9eZ6Vm0KQr7qJw/T0trO1INKmToq+Hu46MQpTB1VSlnJoM7z2//4bP9Pf8w2GupJsr4O6bh7Zy/+239ewbJ126mpa2H6uGEAfPSDB/ONiiBw9oxbV85bwHnHTEzeRkRIXm5OZ2+/CrjozuXccPbhQPD/ZmvDbqp+8jAzJ47grS0N3PNiNW9vaeCdrfW0l43nX//vU4weWsTBQws5eFgRo4cWAlCzq5mykgJy9hyskB7Vt7SxZnMda2rqO3vxR/90MTuaWplSUcKUUaUcEvbi7//3WXv9gOmO59cza3JZ2mrPBOkKfgceMjMHbnT3+WmqI+l6M6TT0NLG029u5ZHXa3hsVU3nV9DDxw/jvGMmMnVUKbk5wZjjJ6aNTmX5EsPMKAuPE3QdGmhf9iKn/uZZvjf3WDbubGbjjibWb2vkube2AnD6L59kV1MbFUMKgh3DsGCHcNOTb1FeWkBFaWHwPKQgtRuVZnXNrbxRU8+azfWsqXkv4Lc3tjKpopgpFaVMCXvxf/rKUYwZVrTXzvPKB16LOywo8aUr+I9x9+rwPr6LzOx1d38idgYzuwC4AGD8+Owbf9ta38KClzcCMPOqxUwbN5TZh1Rw3jEzmVRezMRLH+CcmRPSXKX0Vm6OMaijrdvjG5XzFrDkByfT3NrO5l3NvLujmY07m7h3+btU72hi+fod1NS1UBs+AI79+SNUlBZSUVpAeWmwM7jrhfWUh9MVQwoYWZw5O4lVm+p4Z1tj8AjPpIkN+KkVpUwOA/7OC45izPD3xuABfv7gqqwdikyHtAS/u1eHzzVmdjcwg+BGL7HzzAfmA1RVVXnKi+yj3gzpNO1u56FXN3HPi9UsWbedE99fAcAzl57QOf4s2aswP5cJI4s7z/u+5K4VXPbxfW9dXTlvAX/80szOnUFNXTMAz6/dRk1dCzXhwfsdjcHZWCdc8xhDit47eD+0KPizvvHxN8Pp4DEkPLBfs6uZwQV5DM7P7ffQU21dC3XNrdQ1t4WPoJZrH1pFza6g5pq6FjbvCnZkc29fxvgRg4NHuP13XnAUY4fv3YP/+YOrGD9SAZ9sKQ9+MysGcty9Lnx9CvDjVNeRaPGGdNyd594OLmY646qHOXz8cD51+Bh+/W8fprggLzwVUaEve6ssK6ay7L2f6V9x/6tc89lpe83T2t7BlO8vZP4XqtjZ1Mqu5lZ2xZyeu6W+hTdr64P3mto62z/2q3/S2NJGY2s7heHZWzOvepi8nBzyc4PfTey5Jszsax6jpbWdlrYOmsNngFNveCL83UQ+pYV5lIan/ubmGNPGDaMi/FYyakghM69avM9lCP7z768q4NMoHT3+UcDd4YGYPOB2d38wDXUkVePuNu5+sZr/9/Q62jqCP5bF3zqOitLCNFcm2WJPOO85syvWj+59he9/7NB92ivnLeCF758EBL+RaGpt57DL/sE9c2fR1u60dXjnj+9O/+WT3HRuFQV5ORTm51KQl0NBXi5Tf7CQZT88udtlxzvFVwaWlAe/u78FTNvvjBnsJ39/lb8s20DVhBH88IxDmTV5JBMvfUChLwNKTo5RHP6OIt4B0knlustqNtJPFPvo+kWrqZy3oPNx/aLVQHDwau5ty4Dg6+79Xz+Gm86t4pgpZZ2nmSWrFmCvWvraLiLRovP4+6jrWP7rm3Zx4W1Lef7t7Zx/7EQWvLyRSw/gZ96xB4kr5y3Y5xZ5e9pir8cfr8betHe3vq633zuQdhEZeBT8B+Crf1zKknXbueAjE7nms9MYPCiPny58vVf/Nl5Q9hTayZCInUdP7X3ZkWnnIZIaCv4+qq1r4b/+EYR7VeVwrv/cdIoGxb/yY18DPtv0ZUfWl52HdhQi/afg70HXc/M/MqWMle/u4lOHjwHgy8fu/zrbUQn4ZErE0JWIvEfB34PY8fwpFSU4cNdXjmRyRSk3/fPtveZV4AwMff3WIBJFCv4ebGvYzQ/vXQnAt089hFMOHRX3DB317Ac2DSOJvEfBH8cTq2v57l9e4uPTDgbeu5SvAiG7aBhJokjB30Vzazs/e/B1Hly5ievOmsbRk8v4/ZPvDeuoZx9dPX1rWHT7k4B2CJIZFPzsexB3SkUJCy86lmGDdbcl2b+LT57KxeWNwcT06Z3t+oYgydCX4cl4FPwEf7hlJYP44b2vcM1np/HpD4/BzPSHKwdExxWyXzJ/AAkH/sPNS+LUHfngb2vv4Md/f5Vn3gxumPGZ/zW28z0N60gy6LhCag2UEO5vezJEOvh3Nbcy97ZlmBl/vfBoPnT5Q+kuSWQvUf/WMNDCOVuY+4C/xwlVVVW+ZMmShC7zna2NnHfrCxw1aSRDC/P51aNvdL6X6X8skgbLlwfPMWP8A0l3NwoC9mnrGqqx7XvEu41o1/ZErlP6x8yWunvVPu1RCv6uH67jp5Zzy5dmHPByRQZ68CdCvHBWaA9cCv7QAy9v5MLblnHLF4/g+EMqErJMkSgEv2SeeMEfqTH++1e8yxX3vwqg0BeRyIrMjVjuXV7Nf/z1JbbUBzd/1o1IRCSqItHj/+vSDfzswde5Z+4spo4qTXc5IiJplfXBf9eS9Vz30GpuP38mkysU+iIiWR38dz7/Dr9YvIbbzp+pm0aLiISyMvi7nl523/J3dXqZiEgoKw/uVm9vSncJIiIDVtb1+N+sreex1bX88bwZHDulPN3liIgMOFnV499S38IX//AC3z31EIW+iEgcWRP8Tbvb+fKtSzhz+mjOOmJcussRERmwsiL42zucb/7pRSpHDuYSHcQVEelRVozx//SB19jZ1Mov5xwe92boIiISyPjgP+m6x3mjph6AQ37wIDMnjuBPXzkqzVWJiAxcGR3867c1sr1hN/fOncW0ccPSXY6ISEbI2DH+1vYOvnHni3zt+EkKfRGRPsjY4L/2odUMK8rnS7MmprsUEZGMkpFDPY+vruWeF6tZ8I1jyMnRwVwRkb7IuOCv2dXMt/+8gl+efTgjSwrSXY6ISMbJqKGejg7n4ruWM2fGeI6aNDLd5YiIZKSMCv7fPv4mrW3ON06YnO5SREQyVsYM9Xzuxmd47u1tAEz+/kKdry8i0k9pCX4zOw34BZAL3OTuV+/v3yjkRUQSI+VDPWaWC/wG+ChwKDDHzA5NdR0iIlGVjjH+GcAb7v6Wu+8G7gTOTEMdIiKRlI7gHwOsj5neELaJiEgKDNizeszsAjNbYmZLamtr012OiEjWSEfwVwOxd0oZG7btxd3nu3uVu1eVl+tuWiIiiZKO4H8BmGJmE81sEHA2cF8a6hARiaSUn87p7m1m9nXgHwSnc97s7q+kug4RkahKy3n87v4A8EA61i0iEnUD9uCuiIgkh4JfRCRiFPwiIhGj4BcRiRgFv4hIxCj4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRiFPwiIhGj4BcRiRgFv4hIxCj4RUQixtw93TXsl5nVAuvCyTJgSxrLSZUobGcUthG0ndkk07ZxgruXd23MiOCPZWZL3L0q3XUkWxS2MwrbCNrObJIt26ihHhGRiFHwi4hETCYG//x0F5AiUdjOKGwjaDuzSVZsY8aN8YuIyIHJxB6/iIgcAAW/iEjEZEzwm9lpZrbKzN4ws3npridRzOxmM6sxs5UxbSPMbJGZrQmfh6ezxkQws3Fm9qiZvWpmr5jZRWF71myrmRWa2fNmtiLcxivC9olm9lz42f2TmQ1Kd62JYGa5Zvaimf09nM667TSztWb2spktN7MlYVvGf2YzIvjNLBf4DfBR4FBgjpkdmt6qEuYW4LQubfOAxe4+BVgcTme6NuBb7n4ocCQwN/x/mE3b2gKc4O7TgOnAaWZ2JPAz4Hp3nwxsB85LX4kJdRHwWsx0tm7nbHefHnP+fsZ/ZjMi+IEZwBvu/pa77wbuBM5Mc00J4e5PANu6NJ8J3Bq+vhX4ZCprSgZ33+juy8LXdQSBMYYs2lYP1IeT+eHDgROAv4TtGb2Ne5jZWOBjwE3htJGF2xlHxn9mMyX4xwDrY6Y3hG3ZapS7bwxfbwJGpbOYRDOzSuBw4DmybFvD4Y/lQA2wCHgT2OHubeEs2fLZvQH4LtARTo8kO7fTgYfMbKmZXRC2ZfxnNi/dBUjP3N3NLGvOuTWzEuCvwDfdfVfQUQxkw7a6ezsw3cyGAXcD709vRYlnZmcANe6+1MyOT3M5yXaMu1ebWQWwyMxej30zUz+zmdLjrwbGxUyPDduy1WYzOxggfK5Jcz0JYWb5BKF/m7v/LWzOym119x3Ao8BRwDAz29PJyobP7izgE2a2lmDY9QTgF2TfduLu1eFzDcGOfAZZ8JnNlOB/AZgSnjUwCDgbuC/NNSXTfcC54etzgXvTWEtChGPA/w285u7XxbyVNdtqZuVhTx8zKwJOJjiW8SjwmXC2jN5GAHe/1N3Hunslwd/iI+5+Dlm2nWZWbGale14DpwAryYLPbMb8ctfMTicYV8wFbnb3K9NbUWKY2R3A8QSXe90MXAbcA9wFjCe4HPVZ7t71AHBGMbNjgCeBl3lvXPh7BOP8WbGtZvYhgoN9uQSdqrvc/cdm9j6CnvEI4EXg8+7ekr5KEycc6vm2u5+RbdsZbs/d4WQecLu7X2lmI8nwz2zGBL+IiCRGpgz1iIhIgij4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX6QPzOx6M/tmzPQ/zOymmOlrzeyStBQn0ksKfpG+eQo4GsDMcgh+f3FYzPtHA0+noS6RXlPwi/TN0wSXYYAg8FcCdWY23MwKgA8Ay9JVnEhv6CJtIn3g7u+aWZuZjSfo3T9DcBXKo4CdwMvhpcNFBiwFv0jfPU0Q+kcD1xEE/9EEwf9UGusS6RUN9Yj03Z5x/g8SDPU8S9Dj1/i+ZAQFv0jfPQ2cAWxz9/bwAl3DCMJfwS8DnoJfpO9eJjib59kubTvdfUt6ShLpPV2dU0QkYtTjFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRi/j/rsirPrnndsgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c_obs3.plot_rho()\n",
    "c_obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now redo the error analysis and alter the value of S or attach a tail to the autocorrelation function to take into account long range autocorrelations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 3.27194697e-01 +/- 2.14280114e+00 +/- 2.48970994e-01 (654.901%)\n",
      "  t_int\t 6.41297945e+00 +/- 2.18167829e+00 tau_exp = 20.00,  N_sigma = 1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlR0lEQVR4nO3deXxcdb3/8ddnsk32bmkptDRQiggCFSMCAlIWBQG3y0NFvD+8LogogijXCirlCoILKLigrNWrqFy1gFSWUpZCa4ECbSlUaAstpWu6ZE8myczn98ecDNM2aZM0k0nmvJ+Px3nMnO+ZnPM5yeQ93/meM2fM3RERkfCIZLsAEREZXAp+EZGQUfCLiISMgl9EJGQU/CIiIaPgFxEJGQW/DHlmdp6ZPZLtOkRyhYJfBpSZNaVNCTNrTZs/rz/rdPc/uvsHe7n9z5nZ0/3ZTvDz1WbmZpbf33UMB2Z2ppk9bWZ1ZrbRzG43s/K05UVmdqeZNQTLL8tmvTKwFPwyoNy9rGsC3gTOTmv7Y7brk5RK4BpgX+CdwH7AT9KWzwCmAJOAacB/m9npg1yjZIiCXwaFmR1tZv8KepgbzOyXZlYYLNull21mT5jZF4P7O/Tig8deaGYrgvX9ypLeCfwGODZ4h1HXj1LnBbd1wTqONbPJZvaYmW01sy1m9kczG7FTPQelzc80s2t68Ts5y8wWB/uwwMyOCNo/ZWZvmFlFMH9G0OuuStve183s9aCen5hZn/6X3f1ud3/I3VvcfTtwG/D+tIecD/zA3be7+/Jg+ef6sg0ZuhT8MljiwDeAMcCxwCnARXuxvrOA9wJHAJ8EPhQE1IXAv4J3GCP6sd4Tg9sRwTr+BRhwHW/3jieS7BH3m5m9G7gT+DIwGvgtcL+ZFbn7X4AFwM1mNhq4A/iiu9emreLjQA1wFPBR4PPBeo8PXkh6mo7fzX6/HKxjJDAeWJK2fAlw2N7sswwdCn4ZFO7+vLsvdPdOd19NMug+sBervN7d69z9TeBxYOoAlNktd1/p7nPcPRaE743sXe0AFwC/dfdn3D3u7r8DYsAxwfKvAicDTwD/cPcHdvr5H7n7tmD/fw6cG9T6tLuP2M20y/EPMzuNZA//+0FTWXBbn/aweqAcyQk5fQBLhg4zO5hkYNYAJSSfe8/vxSo3pt1v4e2wGnBmNg64CTiBZPhFgO17udpJwPlmdnFaWyHJdxW4e52Z/R9wGfAf3fz82rT7a7p+rq/M7BjgbuAcd38taG4KbiuAtrT7jf3Zhgw96vHLYLkF+Dcwxd0rgCtIDqEANAe3JWmP36ef29nby8129/M/DNoPD2r/LG/XDskXnr7Wvha4dqfeeIm7/wnAzKaSHL75E3BzNz8/Me3+/sD64OdO2OnMqp2nE7p+KBhuuh/4vLvPTf0CkmP+G4Aj07ZxJMFQkAx/Cn4ZLOVAA9BkZocAX+laEAyfrAM+a2Z5ZvZ5YHI/t7MJmNB14Lg7ZjbDzJ7oYXEtkAAO3Kn2JqDezPYDLt/pZxYDnwlqP53eDQPdBlxoZu8LDkyXBqdYlptZFPgDyRfH/wL2M7Odj4dcbmYjzWwicAnwFwB3fyr9zKpupqeC38G7gIeAi939H93U93vgu8E2DgG+BMzsxX7JMKDgl8HyLeAzJIcLbiMIqjRfIhmoW0keRFzQz+08RrJnutHMtvTwmInA/O4WuHsLcC0wPzgYegxwNcmDqPXAbODvO/3YJcDZQB1wHnDvnop090Uk9/mXJIeNVvL2WTPXAWvd/RZ3j5F8h3GNmU1JW8V9JIfKFgc13bGnbe7km0AVcEfau4H0Hv1VwCqSw0hPAj9x94f6uA0ZokxfxCJhY2aLgVPcfWu2a+kPM3OSQ2Yrs12LDE86uCuh4+5Ts12DSDZpqEckQ8zsih4OsD6Y7dok3DTUIyISMurxi4iEzLAY4x8zZoxXV1dnu4zhraUleVtSsvvHiUjOeP7557e4e9XO7cMi+Kurq1m0aFG2yxjeFi9O3k6dms0qRGQQmdma7to11CMiEjIKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhIyCX0QkZBT8IiIho+AXEQkZBb+ISMgo+EVEQkbBLyISMgp+EZGQUfCLiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjIKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhIyCX0QkZBT8IiIho+AXEQkZBb+ISMgo+EVEQiZjwW9md5rZZjNbltY2w8zWmdniYPpwprYvIiLdy2SPfyZwejftP3P3qcH0zwxuX0REupGx4Hf3ecC2TK1fRET6Jxtj/F8zs6XBUNDInh5kZheY2SIzW1RbWzuY9YmI5LTBDv5bgMnAVGADcENPD3T3W929xt1rqqqqBqk8EZHcN6jB7+6b3D3u7gngNuDowdy+iIgMcvCb2fi02Y8Dy3p6rIiIZEZ+plZsZn8CTgLGmNlbwFXASWY2FXBgNfDlTG1fRES6l7Hgd/dzu2m+I1PbExGR3tEnd0VEQkbBLyISMgp+EZGQUfCLiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjIKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhIyCX0QkZBT8IiIho+AXEQkZBb+ISMgo+EVEQkbBLyISMgp+EZGQUfCLiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjIKfhGRkFHwi4iEjIJfRCRkMhb8ZnanmW02s2VpbaPMbI6ZrQhuR2Zq+yIi0r1M9vhnAqfv1DYdmOvuU4C5wbyIiAyijAW/u88Dtu3U/FHgd8H93wEfy9T2RUSke4M9xj/O3TcE9zcC43p6oJldYGaLzGxRbW3t4FQnIhICWTu46+4O+G6W3+ruNe5eU1VVNYiViYjktsEO/k1mNh4guN08yNsXEQm9wQ7++4Hzg/vnA/cN8vZFREIvk6dz/gn4F/AOM3vLzL4AXA+cZmYrgFODeRERGUT5mVqxu5/bw6JTMrVNERHZsz32+M3s/b1pExGR4aE3Qz2/6GWbiIgMAz0O9ZjZscBxQJWZXZa2qALIy3RhIiKSGbsb4y8EyoLHlKe1NwDnZLIoERHJnB6D392fBJ40s5nuvmYQaxIRkQzqzVk9RWZ2K1Cd/nh3PzlTRYmISOb0Jvj/D/gNcDsQz2w5IiLSHz+b8xo3zV2Rmr/klCk9PrY3Z/V0uvst7v6suz/fNQ1AnSIi0kczZszAzFLTjBkzAPjGaQez+vozAVh9/Zl847SDe1xHb3r8/zCzi4BZQKyr0d13vuSyiIhkWOX7P8Okb783bb7nnn1PehP8XdfWuTytzYED+7w1ERHple6Gbr5x2sGpycxIXuS47/YY/O5+QL/WLCIiu9VTuAOpgK+ePjs1hDNQdvcBrpPd/TEz+0R3y9397wNaiYhIyGQy3Hdndz3+DwCPAWd3s8wBBb+ISC/srmefDbv7ANdVwe1/DV45IiLD157G5Qe7Z9+TPY7xm9n3u2t39/8Z+HJERIa+4RLwPenNWT3NafejwFnA8syUIyIy9A2XgO9Jb87quSF93sx+CjycsYpERIaIoTY2P1D68w1cJcCEgS5ERCRbhvvQTV/1Zoz/JZJn8UDyOvxVgMb3RSRn5GrA96Q3Pf6z0u53ApvcvTND9YiIZEyuDt30VW/G+NcAmNlYkgd39w0+KvxmposTEemPsA3d9FVvhno+AtwA7AtsBiaRPKvnsMyWJrJnPf2D97VdcosCfvd6M9TzA+AY4FF3f7eZTQM+m9myJKz6Gtg9/YP3tV0vIMOT/g7905vg73D3rWYWMbOIuz9uZj/PdGGSGzId5AMl0y8gsnc0dDOwehP8dWZWBswD/mhmm4GmzJYlw0l/rjCYq/+wu9svvSj0X64+X7KlN8G/BGgBvgGcB1QCZZksSoYm9br2joaZ9iyX920o6U3wT3P3BJAAfgdgZkszWpVklQJ+cGXrOEWm9acePccGx+6ux/8V4CJg8k5BXw7Mz3RhknkK+NySrQPa2Tpe4+7E43E6OjpIxFrYvn07nZ2dO0wd29axfPly4vF4qi0ej9P21ivMm1eeao/H48TjcVpWPMOsWe2p+UQiQTwep+mlF7jrrs0kEolUWyKRoOH5pdx886pUu7uTSCSoX/gK11//Uqqtq73u6VeZMeO5VFvXtH3eCq68ckHqG7XSv1lr+5Mr+c53nu7376k7u+vx3w08CFwHTE9rb9T37Q4vP5vzGnPufgqAV8atU8AL0HMAX3rqFC46cRJTpt/Hou+cRCwWY8WKFZw2vp0TzhnHWT9/nLs//x5isVXce+8r7BuL8YPD2rn07uf4wdmHEFv2Gj96YRaxWIz29na2P76cSy+dQ3t7+w7T5hff5OyXf0N7ezsdHR2p246ODtav3co7ZkVT813T9qZWSm6Gjo4OOjs7iUQiFBQU0O4RDpxZRH5+fmrKy8tjc0M7H59XsUNbXl4edesauXLN/Tu05eXl0fTaFn4XX7pDWyQSoW3NBubN204kEklNeXl5dGx9i5UrC4lEIphZalmitYG6urrUF6J3LceTLwTpbWaG5RdSUlICkHxc2m2ksISKiooB/dvv7nr89UA9cO6AblEyZrc9r6oWPnzTUwr4YSaRSJBob2PLli20tLTQ0tJCa2tr8vaNF7n//jitra07TPX/WsL3vrcwNd/W1kZrayubn1/NGUt+SVtbW6qtra2NtzZsZ5+ZRltbG7FYjFgsRn5+Pp2Wz2F/KCUajVJUVERRURHRaJS6jS38YN0/U22FhYUUFRURW7eZZcs6U/Ndt5HiCqqrqyksLNxhmhtbygWfP5aCggIKCgooLCxM3T/71wu57/JTUvNd01HXPsa/rz0rFeSRSASgx85L9fTZ/LuH9qd6aL+3h/a7umn/5/TZ3NxN+/9Nn8313bT/b8dsrr561/Y7mmdz5ZXd/2/+tn423/nOrsuuuOKKbh8PsHZbC29ua+lxeX8u0iZZpiGaocPdaWlpobm5maamJpqbm2lubqZ19WLuvz+emk+fts19mQsuuC813/Xz61du5B2zLksFfEtLC7FYDPIKOeT3ZRQXF1NSUkJJSQnFxcU0rG/mtpZnKC4u3mFKtLdSWFhIZWUlxcXFRKNRiouLebTtZb7+peOJRqOptmg0yik/n8+iqz6cai8qKiISiew2TB/vof3X3bTf0TybSy/dtX364nLOPrv752lh1UYOOeSQXdojRaWpnnFYuTt1LR0U7nMQDyxdz5vbWlJBv3ZbKwCfvnUhE0YW97gOBf8wpIDvv/b2dhobG+mo28hLL71EY2PjDlPD889w3XVLaWxspKmpiaamJhobG9n0/BucMP/6VFtqam5h9M2FlJWVUVpamrqt39jGbc0LKS0t3WXKr6iipqaGkpISSktLU2H+yTte4L7pH0rNl5SUEI1GOfCKB3sM4H900z5r+my+971ueohLZ3PGGWfs0l4wYhXjx48fmF+w7LVYZ5yN9W2sq2tlfV0bANP/tjSYT7bl5xmjPvQ1Zi/dwP6jSjhs30pOf9d49h9VwrSfPsH86ScDcM+F3W9DwS9DnrvT1NREQ0MDHVvWsnDhQhoaGnaY6uYv4lvfepyGhoZkgAft6199iwPvuTjVlkgkKC8vpyGez7lPj6W8vJyysjLKy8spLy+nY+tW6utHUF5ezrhx41LLnuQVrr34ZMrKynaYDr/mCdb8+CO71NxTKAPc3jSbCy7YdVnR7IZue7mSW7Y2xVhf15YW5Mle+kd/NZ/1da3Ut3QwtqKIfUcUs9+IZK/98AmVfOhd+7DfiGLGV0YpjxZgdjq3zPTdbapHCv4hLBfOaXZ3Wltb6WzcyvLly6mrq6Ouro76+vrk7cJnuOKK+W/P19dTX1/P+lffovrPX6W+vp7GxkaKioqorKykNpbHJS/tS0VFRWoqLy+HeCfjxo1jypQpqbaKigo+fddi5nz3zFSwR6NRzIzq6bNZ1k0wP9nD2Ox/vzCbE088cZd2i+Rl5Pcmw1NnPMGG+mQv/Z5Fa1PBvr6uLRXwp9z4JPtWFgfBHmXfINy/f9ah7DeimKryIvIillrnrBfXcd77Jg1onQr+IWyoDOl0dnayfft2Nm3ZyqYtdbS9+RI/v6uZLfVNbGtoZntTKw0t7TS9voX3/L8FtHYkaOuEdjc6PI9IYTHkF3HmjXPJKygkkl9IJL8AyxtJ2ZEf4j4vgsoIPiJCAiOBMf64BCWF+ZRFjLyIkReJEDEobGrHKqI0G7RGjFozIgYjSlqYF60k6hGiLXkUtedR1BChcNxk/vhKjPJonPJoE2VF+ZRHCwB4eX09o0uLGFlaQFG+Alx6J9YZZ+22Vt7clvxW2qvuW8aabS2s2drCurpWqsqKAHjm9W3sNyLKUfuP5MwjkiF/6o3zWPz9D+6yzuse/DfvmTRy0PZBwT8EDEbP3t2Jm7EueCvZ0NZBfWsHq95cT+vqxXz426tpbIvT3OG0xSHm+XR4PhMvmgkFRUQKSzADi7dTMHoiv32pncJIAdH8kZQUjaasIp+i9rG8d+oERpaXMLK8hDGVZYwZUUZlaTHn3/ksf7nwWPIjRkFehPw8Iz8S4dQbn+Th/56WbIsY+cHtYVc9zDNXnkrCIZFwEu4kHN577aPM+upxxBOOOyTciSeck294kms+9i7aOuK0dSaStx1xZi/dwL4jojS2dbKhvo3Gtg6aYsmvk/jmPUvY1tzO9pZ2ivLzGFVaCMAXZj7HqNJCRpUVMqqkMNX+740NjK8opqI4P3WqneSmRMJZV9fKytomVm1OXqHmM7ctZM3WFmobY+w7Isqk0aUA7D+6lBMPrmLS6BImjCwhWpBH9fTZ3PDJI7O5C7uVleA3s9VAIxAHOt29Jht1DBUD0bNvbY/z+pYmVtUmeyEz7n+ZzY1tbGqIsbmxjc31rbSPmcQ5tyygsriAimgBFcUF0N6cPIe4uIT9xhQxuqKYqsoy9hlVwbfve5WHv30640ZVUhbNT/WKq6fPZnkPBxt/c9GHe6zxvdWjum2fMLL7szS6euY7G1/Z/dkKR04csUvbJX9ezBdPOHCX9urps3no0uTQjbvT0NbJtuZ2pv30CT599P5sa46xtbmd2sYYr25sBODiu19kY30bcXf2qYwyvjIKwK8eX8mk0SVUjy5l0uiSHuuWoacjnmDN1hZWbm5kZRDwZ/3iKVZtbqayuICDxpZx0NjkFWou/MBkqkeXsu+IKPl5b59G+oXjD8ha/f2VzR7/NHffksXtD0vxhLNycxOL124H4P/d+SyrNjexpSnGpNElTK5KPkknjiqhpnokY8ujjKsoYuwbr/Ifv5rPP3988S7rrF7VxF9n7BrkVz62hUP3H5vZHRoCzIzK4gIqi5OBfdqh43Z5zN9fXMecyz4AQGNbBxvr29hQ38b8lVtpbOvkgSUbWL21mTVbWygpTL5AXnbP4tSLQfXoUqpHl1JZoheFbGhtj7OqtolVtU2pgD/1xidZu62F8ZVRDhpbxuQg4K/52OFMrird4QV85oLVnHhwVVZqzwQN9Qyi/gzpbG5MHii67sHlLFlbx7J1DVSVFzE16N1+7rhJTK4qY8LIktQBoW57IWsjA7cjIVceLaA8WsCUceUATD/j7TNx3J3axhhH/3Auxxw4mjVbm3nk5U2pF4Wuv9F/3fUs+1QWs09F8p3DuODdQ31rBxVRDSX1V11LOys37xjwJ/z4MTY3xKgeXbpDD/4X576bA8aUEi14+/jOb598PfW/lcuyFfwOPGJmDvzW3W/NUh2DqjdDOm0dcRat3s68FbXMe602dYZAaWE+XznpII6cUMmIkuSY86wX13HyIbv2TiV7zIyxFckQ/2TNxB2WuTvbmtt5zzWPct77JrGxoY2N9W08u3obG4O/8/uvf4yEO/tURNmnMso+wbpuf+p1qsqLGFOWnKrKiwZ3x4aQncffu4Y33/ODObR3Jpg8tozJVWVMHpscg//959/HxJHFqeEZgJvmruCd4wf2MgjDSbaC/3h3Xxd8j+8cM/u3u89Lf4CZXQBcALD//vtno8ZBs6mhjQdf2gBAzTWPcvC4Mk48uIprP344R06o5KArH+Trp0zJcpWyt8yM0cEZH6d2M5xUPX02y67+EI1tHWxqaGNjfYwN9a38/cV1rK9rY8lb9WxpjLGlKUZtUwyA9/3w0R1eDMYE679v8Tqq0tpGDMMhplW1Tby1vZW3trfw1vZW1gaXIDjsqocZUZIcf59cVcah45PvvB685ASqyot2eLf044de5YAxpVmpfyjLSvC7+7rgdrOZzQKOJvlFL+mPuRW4FaCmpqZ/n1LIkt4M6WxpivHgso08sGQ9yzc0pIJg/rdP1jhwyHUNJR00Nhlol/91Kd8/+9BdHlc9fTb3fvX9bGlsp7apLbhNviDMeWUTtcGLxJamdlrak2cynXzDE6kD+xXR/OQBfpIHqNPbKoLx7bXbWigtyqekMI+i/Ei/h6Aa2jpoauukKdZJY1sHjW3Jeu58+g1qm2JsaUy+mNU2xtjcmNyHz898jokjS5gwspgJI4s57dBxPLB0A89991TKinaMru/d93LqnZbs2aAHv5mVAhF3bwzufxD4n8GuI5N6GtJpae/kgaXJnv20nz7BtHeM5QvHH8CJB1cRLcjj7y+sU+hLn4yvLA7OcqpMtf3k4Vf55WeO2uFxsc447/juQ9z6n++hoa2ThtaO1C1AY1sn6+ta05Yl2z9960Ja2jtpbo/TGU9+tgLgmB/OJS9iFORZcJscRjn5hieIdSSIdSaIdcaJdSYAOO66x4LPUORTFn37sxRvbmthTFkh1dUjqSovoqosSlV5EcdcN5cnL5+2y/5e8ufFu4S+9F02foPjgFlBzyEfuNvdH8pCHYNm+YYG7n7mTe5fsp73Vic/pPHsFadSXKgPDcng6DoVt+tdRLrv3rtshwPUXaqnz05d8wWSn0pt6YhzxIxHmPXV4+iMO50JpzOeoDPhnHHTU9z6nzUU5UeCKY+iggiHfO8hll39oW7XP+Mjhw3gXkpvDXrwu/vrwND9ZMMAaW2PA/DxX89nQ10bn3rvRB685AT2HVFM9fTZCn0ZdvLzIlQEPfuePkvRdcaMDG16z7QXuhvLP/+4amYuWM0fFq4B4KKTDmLaO6p2OKNgMOqpnj57h29I6voilp3bu3u8iOQ2Bf9eSB/LXzD9ZG576nWm/fQJPnz4PvztK8cx7adPdPthoN7aXZB31941dVtnVfClDFOn7lL/3m53dy8genERGXoU/HtpVW3yQyJn3PQUn6yZwMOXnsg+lX07u6BfQZ7B8OzrdndXz2C8uIiEyYwZM7j66quB5CnCV111FTNmzOjTOhT8/bS5oY0b57zGnFc2AfDk5SelPljVk74GfC4bqBcXvVBI2MyYMYMZM2bs1bW9FPx7sPM4/kUnTaYoP4+7FrzBp2om8ti3TuLIqx/ZY+hD5nvqYTRQLxQiYaLg34P0cfwfn3MENz7yGjXVI/nH145n4qjuryqpcBm69IIgouDvlSVr6wD4y3Nr+fVnj+Ko/Xf/hQnq2Q8/GkqSMFHw70asM87Nc1fwl+fWAvDXC4/d4SPr+ufPfRpKklyk4O/BsnX1fPOeJUwcVcI/v34CR/9w7i7XKVHPXna2u+eEXhRkqFDwB3Y+iFtckMcPP/EuPjZ1P10bXQaE3iXIUKFv5wh847SDeTT4hqVTDhnLE5efxMffPYGfP7qC6umzgeQ/5c/mvJbNMiUHfeO0g1l9/ZmpKf3Db90993pql+Fnd3/Lvv79+/K8UI8/MOeVTUz/21IAbj+/JtXL13COZIuOL2TfQH3QsD+f4RmIz7pc1sN+hT74Ewnn5sdW8Odn13Lb+TV84tcLNLQjw5LOTHpbpoN5ID/dng3mPvS/46SmpsYXLVo04OttbOvgsnuWsK25nXdPHMHtT7+RWpZL/wQALF6cvE27Vo9Id3r6IqG+tnfp6ROm3bUP9LbDzsyed/eaXdrDFvw7P1EO36+Sv33lOArzc/xwh4JfBpnCOvsU/Gnmr9zCebc/w3WfOJxzj87t7/NNUfCLhE5PwZ/j3dxdPbRsI1//04sA4Ql9EZE0oQr+vz3/Fpfds5itze2AToUTkXAKzVk9M+e/wa3zXuf+rx2vr4cTkVDL+eB3d26eu5J7F6/jnguPZcLI7q+oKSISFjkd/ImEc83s5SxYtYV7vnwsVeVF2S5JRCTrcjb4b3zkVW5+bGVq/g8L1+iUMRERcvjg7nOrt2e7BBGRISkne/x/WLiGdXWtPHvlKYwt79sXn4uI5Lqc6/HPXrqBXzy2gv/9wtEKfRGRbuRU8D+9Ygvfv28Zd33uaCaNLs12OSIiQ1LODPUsWVvH1//8IrecdxSH7luR7XJERIasnOjxv17bxBd/v4gf/ccRvO/A0dkuR0RkSBv2Pf5zblnAojXJM3i+9PtFvO+AUfzly8dmuSoRkaFrWAe/u7NPZZRP1UzkR+ccke1yRESGhWEd/HfNX83qrc389cLjsl2KiMiwMWyDf9Hqbfz6iZXMuuj9RAvysl2OiMiwMSwP7tY2xvja3S/yk3OOZOIoXXRNRKQvhl3wd8YTXPynF/hkzQSmHTI22+WIiAw7wy74f/LIqxTkRbjkVF1wTUSkP4bVGP9DyzbywJIN/OPi48mLWLbLEREZloZN8H/qt//imTe2AXDUD+bofH0RkX7KSvCb2enATUAecLu7X7+nn1HIi4gMjEEf4zezPOBXwBnAocC5ZnboYNchIhJW2Ti4ezSw0t1fd/d24M/AR7NQh4hIKGUj+PcD1qbNvxW0iYjIIBiyp3Oa2QVmtsjMFtXW1ma7HBGRnJGN4F8HTEybnxC07cDdb3X3GnevqaqqGrTiRERyXTaC/zlgipkdYGaFwKeB+7NQh4hIKA366Zzu3mlmXwMeJnk6553u/vJg1yEiElZZOY/f3f8J/DMb2xYRCbshe3BXREQyQ8EvIhIyCn4RkZBR8IuIhIyCX0QkZBT8IiIho+AXEQkZBb+ISMgo+EVEQkbBLyISMgp+EZGQUfCLiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjIKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhIyCX0QkZBT8IiIhY+6e7Rr2yMxqgTXB7BhgSxbLyQbtcziEbZ/Dtr8w+Ps8yd2rdm4cFsGfzswWuXtNtusYTNrncAjbPodtf2Ho7LOGekREQkbBLyISMsMx+G/NdgFZoH0Oh7Dtc9j2F4bIPg+7MX4REdk7w7HHLyIie0HBLyISMsMm+M3sdDN71cxWmtn0bNeTKWZ2p5ltNrNlaW2jzGyOma0Ibkdms8aBZGYTzexxM3vFzF42s0uC9lze56iZPWtmS4J9vjpoP8DMngme438xs8Js1zrQzCzPzF40sweC+ZzeZzNbbWYvmdliM1sUtGX9uT0sgt/M8oBfAWcAhwLnmtmh2a0qY2YCp+/UNh2Y6+5TgLnBfK7oBL7p7ocCxwBfDf62ubzPMeBkdz8SmAqcbmbHAD8CfubuBwHbgS9kr8SMuQRYnjYfhn2e5u5T087fz/pze1gEP3A0sNLdX3f3duDPwEezXFNGuPs8YNtOzR8Ffhfc/x3wscGsKZPcfYO7vxDcbyQZCvuR2/vs7t4UzBYEkwMnA38N2nNqnwHMbAJwJnB7MG/k+D73IOvP7eES/PsBa9Pm3wrawmKcu28I7m8ExmWzmEwxs2rg3cAz5Pg+B0Mei4HNwBxgFVDn7p3BQ3LxOf5z4L+BRDA/mtzfZwceMbPnzeyCoC3rz+38wd6g7B13dzPLuXNwzawM+Btwqbs3JDuDSbm4z+4eB6aa2QhgFnBIdivKLDM7C9js7s+b2UlZLmcwHe/u68xsLDDHzP6dvjBbz+3h0uNfB0xMm58QtIXFJjMbDxDcbs5yPQPKzApIhv4f3f3vQXNO73MXd68DHgeOBUaYWVdnLNee4+8HPmJmq0kO1Z4M3ERu7zPuvi643UzyBf5ohsBze7gE/3PAlOAMgELg08D9Wa5pMN0PnB/cPx+4L4u1DKhgnPcOYLm735i2KJf3uSro6WNmxcBpJI9tPA6cEzwsp/bZ3b/j7hPcvZrk/+9j7n4eObzPZlZqZuVd94EPAssYAs/tYfPJXTP7MMkxwjzgTne/NrsVZYaZ/Qk4ieTlWzcBVwH3AvcA+5O8PPUn3X3nA8DDkpkdDzwFvMTbY79XkBznz9V9PoLkQb08kp2ve9z9f8zsQJK94VHAi8Bn3T2WvUozIxjq+Za7n5XL+xzs26xgNh+4292vNbPRZPm5PWyCX0REBsZwGeoREZEBouAXEQkZBb+ISMgo+EVEQkbBLyISMgp+kT4ws5+Z2aVp8w+b2e1p8zeY2WVZKU6klxT8In0zHzgOwMwiJD9vcVja8uOABVmoS6TXFPwifbOA5OUVIBn4y4BGMxtpZkXAO4EXslWcSG/oIm0ifeDu682s08z2J9m7/xfJK0oeC9QDLwWXDhcZshT8In23gGToHwfcSDL4jyMZ/POzWJdIr2ioR6Tvusb5Dyc51LOQZI9f4/syLCj4RfpuAXAWsM3d48EFtkaQDH8Fvwx5Cn6RvnuJ5Nk8C3dqq3f3LdkpSaT3dHVOEZGQUY9fRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZD5/wMvkbuX2lVYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c_obs3.gamma_method(tau_exp=20)\n",
    "c_obs3.print()\n",
    "c_obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jackknife"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison and as a crosscheck, we can do a jackknife binning analysis. We compare the result for different binsizes with the result from the gamma method. Besides the more robust approach of the gamma method, it can also be shown that the systematic error of the error decreases faster with $N$ in comparison to the binning approach (see hep-lat/0306017)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binning analysis:\n",
      "Result:\t 3.27194697e-01 +/- 1.81819841e+00 +/- 3.98347312e-01 (555.693%)\n",
      "Result:\t 3.27194697e-01 +/- 1.66475180e+00 +/- 5.21149746e-01 (508.795%)\n",
      "Result:\t 3.27194697e-01 +/- 1.41273466e+00 +/- 6.28627238e-01 (431.772%)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABHLUlEQVR4nO2debwcVZn3f8/d9/0muVlvVpKwBbhsCsgWQEDR0XFwHJcZEHF5R3D7BMdXHed1BseFYUYHJioDKuLCgDIGhRB2JUBAAiH7Svabu+/78/5RdapPd5/qrq6uXu/z/Xzu53ZXVVef6q4+z3l2YmYIgiAIQqIUZHoAgiAIQm4iAkQQBEHwhQgQQRAEwRciQARBEARfiAARBEEQfCECRBAEQfCFCBBBSDNE9CEiejzT4xCEZCHJAxEEd4hoQHtaAWAUwKT9/BPMfH+K3/9jAG5k5gtS+T6C4IeiTA9AELIZZq5Sj4loP6zJ/InMjUgQsgcxYQmCD4joHCJ6gYh6iOgoEX2fiErsfa1ExERUpB3/NBHdaD/+GBE9r+1jIrqZiHbZ5/sBWawAcDeA84logIh60nyZghATESCC4I9JALcCaAJwPoDLAHwqifNdC+BsAKcB+ACAK5l5G4CbAbzAzFXMXJfUiAUhYESACIIPmPkVZt7IzBPMvB/AfwF4RxKnvJ2Ze5j5LQBPAVgVwDAFIaWID0QQfEBEywB8D0AbLOd6EYBXkjjlMe3xEIAqtwMFIVsQDUQQ/HEXgO0AljJzDYAvAyB736D9v0I7fpbP95EwSSFrEQEiCP6oBtAHYICIlgP4pNrBzCcAHAbwN0RUSER/B2Cxz/c5DmCuctALQjYhAkQQ/PEFAH8NoB/ADwH8MmL/xwF8EUAngJMB/Mnn+zwJ4E0Ax4iow+c5BCElSCKhIAiC4AvRQARBEARfiAARBEEQfCECRBAEQfCFCBBBEATBF3mVSNjU1MStra2ZHoYgCELO8Morr3Qwc7Of1+aVAGltbcWmTZsyPQxBEIScgYgO+H2tmLAEQRAEX4gAEQRBEHwhAkQQBEHwhQgQQRAEwRciQARBEARfiAARBEEQfCECRBAEQfCFCBBBEATBF3mVSCgIQjh3rN+JOzfscp5/9rKluHX1sgyOSMgn8qofSFtbG0smuiBE07pmHfbffk2mhyFkIUT0CjO3+XmtmLAEQRAEX4gAEQRBEHwhAkQQBEHwhQgQQRAEwRciQARBEARfiAARBEEQfCECRBAEQfCFCBBBEATBFyJABEEQBF+IABEEQRB8IQJEEARB8IUUUxQEA1KEUBDiIwJEEAzcunoZbl29LOEihCJ4hOlEykxYRHQPEbUT0RZt29eJ6DARvWb/Xe3y2quIaAcR7SaiNakaoyAEza2rlzkCZ//t14jwEPKaVGog9wL4PoCfRGy/g5m/4/YiIioE8AMAqwEcAvAyET3CzFtTNVBBSAWijQj5Tso0EGZ+FkCXj5eeA2A3M+9l5jEAvwBwXaCDE4Q0INqIkO9kIgrrM0T0um3iqjfsnwPgoPb8kL3NCBHdRESbiGjTiRMngh6rIAiC4EK6BchdABYDWAXgKIDvJntCZl7LzG3M3Nbc3Jzs6QRBEASPpFWAMPNxZp5k5ikAP4RlrorkMIB52vO59jZBEAQhi0irACGiFu3pewFsMRz2MoClRLSQiEoAXA/gkXSMTxAEQfBOyqKwiOgBABcDaCKiQwC+BuBiIloFgAHsB/AJ+9jZAH7EzFcz8wQRfQbAYwAKAdzDzG+mapyCIAiCP1ImQJj5g4bNP3Y59giAq7XnjwJ4NEVDEwRBEAJAamEJgiAIvhABIgiCIPhCamEJQhpoXbPOeSwZ6UK+IBqIIKQByUgX8hERIIIgCIIvxIQlCDmAFGYUshHRQAQhB5DCjEI2IgJEEARB8IUIEEEQBMEXIkAEQRAEX4gTXRDyAHGyC5lABIgg5AG3rl6GW1cvQ+uadY6zPRIRMkLQiAlLEKYJEsklBI0IEEEQBMEXIkAEQRAEX4gAEQRBEHwhAkQQBEHwhQgQQRAEwRciQARBEARfiAARBEEQfCECRBAEQfBFygQIEd1DRO1EtEXb9m0i2k5ErxPRw0RU5/La/UT0BhG9RkSbUjVGQRAEwT+p1EDuBXBVxLb1AE5h5tMA7ARwW4zXX8LMq5i5LUXjEwRBEJIgZQKEmZ8F0BWx7XFmnrCfbgQwN1XvLwiCIKSWTBZT/DsAv3TZxwAeJyIG8F/MvNbtJER0E4CbAGD+/PmBD1IQhNSQj8Ud8/GaYpERAUJE/wBgAsD9LodcwMyHiWgGgPVEtN3WaKKwhctaAGhra+OUDFgQhKRwm1jjVRDONfLxmmKR9igsIvoYgGsBfIiZjRM+Mx+2/7cDeBjAOWkboCAIgSOVgPOTtAoQIroKwJcAvJuZh1yOqSSiavUYwBUAtpiOFQRBEDJHykxYRPQAgIsBNBHRIQBfgxV1VQrLLAUAG5n5ZiKaDeBHzHw1gJkAHrb3FwH4OTP/IVXjFAQhOaab3V8IkTIBwswfNGz+scuxRwBcbT/eC+D0VI1LEIRgmW52fyGEtLQVhCxAVvFCLiICRBAyhElo3Llhl6zihZxBBIiQVmSlHcJk+tE/G0HIdqSYopBWJJxTEPIHESCCIAiCL0SACIIgCL4QH4ggZBniJxJyBdFABCHLED+RkCuIBiIIQlYimlj2IxqIIAhZiWhi2Y8IEEEQBMEXYsISBCGttK5Z5zz+7GVLMzgSIVlEAxEEIa2IWSp/EA1ESBpxdgrC9EQ0ECFpxNkpCNMTESCCIAiCL0SACIIgCL4QASIIgiD4QgSIIAiC4AsRIIIgCIIvJIxXEIS8QsLK00dKNRAiuoeI2oloi7atgYjWE9Eu+3+9y2s/ah+zi4g+mspxCkIuccf6nWhds875u2P9zkwPKauQsPL0kWoT1r0ArorYtgbABmZeCmCD/TwMImoA8DUA5wI4B8DX3ASNIEw3ZIIUsoWUChBmfhZAV8Tm6wDcZz++D8B7DC+9EsB6Zu5i5m4A6xEtiARBEIQMkgkfyExmPmo/PgZgpuGYOQAOas8P2duiIKKbANwEAPPnzw9wmIKQ/4i/QEiGjDrRmZmJiJM8x1oAawGgra0tqXMJ2YtMdKnh1tXLcOvqZWhds84xiwmCVzIRxnuciFoAwP7fbjjmMIB52vO59jZhmiJ2f0HIPjKhgTwC4KMAbrf//9ZwzGMA/llznF8B4Lb0DE8QUodJkxKEXCXVYbwPAHgBwElEdIiIboAlOFYT0S4Al9vPQURtRPQjAGDmLgD/BOBl++8b9jZByGlEkxLyiZRqIMz8QZddlxmO3QTgRu35PQDuSdHQhDTiZdUtPg5BCI50/Z6klImQcrysumVlLgjBka7fkwgQQchz9p4YyPQQhDxFBIgguHDvH/dlegiB8PjW45kegpCnxBUgRPR2L9sEId94fndHpocQCKPjU5kegpCneNFA/sPjNkHIK4bHJ329rmdoLOCRJMfohL/rEIR4uEZhEdH5AN4GoJmIPqftqgFQmOqBCUKmGRrzN/F+/lebAQAv7OnE+YsbgxySL0YnLA1kYHQCVaXSwWG60bpmnfM46GisWHdTCYAq+5hqbXsfgPcHNgJByFKGfQqQ/pEJAMAXfr0Z5y5sCHJIUUSGa5pQGsiBzkGcPLs2peMR0ouXcN39t1+TslI1rgKEmZ8B8AwR3cvMBwJ/Z0HIcvyasEbsCXv95y7CR378UpBDiiKylpW+2lQoH8hbnUPTToDke35RpmuZefGBlBLRWiJ6nIieVH8pH5kgZBi/GsiILXgqSorw2cutpMmpqczV+VQmrANdQxkbQ6aQ/KLU4sUg+msAdwP4EQDxxgnTBt8aiBb1dMGSJgDAhu2mmqHpQQm0A52DGRuDkJ94ESATzHxXykciCFlGshoIABARAOA/n94dyJj84GggnfmhgeS7WSqX8GLC+l8i+hQRtdj9zBvslrOCkLeMT05hwqfZacSgufQMjSc7JN+EnOj5IUDELJU9eNFAPmr//6K2jQEsCn44gpAYqVqN+g3hBYCRiejEvU9ctAhrHnojmSH5RmkgJ/pHJSdECJS4AoSZF6ZjIILgh1RFoZi0CC9MTTHGJ6MFyGUrZgLIkACxfTItdWU41D3s6TViJsocufTZx0okvJSZnySivzDtZ+aHUjcsQcgsfv0foxNTKCkscFb9iobKEgDA5BSjsICSHl9iY7KuZX5DBd7yaMZyE8y5NLnlKpkOzU2EWBrIOwA8CeBdhn0MQASIkLf4NWGNjE+irLgwSoAoodE9NIamqtKkx5cIaiwLGiuwP8lIrFya3CIR4Rc8sRIJv2b//9v0DUcQsoPh8UkUFVDCjvSRiUmUFReg18VS1DngX4A8s/OEr9cpAdLaWOk40r1ksOcbuSz8spW4PhAi+qppOzN/I/jhCLlCvvf2Hh6bRF1FMToGEiuMODI+hbJi91JxnQOjCK8M5J097f76eoyOh0xYf9rTCcBbBrsgxMNLFJau85YBuBbAttQMR8gVTKu5XFnRejFlDI9Poq6ixIcAmURZUQwBMui/Uq//0irKhFWZtAlrupPNZrBMLOq8RGF9V39ORN8B8FjKRiQIKcaLKWN4fBL1FcUJn9vygbinV1kaiD/8OPaZGWO2AJnfUOE5Cssv2TzBBkE2m8Eysajz05GwAsDcoAciCNnE8NgE6ipKEn7dyPgUSmOZsJLQQJRjf8yQZwIAE4bw4dGJKZQUWT/z8pJCX0IxESTJb3rhpSPhG0T0uv33JoAdAP7N7xsS0UlE9Jr210dEt0QcczER9WrHGP0wgpAqhsd8aiATk7F9IEmZsKwy8e39I8b9Xfa5dUEyOjGF0qLQz3xBQ6Xv9xeESLz4QK7VHk8AOM7ME37fkJl3AFgFAERUCOAwgIcNhz7HzNcatgtCyhmyfSCAZQZSNa3iMTo+ibKi1JqwjveNYG59RdR+5a8ZHp9EdaE1htGJSZQWFaIf1k92Tn05sN/3EIQAyCczX1wNhJkP2P1AhmF1IpxNRPMDev/LAOyRfiPTix89tzfTQ4jLyNgkKkosTSIypyPm6+JGYSVvwjreZxZCHbZw0p3to+PhGki9D7OcECz5ZObzEsb7bgDfBTAbQDuABbCisE4O4P2vB/CAy77ziWgzgCMAvsDMb7qM7yYANwHA/PlByTUhlby8vytt7+V3tTc8Punka4zGEQo68ZzoXQFEYR3rNZuwOgctATIyFmHC0sbTUJlaH4gwvfBiwvonAOcBeIKZzyCiSwD8TbJvTEQlAN4N4DbD7lcBLGDmASK6GsBvABhj0ph5LYC1ANDW1pa5rj2CZ3b5zGfwg9+omaGxSZQ7GsgkAG8Tr8pEd6MjIBOW8dz9IROWQpmwFH4CA9LNdExyzFW8RGGNM3MngAIiKmDmpwC0BfDe7wTwKjMfj9zBzH3MPGA/fhRAMRE1BfCeQoYZnZhMWVnx1jXrnL871u9M6lzD45MotwWB3iAqHiMTsbWVobFJ1yiqeAzFEyC2BjI0FnJRRjrRVU2ubCbSxJMK7li/M9D7ZbriRQPpIaIqAM8CuJ+I2gEEsYT8IFzMV0Q0C5aznonoHFiCrjOA9xQyzP6OIUymqL2ryqgOYtIZGY/UQLy/LpYTvb6yxLcZyzFhJaCBjIxPhgmQuhSH8eYK2ZzPkUt4ESCbAQwBuBXAhwDUAqhK5k2JqBLAagCf0LbdDADMfDeA9wP4JBFNwHLeX8/MYp7KYu77035Px+1q78eM6lK09/s35aSDIc2JrjQQL/6UkfEpVJe5/6waK0scX0XiY7I0CzcnuuMDCTNhhWtEuaCBCLmDFwFyCTNPAZgCcB8AENHrybwpMw8CaIzYdrf2+PsAvp/Mewjp5bevHfZ03K7jAzhzfj3+8OaxFI8oOYbHQr4MpYF4WbWOjE+iudq9WGJTVanvSCzdB2JaT6nz6pWEJQpLSCWx+oF8EsCnACyOEBjVAP6Y6oEJuYXXEhm72wdw5oK67BcgPn0goxOxo7AaIjSQLYd7AQAf/8kmdMcxbSnTVAER+kejU7GcMN6xCCd6se5Et0xYieS2uKEXYMy3YprZRrbmjsTSQH4O4PcA/gXAGm17PzOnLw5TyHpGxic9m6R2tffjb9/emtoBBcDw2CQqSqyfRyLdCUfGp2IWU2ysKgnTQO56Zg8A4H1nzsHQ2CQ2HejG/o5BtDaFZ4yPT045vqMZNaU4bgjlVeeNNGHpGoiKyBoYnUB1Wbg/JNFifJE+p3iRU/s6pJCjX7LVZ+O6VGLmXmbez8wfVMmE9p8IDyGMIz3eC/Qd6BzCKXNqAcDY+jVb0DWQxBIJY4fxNlWVhpUz2XakDwBw1Skt+IszrRJzN9z3MnqHx6PGowTarJoyox9ECYowE1aEAFF0D45HbUsmwS1eZBkz48o7no3a/r+bj3h+DyH78OIDEXKMdKu7BxOo8NpSW+ZMsP0jE1nr1B3W8kAS00Dim7AO2CXVB0cncMTQeeqCJU34+wf+bBzPwOgEZtWUGSOxGqtK0D86EZGJHp4HougeGsP8xuhyKH55LI5Jcnh8EmOGBcNv/uzNdyZkJ36q8QpZTrpLJRzqHsKiJm9F+pbMCDVT6h+JXgVnC8NhYbyJ+EBi54E0VoZMWNuP9WOp9nko/u+1K6Oy9fWosBk1ZcZcEJU5PxxpwjIItK4h/xnxJk7EMWH2j0T7bIbGJvD87o5Ax5FNTIdcE9FAhKQ51D2MU+bUYq8HG/fSmaEIcNOkopOprofMHOFED04DadRMWFuP9mFlSw3esB3piqLCAtSUFYeZoobGJpzxzKopxZ4T0Z91Y5WlzQ17MGH1BCxA4uW29I+Mo6W2DEd7R9BnLxye3dmBipLo/vG5RuR9esPbFwLIXr9FkIgGIiTNoe5hnGr7NeKxdEZIgPTF0UAyVXRudGIKxYUFKCwg57lXRsanjCYjhZ4HsvVIH1bOrjEeVxWRS6InNs6qNWsgjUoDGYtvwuoy+ECSQQnFKZck0b6RCcyww5uf22lpHeu3Hse7T58d6DgyQeR92to8fUrmiwARkuZgV8gxHi/fc2mYCct3V4CUMjwW0j6ABDWQOP1A9CisrUfdBUhkMmIyJiyTRhS8BmIJRTfTWP9IKOprw3aretFTO9rx3jPzrzddTxIFM3MNMWEJSXOoexiL7FVX38gEasujy2WoJkeLZ4RWZ+kWIKYifSYz2QfOnudM1kCCPpBx84StqCotwoS9St91vB/LZ0X7QNRxOkOaUHNzojfZJqyRGMUUFd0pMmG19406gkynf2TcEYpP7zgBAJhZU+Z6/X7IliKMQfuXshnRQISk6RsZR7M9abglw6lILRWKCiTuRP/u4zt8jtDCVKTPZCZLSgOJE8ZLRGi0I8+aq0ujcjEUNRHbrSisIud1Jp+DmrhjZaIrTGG8yaBMWG7dEi0NxBq/MmWtXjkzKts/GdJRhNELPUPZGxwSNKKBCEkzp64cBba/4K6n9+CXmw46+1QI8R5DCfdIDWTvidg1Once7/c0nmSd73oZEyDYPBDAMmMd7R3Byhaz+QqI1kCGxydRYZ+3uLAAteUlUaXhlWDyEoUVtAaiFg5u0ViWBmIJxUuXz8D2Y/24YuVMbf8ESquiPzeVXxJE5ny6UMI9l8bsF9FAhKSZW1/uPL7qlFlhq787N+xC65p1uP/F6KaTugbCzLjm35+P+T5HXRopRZKs891K2tMESEI+kKmY1XgBoLHSWoHHEiAmH0i5NqZZtdFmoqZqgxPdxYSVTGOrSCanGH0jql+7mwAJaSBXnjwLAHCy5v/pGzav2pVGk0g5mUyj/EuDY8lrVdmOaCDTnCCSDnUBoq9s9VIX33lsB56ybd8KXQOJTIAz4VWAJIueAwJ410DGJ6dAsMJwY6E0BTcHOhAdhTU8NhEuQGrKsOVwX9gxTZUueSDGMN7gzCzdQ2OoKStC99B4DA1kAvMbrMTF0+fVAUDY6tzNH6a6L/aPjoddfzbTbX+23YNjjiaZqZD0VCMCZBrhJiySjVWfWx/KaHZb2ZqcvvqkES8RbWxiKqlufl6I/HxU4pdXH0hk7w03VL5GLAES6RvRTVgAMK8h9Jmr8dWUWz/nYQ8+kK6hsbgRc17pGhxDQ2UJuofGXX0gfZoT3YSbAFGLhsHRSauMaw6gzHldg2PO92T6nWWDwz9ZRIBMI1KV2KRrIG4rW1PYqZ4HEk+AWCXMfQ7QI/rnc92q2bh19TLcuWGXZw1kxGPvdJWvMaumzPUYkwlLL8V+8Ukz8N9/3A8gJLTVij6qpa1hTAWEuBqfVzoHxtBYaSU3trv0KtHDeM37zfeN0kAGsjDk2xT1NT455XyuQfiZsiWyzA3xgQhJo2sgbj+aYwbzUyIayLG+ERQXps8h6ScKy4sDHQiZsGI5WKsjnehadWAAOG9RAwCrhHtkfxFdAxlx0UDqK/x3Royke2jMqWnm7gMZR00SGkj/qFnAZLJciCnqq2do3CmZH4QAyZbIMjdEAxGSZp6LD0THpIEc6BwM6ykRi6O9I1jYVImdx4Pophwf3d4+4jHE1Frtx1+TXXnKLHzxwdg92SJX65YTPbos+5Pb250QagXDWgkXFxbY/UmihVp9RYmjLeo91P3QOTiGekeAmJtdxdNA3KoSHOuzwr8HR83fQRCmoXh+wET8F91DY6irKEHHwFjg2f7ZiGgggm/UxKMnjpnyC4bHJo1moOLCAudH/8mLF8fUMI72DGPZzPQZwXUNZNRjBFC8XiCKyBwPE1FO9PFJlBdHr/fWbz0e5RsqLy50zChuTvT6ymJH2L96oMd6D59RQ10DY45WVejS7EqPwjLhpoEc6bFNWC4aiBt6+G884kXtJRLV1z04hnqlgSSg4eVq4UURIHlCJm7Aw3ZyoMoBAcwayLG+Ecw02PsjfSCtje41hI72jqRVgFT40EDiFVJMhMjJdngsPLRY8cKezqhukGXFhRgZiyNANBPWKwe6AQAHu4d8jbVrcNQxYc2oKTP6Qfp9OtEdH4iLBuLGYbtHzS5D/lEq6R4ac3xViZiwMlX3LVlEgOQJmbgBTW1sTU70Y70jmFkTnbcwPslOU6mOgVEs0Qotms6hyqVMpKERVZlfDcSDD8QLkT6QobEJowA5eXZNVFOmipJCJxt9dNzsRNdNWJsOWKXj3+r0J0A6B8ecyLLmqtKoSCxm9u1EV4UnE3WiH7X7rGzY1p7Q65Kle2jclwDJVcQHksVkax9kxSHDitUUHnrcRQOpKi1yVp4n+kdx6fIZru91tHcYLbWWr2VgdAJ19o+UmfHlh9/wfQ1uhPlAAnaieyE6jNcsnFavnIn/t25b2DZvJqyQBvLaWz0AgLe6/AkQ3YneXFMaFRAxMj6FwgJCSYwQZzcNpL6iBO39oxg0mMVioXxuT21vxycvXpzQawH/eRvdQ2Ooq7S+uyCTNbOVjAkQItoPoB/AJIAJZm6L2E8A7gRwNYAhAB9j5lfTPc5Mku39BEy5HYTo8NDjfSPGkNXqsiJn5XmiP7YGcrR3BLPrrHP0j4QEyPD4JB546aDr6/zip5iiVYk3GKVenWdsYgolRQUYdtFATAKkrMQSIBOTU5hiRlFBtG+pvqIY++3+LU3VpegfnfBtwuocCJltZlSHBIhbXo2JN4/0hgVUqAm7pbYM7f2jGEhQgKjora1H+xKqPBw5hjs37ErIOd89OOb4BKdDTaxMm7AuYeZVkcLD5p0Altp/NwG4K60jEwDE9q1Eho8C1oqxO+KHc6xvBLNqTQKk2Fl5dg2OYVGTuwDpHhpzoo301WrkewWF7kQfm5zy5Iz16kT3ggrxVRPnUEQYr2JBYyXOWdgQtq2iuBDDY1YL2dKiQmO4cENlCbrsz+6sBfUArLL8fujSTFgzqsucUF7drLqoqTKm9lxXUWI0war7JtHKzcp3ct6iBjyz80Sco0MkawbWTVjTQQPJtACJxXUAfsIWGwHUEVFLpgdlIlcjKLwQy7diygyvryyJij5xM2FVlxU5jvSa8mJnEjLRWFnqlAjR7eXdg2OOczaozGoATuVbAHY4bHwtZMTF35AM6lojKwTr/OoT54c9Ly+xBMjouLmQImBN2GplrgRIECasGdWlaDdoprEc6IC7D0SZLRM1YSkN5NLlM/Hk9vT5Qbq1kObuALP9s5VM+kAYwONExAD+i5nXRuyfA0C3TRyytx3VDyKim2BpKJg/f37qRhuDbDc1mTCZCxLlhFEDKY5yHh7rNWsgNWUhH0hzVamTgGWipS70et2c0TM0jlNm1+KFvZ1RtaGSQZ+sS4sK7D4fsYVDkFFYCvX5RNbnioXygYxMuJdWqa8odlbIbZoA8TPhlRYVOnkpM2pKjcmEsRzogLuGoe6bRE1YygdyyfJm/Otj2xN6bTJYUVjWtRYXFuR9QcVMaiAXMPOZsExVnyaii/ychJnXMnMbM7c1NzcHO8IkyHatJIiIrU6TBmIwYR3vG8XMapMACZmwmqtLnUnb1BuiRRNA4SasMdTbTsv1W495HrtyjE+6tGDVBUhZcaGnfhVW97+gNZAJJ6chlhNaR9dA3MZTX1Hi+EAWN1umw4qSIpzwUW9MaR9AuAlLJ74GYhYQLT4FiNJAWmrLMbu2PM7RwWFlolufR31FtDaeTtIx52RMA2Hmw/b/diJ6GMA5AJ7VDjkMYJ72fK69LSfIRa0kkagvZjaasOoqiqOclu39I5hhCOPVnejN1aWOrb53eBwzqsMnvlk1oUlAN3f02Jm/APD41uMxr0/nz3bk0c7j/VhhKKuur/bLigs8lRMfGZ8MzAeiGBidsJtJeT+v0kDcIrAAy9SoVscqj2deQ4XRDxIvIqk+TID4M2EB5mg3PfIuEfR78NLlM7D1aHDaaSy6NHNefWVxxvwgh7qH8B9PWt+Z3wRRL2REAyGiSiKqVo8BXAFgS8RhjwD4CFmcB6CXmY9CSBmJ5JIMudyUDZUlUdnoVaVFxpVwdYQGojD1hpitmbD0TGfLaWlpIG51mExs2m/lPmyyk+gi0Sfs0iJvGkhqTFjjdhZ6AgKkRAkQcy8QAKgsKURJRNn5+Q0VYX6QR9+wfm6R98WnL1kS9rpGTYDUVRQbhW08E1a1Zs7UURpIpA9kn609uaF6rgDANadZrtN0aAP9WktnSxvPjADZsK0d7z3D6jf/hQc3p+x9MmXCmgngeSLaDOAlAOuY+Q9EdDMR3Wwf8yiAvQB2A/ghgE9lZqi5yT/+75spPX/HgLn3dZ3hR2NyoAMRGoh2rl6DAJkVy4RlayBueSQv28JC5yV726tuAqTYjwYSvAlrYHTCNYnQjfJiK5EwlgZCRFE+p/kN5XirM5Qc+sBLbxlfe8l3ng57rpuwiChsMaCIp4Ho9wIATNmmRaW5RiYSKuHmhn6/KA3zn363NeZrgqCqtAiFtkZnlbhPvQAxaY1PbDuO1Sut38NhQ8JvUGREgDDzXmY+3f47mZm/aW+/m5nvth8zM3+amRcz86nMvCkTYw2adPhGmBnrEzDn+MFNgJic6CYHOuCugZji58N9IONhxyoT1mqtRapOZJ4EEEqeU1nYkVRkjQYyYRdS9G5tLi8pxMh47CgsAPi361eFPZ9XX+HkggyPTeLFfdGfzeDohFMmRKFrIABcBEg8DaQ4bGHQYWegKw1qcGwizMEfzzTUYrjnXj7Qhad2pDYiSxemVrmY1OeC/MpuIa1/Pn9+qwcXLrV8wndGfM9Bks1hvHlJOkqOHOgcMpYZCZKOgTE0GcJuTU50t74XutlCn3RMGkiL5ggdiNJArMnpIvsHo9u/OwdGseOYZf/WHeYztfwC3WavjtFX7p41kBQ40b/92A5c+x/PY9vRPs+LjXI7DySWCQsA3ra4Key5bsJ6cV9nVDkVADjSE31fNXgSIF40kND3Gln+v6y4MMxs2q31HTdhWrT8y3tPw1cejrSUB4uu1dVrodKppMOOhtQF/pkL6lFpf3+z66zfzpRLwEgyiADJQ57b3YFFTe6FCYPAVQOpjP7RzIghQPo0J7pCFyCqVpa+PzKRUGkgym+hrzKf29WBC5ZYgkU3ZZ3daoWunjW/PswPorLo9eS70qLCKAevSZO0OhIGK0A+dK4Vmn7BkibPiw0vTnQTuhP9uV0deN9Zlg1dn6QP9ww7n12vvVCojxAgMwwCJFYvEGt/cZhmGdm+uLK0KMyR3mkLELeoMdOi5YKlTXjb4saY40gWvelXQxJO9ESCBlQ05N3P7HG2Xb4iZM4ttn1dA0mW7TchAiQPeW7nCXzg7HnxD0yCzoExY+Kfnl+gcNdAisPyQBS6AFFlMYo1h68uQHo0DUTx+Jsh893TO9px8UmWAPm9ZjdvW2Blb5+5oN6pRguYI1bKikOJhErgmTTJVOaBJBSFVaJrIN7H01Jb5lQXeHbnCVx9quV87tM+7yM9I07Y79M7LUEdacLS2+0qEnWiR2og1RECRJlJ950wO9PdzKafvdyKIHPrP5IsugAx+QO9cttD3uu7KWG69UgfthzuBWD2B/amoGqDCJA85IW9nXjvGXMAhPoiBE0sJ3qkD2NWbfRxgLUqVdVWVeQKEC5AVFVVnbAorMGxsB8tADy/q8PRGJ7d1RESIFuOOWq8Kv/RtiBCAzEIEF0DURVrTROQl2TDRFGTZqJO9OHxSbsboffXFRUWOAmbHQOjOHVOLYBws9XhniHHJKL8bJEmrCXN0SVp4puwisM+U6MGogkYtUhxi8Zqccn9UN0zH3jRHCDghTufcK+JpS9mTBGJXtmeQNix0kD+9u0LcesvXwMQ3iVUYTINJ4sIkDxkbn2FE/lk6gQYBG4CpKasKMrc4x6FVWzsKaLf6KqhkI4ydUxMTmFwbBI15eGr2xUtNfjTng4A1o9Y/ZjqK0rw8J+tVCLVx/30eXXYeazfea2pT7iugSgHsl76XJmzXtrfhevXbgw0MEJda0ICxEkkTFwjmm9rD29f0uREE+lC/EjPiCNAnrVrTOkhswCw2FAUM14TLcucqWsg4QuHqtKisFDe7rgCxL3fPADc88d9ngIjTNz/4gHXffURTnQ/Gsj45FTcMGUdpTV+6Lz5xtbRClN4fLJIOfccJDKxS6ESvC5cGnKOHu4ZDjMp/O71I1GvM63y49HhYsKywkNL0DEw6qz2YznRTX69h149jIdetSb68xdF26yVqaN32GpSVBhRbfaKk2c6ZqyLl4WqE7zz1Fn49mM7nHEClnN22axqbD7YA8AsQHQNRK3G93UM4hR7ha4njf7hlguxfFZNwm1V3VDXmohmU15SiCHlA0lQI1LCVgUkAMBhTYgf7hl2cnIWNVfhtYM9aIi4D1SLY728vRcnuh7dZdJAlOY5Mm4VigSAvRET7WRE+K8by2ZW47evRf8W4tHePxIz3yjcB+Ledz5WyZgDnYOoKbdMwcxsLIapGLFL1gCWkL7jr1bhxp+YA1Z7RAPJf7yE+Uba39V/5WS9YElIgOjmh86BUXzm5382vmeiY3tpXxcefd0ci/+v7z8VAHDVnVZhgUgTk6KipDBq8gdCDu79t19jtGUPjFohnXrlU53VK2fiCbuR0MUnhWzBV5/aYixBr2pBAeb+4KWaBqI+zwOd5hVi0JnoSoAkasIaiZMH4obSQC5cFrqHjuomrO5hzK2zjlGO2sqIsamil3tPDDoTZTwfiF7WhpmxO6KTYHVZyISl5/5ErtSVOSee6e4TFy3GD5/dG/MYE28e6TNq3grdhFVnCGkHLK3ybbc/6XqO3e0DOHN+HQDgSAyNArArIWsa4OUuoeyAmLCymp9tdFdr4wkFff+dG3Y5mkSiYb7KhqyX99ZXcjuPR7f3HBmfxGP2al13skWOGQgXXLXlxfjcFScZx3Hpcusm/so1KwGEm6d0iAhVhlBR/UbfeyJ6zMWFhOHxSbuMSfTEtKCx0nHsnr0wJByWzqjCOa0NUccrAcLMxnIa4RqI9Xnud+nelzofSAJ5IMV6JnpiP/EFjZZw0H0I+iKkvX8EM22f1pUnzwIA1xXynhMDTvhzvDpeeiLhkd6RqHNWlRZhcCxU+l/5XQ52DYWFZ5sWCCbevqQxLDDDK1uP9OFdp1vBBaZIKd2EVVZcaHyPF/Z0RmlYOruOD2DJDKt9czxfiF5KPx4iQDKMmyCYnGJ805CspoiX+xFrvzKteGHjnk4A4ZPY77ccdcb7wR9udMwLKjx2w7Z2x1n67K5Q3wSTlqMzODqBuvLYq8qLNPORGybThrrRmRl7DFE2VaXFGBiZcNVAAMuMBYSvRIkIv7r5/KhjL7EjVn668YCxRMum/V24c8MutK5Zh3V2JJerBhJwFJaapBIpZaJa2o4m6EQHzNE7+iq4vqLEOefSOD3qd7cPuJZpj0SPyHvjUC9Om1sbtr9S616pC5CmqtKwTOtYE7MOETklThLhzSO9OH1uHQA4EU86kfej6f58bldHzPfY1T6ApbYfabvmnzPRMTAaFcTghgiQDOBFOzjQOWi0nQfRC+Cup/fEP8jmhb2dUduaqkodAfChc+fj7y5YCMCKVAKAh/98GO+xI7aeSqBvQkNliatmkQgm04aK4jrRP2pcudbYDtduFw0EAD5z6RLjdhNK4N75xC5sNHyGl62YiY+9rRVAKB9lX0d6NJDiQuszTiSMt0xlok/EzkQ3vtYwft1HNqfee2XbPScGwhzjsdDDeN843OP4l/T9yonepfXcWNRcib0dIS01lhM5Et3U65U3j/Th5NlWaZQ3DpkESPj9qCpF6zy76wRaG6OjpBS72wewdKY3AdI5MBbTpKbT42JhKJm15CxPJzAgAiQOXjLHtx7tM64Qk21kMzXF2LgvekJzY+Pe6NITuvlh1/EBnGSvGlU00ot7O3GlvVp/eucJ1/LmkTR6vGnjYdJApmzBu+fEoDEhsso2d/QMRYfwKvwk9P3zX5yKn22MDu/Uo7CUmW9gdDzMhDFlyGAPAiVg/YXxJm7CMnG8NxQQoSKwvLDnxGACGkjIhPXG4T6cFiFAKksKnc+7e3DMMVEubKoM84N41UAAOEIqEaHT3jeKRXaY8uZDPc52tVis86CBDI9NxtTO93YMOLk2qoqCjj75f/7Xm7HfRRuORI/C0ue1sWO7X/F0AgMiQAJg29E+Rx3WTVz/8eRuAIndoGHnPdaHogLrK4qnzfQMjRmLqh3VImh2tvc7ZgelbVy8fIYzSTVXleI1jyYzUxkTP3zOIJBVTsjejgEsao4WINVlVlKZXok3CK48eRZuuTy6uVZpUSFGbQ1TRfcsaKgMM2MpARMrYsYPqpxIIias4sICFJBl/goiM76mvMgp3T8nAQGyr2PAs9kk3ITVg1MjTFhVZcWOAOnScn8iBUgiYesqgOP53bFNSjonzap2XveGZsJSEWKRGrPJvHTh0mbXPBXACotWZUgOdA454cbbbWHy95ctdSb/T1y0CFesnOVp7GLCylK2He13Cvm9+OXLAFjaikri86uJvLCnE1edYp33eF/sUuUv7uvCmVo0kWKK2XGuE0ITv3K0v/eM2c6xl66Y4dmM5VVtjsd5hjBdldex98Sgs9rTqS61Jhu9F0hQ3HJ5tEDTNRC1Al/QWIEDmiPd5HwPAqWhJWLCumP9ToxPMn6/5Rhu+eVrSeeltNSWOyG2s+PkV+jUV5RgRxwTjEI3YRUXFkTlDlWVFjpRWHrPjWgNJPGQ9Oc13188TpkT6h3TOTDmlO0x9UABzEmVFy1rcpJrTT7VJVoezbyGCuxpt65PBerogQJu4fQmUiFAJA9Ewy2/wsuxKvxR2UV7h8adFeqGbcfx1+eG2u2++pa5hHgkf9rTifedORc/2/gWbv/9NvxGi1uPbOqzcW8nzlvU4CR3KWbXlTtmrKUzq50V8nvOmIMN29udip2A5UD92m+jy8Afsiu09g6No9Ze8QelgZhwNJATAzh3YXTU1IGuQXzq/lftZwedciepQo/CUivwyIlrxGdSWjyqbAGSiAnr1tXL8MBLb6G9fxT//bGzccnyGUnlpcyuK3NMQ3MMGc5uLJlR5TTuikdZcaG1wgGcoA6dqtJizYQ17giQRU1V2GsHWhzrHfEssHSe393p2V958uzQ2FbOrnG0EGVtiOTGCxfhu+t34omtx/EOuyLCBUuanIjI/bdfE9Z07s4NuxwHOgAsn1WNHcctzeMR+/d/oHPQuQ87B0c9/xZ7hoMv7CgaiEakvyNWFJJ+bHVZEb72Litk9XX7hnp5fxdW2bHcL+7rckpkDI5O4IZ7X/Y0npf3d+G8RdYEunJ2TUxfzMa9XcbVfEtduWPGWjYzdGNeY9c50sMMz5hXhyOGFdzTOyyhpGofAcFpICZCJiyzBnLeokZ85ZoVAICf33huSioa64RrINbKeEFjuAnLS7VeP1SXWp9FeXFiaz0lcILwgbTUhhYhemOveCy2Ew29ogouRpqvAEuQ6k50JUDm1Jc75rUb7nsZH79okef3c85dWhjXWa1QDnQAOG1OLV63F4yvxGlM9pXfbMEz9u+osarUtVYXAMeBDlgCRI3tAjtBWK+C0DkwFlUJwA2phZWlrJhV46zs37Aday/u68R5C60J/eTZodIaP914wJmM4jGnrtxxVse7wQ92DRlXbnPqyhzzwzIt7NIUQVVUWBCWua1QAkQl5wHhTnTl1AMQt8eJl2NVePDR3hFHs9PR7eVBm7BMlBWHNBBlwmptrAjLBUmVCcuPBgKEoqkSjcIyMaeu3Ml/ScQHsri50nNeBhAKGDDdx//zyiFsPtSL1jXr8MLeTqyzk1gLC8iptLCypQaffMdiz++nuGBpE56LY8ZS36/+GzptXp1T4fkb150c8/WXr5zh1KkCQtUZTJrPkjANpAbbj1q//b85bwGA8BykzoFRzyaswbHJwEu6T3sTVrx+z15Yqa1KlEr74r4uJ5Hu8hUzscH2LfzouX348tUr8JXfxO9LcL5Wenrn8dgC5MwF9cakpdna6nFZnLh9ALhu1ZwwUxlgRWoBVu0jlT+iq82qlIeOm8nEy7FKA5lbV24M460uLXJMKqYwyaApLSpwTFSOD6QpUgPJHh+IfnwQTvSWujI8Y5tGa+Pk/uiYamLFQl2rSYB8+pIleHzrMRzvG0VzdSk+e/lS/NzumLioqRK72wfwzfee6iuI4YIlzTHrWwGh358e5nzanFr8/QNWZQeVPOvGmneuwJ92dzrOdvX96L1zlDBZ0hz6nZ40qxp/tJ38qqzPW12h+65j0LsGUlFSaLXcDTDwZNprIEE0eFrRok/M1g28u30Ap8+zfgiXrpiBJ+3V+zkL63HdKstxHS9kVm/4s7t9IObxytQVSUtdYgJE1dHSM8DVRNDaVImX7aY16TBhmSKwgPCQT7cw3mTRNaV3f/+PUSvwlpqysLj61JmwfAqQ4uBMWLPryvFn22+XyARtciDHQgkQU/+YSs2JHlmB4P+95xQA8bPd3Th/caNra2PFm0eiw2kXNFbgqpO9RUBVlRbhic+9I2q7HqGpFkX6BD+3vtwRWuqz14M3igvI871RW14cuB9k2gsQL6iV5vVrX8Dl33smav/KltCKSWXQnja31ln9LW6ucr7kz1yy1FHV3SpuqrA9vSRJc3VpzHhvk/8DsGzWr9qOTC8Zq6qOkcoTAeCUQ798+Qys32aVPUmlAKlxBIh5AqoqK3IK2gWduKdQC4v9t1+DZ754MSpLrfdRVV4LCijMvJYqJ7qTB5LgdVYEqIHMri3HoCFDPx7N1aXGroZuKH+PicqSIidZt7SoMOy63BqWeaW2vNgJcjFFKr12sMdYwp2IcPeHvefgmczGx/q0PK326LI9RITv/OXpYdve6hxytJVE8rFqy4sDj8QSAeKBj/235fT+zCVLcccHVgGwQmzVRK87vZT6fe7C8Aldhfnq5q7vPLYjqt7UF369GZ/6mRVhpJsLTppZ45QdZ2bnhv7yw2+EvW8ks2vLnTalifDQq4cde+k7bL/I5StnuvaACBJHA3HpqlhdVmzMeUkVZcWFThi1njm/oDE0vtEUm7CKEqzbFKQPpLm6FEU+qg4QEZa3xNd6FTXl7sKmoIAcrSoV996Xr7aCMj784xej9t1w78txfRx+OdYbiiA0JQ0CwFWnhGs5RCHTl1f/B2AVdxQBkkKYGQ++cgiAlVms7NqqaNwFS5ucCJHbHnrdqYUTZhe1958bYVL6/BXRprGFzZWO+eyHH2kDAJzT2oB3nT476tiTZlU5jvTfvX4Uv99iORFXtFgCya0wXKxoj1hUlxU5Wshpdu2f5bOqoXx+fs0FXqiNo4FUlxWlvOe7TmlRgbExl16OIlUmrKo4ZdDdCNKEVVhArj1d4nHPx872fOynL4ldfkZ9FpEtdINAmYfOtottXvqdp3HWP60HADxw03m4wqOpKlH0viebD0aXRjHRqvnfvPo/gDzRQIhoHhE9RURbiehNIvqs4ZiLiaiXiF6z/76ajrF9/ZE3nRLPl3z3aXzkxy8BAL50ZXTV2VPn1uGLv37dsN0SIGfOD0/qM5kStmp21fVbjwEAPnD2PKc2lc5Js2qcGPd/fWw7vmqHDX/Yjsxwo6y40Je56X1nzsU3frcVQChjl4gcTSqVxPWBlBY5/SDSgZuZrFXTkNRiI5GINC/EK4PuhjJhBWXiSyR8VyeR8esanQmVnd0QoBM4EhUevvYjZ+H3t1wIwJvv0C96lJrXkOf5DaEk1kTysWrLi6O6hSZLJjSQCQCfZ+aVAM4D8GkiWmk47jlmXmX/fSOVA1LO6S1H+vDrT1rVWu+8/gxcvtKqTGqyXX7tXSvRZaj1P6Pa+qF5+eFu00o1x8tWXz6r2okEWTajOszBHg89/8Mr162abaxjdK2PCqaJTqrKiRjZa1vhd1L1S4mLdtdm9y15eke7I0B038lnL1vqRJip608UU7l7L5TZAsRt7InyeZfS/elE+VMaElh1J4rSRJbMqHZ+y6lEOc7b+0eM5eFN6FUQEjHn1eSDBsLMR5n5VftxP4BtAKKX3GnkO49bXep+esM5TuvNVfPqcNNF7jHlTVWluP/GcwFET5BeGRmfdLKo46miC5sqnXyONe9c7vk9AODnHz8voeMByzH57fefHrW9zdBPIx76pOol0m2p3QvBLeInXne7oCkoIKPJbvksy3z4hV+/jkc2R3e3i7xuPyxJMBRWoUxYQVRMBtyDNNKJMmE1pCF0O12o2l2vvdWD0+fVeXrNgoZKHLBDeRNxoteVlwTe1jajPhAiagVwBoBozxVwPhFtJqLfE5GrB4uIbiKiTUS06cQJ7zVtFCPjk/iFHU+eSNMeIFTWwO9EsXJ2jaOFKG3HjeLCAiy0TSbx+jAEZUZ531lzU/4epvMt/vKjMfef/LXHEMy06J1YvoRffeI8dAy4h0cm81klknehk2jiYS5QWZI6H0imUBrI5kM9WOVRgMxvrHCy0RM1YQWtgWQskZCIqgD8D4BbmDky/OBVAAuYeYCIrgbwGwDGDD9mXgtgLQC0tbV5SrOMTB5c0FgRltCTLla0aAJkxUz84KnYvT+++d5T8b67/hT3vIkk9vnF9B5Bni9Sk4vcf9rXH/PcayIIyooLnez3SBY1V+EPt1yIlV99zLg/Hd9HJIlU780VlAbiZtrMRVRwxmsHe3DjBd7KsCxorMABOwoxUSd6PvhAQETFsITH/cz8UOR+Zu5j5gH78aMAioko8e4vLujJg2e31jshfEHiZdW5sqUGj26xnOeqy1kszjJU252upNsPcpFWdNJEotprqklVfkwmUf6gVCWPZgKVV/T6wV7PJqyZ1WWOfzLTYbxpv+vJMmz/GMA2Zv6eyzGzABxnZiaic2AJOu+dlRLgQOeQsY1nsnhZda6cXeO0rA3KVj1dSLcf5LsfOB3/8+qhtL5nMmSbQAuCKseJnj8CZGZNGfZ3DqG+ssTzdRUUEObVV2BX+0BCAuR3m4/ghb2dzsLWT9mmqLEkfYbEeTuADwO4VAvTvZqIbiaim+1j3g9gCxFtBvDvAK7nIPrDGvhA2zzXHIpUs2RGVcwEraD9DPnEZSuCF/r5RHmJdU/n0z1UmYcCRGkgXv0figV2DlJDAtrYpy5Z4pTi8Vu2KZK0L1OY+Xkgtg+Umb8P4PupHIcKu/yrs+c523TfSOuadY6ENkVY6fsTQX+Pk77yB1yxciYet7O7I0nEdm4ae+RrI4/JZb545fK4PqN8wu3edEOZsDLhf0kV1WX5J0Bm+hQg8xusgJpEKhSkIow3//RcFyId56ps+TytnpHbjy1WZFWiP0bTewQxmXtxapsc1YlOTNlGPgnFWCQqCBKdkHIBZcKqSbP/K5W02Bn+qneQVxY0Rrc5iEd1aaieGOC/1bbOtBEg6gfYumYdHr/1Ily/diMA95W7lxV9PpDrK9R40VvZjJvw8yvIdU051xYCXlAmrHzyF6pSQytbauIcGc5SH/lBBQWE6rIiJxLrzg3JmzOnjQDR+czPX8Wady7Hlx583XXlHmSY6nQRRsD00QiCwO0e8/sZRmrKubQQ8EIilX1zhSV20myiUXN6r6BEUKG8U1OMDdtiV7/wQv59Ix5Y0VKDvzxrLr70YHQtK6+4mX5MgsKLMErn6jGVAk1da7YLzWxeradKq8p1c+WyWamrSZUp/FYa8NM4Cwglpm450hvWFsEv00qA/PY1q7qs385lOrE0Fz+o1WM6Jt5EtSs/k20QGlwin4Xb5Oi2PVOr9UwK1lw3V6ayB810QQmQJ7a14/IVM7H2xN6kzpf3AiTSeQ74L1CXDrxMvCYTRxCTUbZNtokIIbfJMdsmzaAz+AUhEZQA2bDtOL567UqsfVYESExuXb0MN164EKd+/XF87wOn43O/2pz15pV4pGoSSvdkm+vfQ6bIdVOUkDmUADncMxxIZYu8FyDM7Pg6/uLMufjcrzbn5CrQ72SbzZN0qr6HVEbWZYPvJNu0KiF3UALk4mXNCXe5NJH3AuSpHe3YYffRCNLck278Tra5KCyTJZHIukRX8/ke6STkN3V2r53LA2oMl9cCZGJyCv/86HZ8+Z0rcONPNvnuyZAqslk7SCTKLJfJ9dV8NmhVYlLLHZQGctGy2MVBvZLXAuQXLx9Ec1Vp1tZNymbtIOgos1wgmwW6G0HcQ8lqVV6EsOQHZQeqt31Q2fx5KUAiI68W3mY1KMqVSWE6k21hrkIw5GMZnVzkHQFpHoq8FCC3rl6Giakp/OCpPVlnthJiI5N4eognqFM1uevnvXPDrqj3zSXzYS6SbP5bJHkpQI72DuP+F9/K9DAEIWuJJ6jdzFLJmqLSkcQqpI+8FCD/tn4Xrj97Pu5+ZvqU+55OZIOvIh/NL14+13QXr5Sot+wm7wTIruP9eGTzEadssfg98o90m7ncJlY1hmwQaEGQyOeaSK6NkL/knQD51h+24/NXLMONF3prUC8I8fBj7gmCbBZMieTaBK01eBFSYvpKD3klQN443IuObe14Yls7+kcmsubHJgh+mK4BBfEEpxchFas4qRAceSVA5tWXYwLRdlNBmE4Erbmku4xOkIIznqkxaKabwMorAXKwexgtyD51XxDSSdCaSzrK6AQlpBTxNJdUTfS5XtkgUfJKgJw6pxabRPsQhJwj3bXeYgmTWEw3DSMeyZdj9AERXUVEO4hoNxGtMewvJaJf2vtfJKLWDAxTEIQs4Y71O8Mc43esT76ft+LW1cuw//ZrwoRB65p1Ye+nHqvkR3W8qS+PGpvbmCO3p/NagybtGggRFQL4AYDVAA4BeJmIHmHmrdphNwDoZuYlRHQ9gG8B+Kt0j1UQhOwgHQEFib5HItF5sTLwI8u66PtjNY+LPG8mNCJi5vS+IdH5AL7OzFfaz28DAGb+F+2Yx+xjXiCiIgDHADRznMG2tbXxpk2bUjd4QRCEgIis2efH76O/ztR9NRKlNalAo9Y163DgW9e+wsxtfq4hEwLk/QCuYuYb7ecfBnAuM39GO2aLfcwh+/ke+5gOw/luAnATAMyfP/+sAwcOpOEqBEEQcgs3AXP0vlswenSXryJZOe9EZ+a1ANYClgaS4eEIgiBkJW4mN/rWta/4PWcmnOiHAczTns+1txmPsU1YtQA60zI6QRAEwROZECAvA1hKRAuJqATA9QAeiTjmEQAftR+/H8CT8fwfgiAIQnpJuwmLmSeI6DMAHgNQCOAeZn6TiL4BYBMzPwLgxwB+SkS7AXTBEjKCIAhCFpERHwgzPwrg0YhtX9UejwD4y3SPSxAEQfBORhIJBUEQhNxHBIggCILgCxEggiAIgi9EgAiCIAi+SHsmeiohohMA8jEVvQlAVBZ+HiHXl7vk87UB0+P6Kpm52c+L80qA5CtEtMlvrZpcQK4vd8nnawPk+uIhJixBEATBFyJABEEQBF+IAMkN1mZ6AClGri93yedrA+T6YiI+EEEQBMEXooEIgiAIvhABIgiCIPhCBEgWQUTziOgpItpKRG8S0Wft7Q1EtJ6Idtn/6zM91mQgokIi+jMR/c5+vpCIXiSi3UT0S7vMf05CRHVE9CARbSeibUR0fj59f0R0q31vbiGiB4ioLJe/PyK6h4ja7S6oapvx+yKLf7ev83UiOjNzI/eGy/V9274/Xyeih4moTtt3m319O4joynjnFwGSXUwA+DwzrwRwHoBPE9FKAGsAbGDmpQA22M9zmc8C2KY9/xaAO5h5CYBuADdkZFTBcCeAPzDzcgCnw7rOvPj+iGgOgL8H0MbMp8Bqx3A9cvv7uxfAVRHb3L6vdwJYav/dBOCuNI0xGe5F9PWtB3AKM58GYCeA2wDAnmuuB3Cy/Zr/JKLCWCcXAZJFMPNRZn7VftwPa/KZA+A6APfZh90H4D0ZGWAAENFcANcA+JH9nABcCuBB+5CcvT4iqgVwEax+NmDmMWbuQR59f7BaQJTbnUIrABxFDn9/zPwsrJ5DOm7f13UAfsIWGwHUEVFLWgbqE9P1MfPjzDxhP90IqyssYF3fL5h5lJn3AdgN4JxY5xcBkqUQUSuAMwC8CGAmMx+1dx0DMDNT4wqAfwPwJQBT9vNGAD3aDX0IltDMRRYCOAHgv20T3Y+IqBJ58v0x82EA3wHwFizB0QvgFeTP96dw+77mADioHZcP1/p3AH5vP074+kSAZCFEVAXgfwDcwsx9+j67tW9Oxl4T0bUA2pn5lUyPJUUUATgTwF3MfAaAQUSYq3L8+6uHtUpdCGA2gEpEm0fyilz+vuJBRP8Ay2x+v99ziADJMoioGJbwuJ+ZH7I3H1eqsv2/PVPjS5K3A3g3Ee0H8AtYpo87YZkCVHfMuQAOZ2Z4SXMIwCFmftF+/iAsgZIv39/lAPYx8wlmHgfwEKzvNF++P4Xb93UYwDztuJy9ViL6GIBrAXyIQ8mACV+fCJAswvYH/BjANmb+nrbrEQAftR9/FMBv0z22IGDm25h5LjO3wnLWPcnMHwLwFID324fl8vUdA3CQiE6yN10GYCvy5PuDZbo6j4gq7HtVXV9efH8abt/XIwA+YkdjnQegVzN15QxEdBUsM/K7mXlI2/UIgOuJqJSIFsIKFngp5smYWf6y5A/ABbDU5dcBvGb/XQ3LT7ABwC4ATwBoyPRYA7jWiwH8zn68yL5RdwP4NYDSTI8vietaBWCT/R3+BkB9Pn1/AP4RwHYAWwD8FEBpLn9/AB6A5c8Zh6VB3uD2fQEgAD8AsAfAG7Ci0TJ+DT6ubzcsX4eaY+7Wjv8H+/p2AHhnvPNLKRNBEATBF2LCEgRBEHwhAkQQBEHwhQgQQRAEwRciQARBEARfiAARBEEQfCECRBDiQEStejVTbfuP7AJ0iZ6vjYj+PZjRCULmkDBeQYiDXZfsd2xVoBUEwUY0EEHwRhER3W/3+HjQzsZ+mojaAICIBojom0S0mYg2EtFMe/tf2r0zNhPRs/a2i7VeKI8S0Wv2Xy8RfZSsfinfJqKX7Z4Nn8jcZQuCOyJABMEbJwH4T2ZeAaAPwKci9lcC2MjMpwN4FsDH7e1fBXClvf3dkSdl5quZeRWsDOEDsLLXb4BVJuNsAGcD+LhdWkIQsgoRIILgjYPM/Ef78c9glZ3RGQPwO/vxKwBa7cd/BHAvEX0cVgOmKIioCVZZkL9m5l4AV8CqufQarHL+jbDqEglCVlEU/xBBEBBd0jvy+TiHHIqTsH9bzHwzEZ0Lq4nWK0R0lv4iu+PbLwB8g5mVo54A/B9mfizICxCEoBENRBC8MZ+Izrcf/zWA5728iIgWM/OLzPxVWM2m5kUccjuA15n5F9q2xwB80i7tDyJaZjemEoSsQgSIIHhjB6we9dtgVdj12g/720T0hh0G/CcAmyP2fwHAFZoj/d2w2v1uBfCq/br/glgLhCxEwngFQRAEX4gGIgiCIPhCBIggCILgCxEggiAIgi9EgAiCIAi+EAEiCIIg+EIEiCAIguALESCCIAiCL/4/bSNbt3Kdl9wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result from the automatic windowing procedure for comparison:\n",
      "Result\t 3.27194697e-01 +/- 2.02097394e+00 +/- 3.22723658e-01 (617.667%)\n",
      "  t_int\t 5.70449936e+00 +/- 1.53928442e+00 S = 1.50\n",
      "Result\t 3.27194697e-01 +/- 1.96872835e+00 +/- 3.38140198e-01 (601.699%)\n",
      "  t_int\t 5.41336983e+00 +/- 1.59801329e+00 S = 2.00\n",
      "Result\t 3.27194697e-01 +/- 1.89700786e+00 +/- 3.67353992e-01 (579.780%)\n",
      "  t_int\t 5.02613753e+00 +/- 1.69573607e+00 S = 3.00\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeXUlEQVR4nO3deXwddb3/8dcnS5M0Sdckha6pXVC42uIvtJQWoewiiterCBev/BRBpV4RXG5xA72K6GVz+ymVi/DzsogLmwUEylahLF2hLG1ZWtrQNumevVk+94+ZhtM0JyTpWXLOvJ+Px3mcM98znfkMnLzne74zZ8bcHRERiY6cdBcgIiKppeAXEYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfBL5JnZuWb2ULrrEEkVBb9kFDOrj3l0mFlTzPS5/Vmmu9/q7qf0cv3/18z+0Z/1xCzjfDN71czqzGyrmd1vZqV9XMY3zWx1uIw3zeybB1OTREteugsQ6Qt3L9n32szWA19w90fSV1HfmNlxwJXAae6+wsxGAB/tz6KAzwIvAJOAh8xso7vfkbhqJVupxy9ZwcxmmNkSM9tlZpvN7FdmNih8r9LM3MzyYuZ/3My+EL7erxcfzvslM1sXLu/XFngf8FtgVvgNY1c/Sj0KWOLuKwDcfYe73+LudX1ZiLv/zN2Xu3ubu68B7gFm96MeiSAFv2SLduASoAyYBZwIXHQQyzuDIKQ/AJwFnOrurwBfIgjuEncf1o/lPgucamY/MLPZZlYQ+6aZzQ93Nt0+ulugmRlwLPBSP+qRCFLwS1Zw92Xu/kzYA14P3AAcdxCLvMrdd7n7W8BjwPQElIm7LwY+AXwQWAhsN7NrzSw3fP8qdx8W7xFnsVcQ/C3/PhE1SvbTGL9kBTObClwLVAGDCT7byw5ikVtiXjcCJfFm7Ct3fwB4wMxygLnAn4A1BDurPjGzrxCM9R/r7i2JqlGym3r8ki1+A7wKTHH3IcC3CQ6AAjSEz4Nj5j+kn+tJ2OVs3b3D3RcBjwL/BGBm3+5y5tJ+j9h/b2afB+YDJ7r7pkTVJdlPwS/ZohTYA9Sb2XuBL+97w91rgWrgM2aWGwbmpH6uZyswdt+B4+6Y2RVm9nic9840s7PNbHh4wHgGwZDUM2GtV4bHD7p9xCznXIKzg0529zf6uS0SUQp+yRbfAP4VqAN+B/yxy/sXAN8EtgNHAE/3cz2PEhxE3WJm2+LMMw54Ks57O8Na1hHsqP4H+C93v7WPdfwIGAk8H/ON4Ld9XIZElOlGLCKJZWYrCYZftqe7FpHuKPhFRCJGQz0iIhGj4BcRiRgFv4hIxGTED7jKysq8srIy3WWIxNfYGDwPHtzzfCIptGzZsm3uXt61PSOCv7KykqVLl6a7DJH4Vq4MnqdPT2cVIvsxsw3dtWuoR0QkYhT8IiIRo+AXEYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRiFPwiIhGj4BcRiRgFv4hIxCj4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hETNKC38xuMrMaM1sd03aFmVWb2crwcXqy1i8iIt1LZo//ZuC0btqvc/fp4eP+JK5fRES6kbTgd/cngR3JWr6IiPRPOsb4v2JmL4RDQcPjzWRmF5rZUjNbWltbm8r6RESyWqqD/zfAJGA6sBm4Jt6M7r7A3avcvaq8vDxF5YmIZL+UBr+7b3X3dnfvAH4HzEjl+kVEJMXBb2aHxkz+M7A63rwiIpIceclasJndDhwPlJnZJuBy4Hgzmw44sB74YrLWLyIi3Uta8Lv7Od00/3ey1iciIr2jX+6KiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRiFPwiIhGj4BcRiRgFv4hIxCj4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRiFPwiIhGj4BcRiRgFv4hIxCQt+M3sJjOrMbPVMW0jzOxhM1sXPg9P1vpFRKR7yezx3wyc1qVtPrDI3acAi8JpERFJoaQFv7s/Cezo0nwmcEv4+hbg48lav4iIdC/VY/yj3H1z+HoLMCrejGZ2oZktNbOltbW1qalORCQC0nZw190d8B7eX+DuVe5eVV5ensLKRESyW6qDf6uZHQoQPtekeP0iIpGX6uC/FzgvfH0ecE+K1y8iEnnJPJ3zdmAJcJiZbTKz84GrgJPNbB1wUjgtIiIplJesBbv7OXHeOjFZ6xQRkXf3rj1+M5vdmzYREckMvRnq+WUv20REJAPEHeoxs1nAMUC5mV0a89YQIDfZhYmISO9d9/Bafr5oXef0xSdOiTtvTz3+QUAJwc6hNOaxB/hkIgoVEZHEuOTkqay/6iMArL/qI1xy8tS488bt8bv7E8ATZnazu29IeJUiItJn3fXsewr57vTmrJ4CM1sAVMbO7+4n9GlNIiLSa/ECft+jcv7Czh5+X/Um+P8E/Ba4EWjv11pERCSunkL+YAI+nt4Ef5u7/yahaxURkU7JDPnu9Cb47zOzi4C7gJZ9je7e9ZLLIiLSg0SMzydCb4J/37V1vhnT5sB7El+OiEjmS+b4fCK8a/C7+8RUFCIikmkGesDH09MPuE5w90fN7BPdve/uf01eWSIiA99AD/h4eurxHwc8Cny0m/ccUPCLSCQMlLH5ROnpB1yXh8+fS105IiLplepTK9PhXcf4zez73bW7+w8TX46ISGpk6vh8IvTmrJ6GmNeFwBnAK8kpR0QkNaIQ8PH05qyea2Knzexq4O9Jq0hEJIGybXw+EfpzB67BwNhEFyIicjCiPHTTV70Z43+R4CweCK7DXw5ofF9EBhQFfO/1psd/RszrNmCru7clqR4RkXel4ZuD05sx/g0AZlZBcHB3tJnh7m8luzgRiTYN3yRHb4Z6PgZcA4wGaoAJBGf1HJHc0kQC8f74491qTj3B7KGAT47eDPX8J3A08Ii7H2lmc4HPJLcsyWZ9CfKeenc9tXdt6+s6NZSQWvrvnVq9Cf5Wd99uZjlmluPuj5nZ9ckuTDJLf8K8t0GeCP3ZeXRtVzgdPA3dDAy9Cf5dZlYCPAncamY1QH1yy5KBKlG98kwUb1uue3gtD9+2GICXR1Vrh9CDbPo8ZLLeBP8qoBG4BDgXGAqUJLMoST/1zHrvkpOnckl5I6f/fLGGl2Jk2/Zkk94E/1x37wA6gFsAzOyFpFYlKaOAT55kDi+lI1QTdWxG0q+n6/F/GbgImNQl6EuBp5JdmCReFK46mG36uvNIxBlQUR7Oy1S7G1t5rbaO12qCUfjP/f45XquNPyLfU4//NuAB4CfA/Jj2Ot1vd2BTDyy6EnEGlD4nA1NLWzubdjaxYXsD67c1AnD2giW8VtNAc2s7k8qLmVQRjMKfM2M8kytKmPQf3S+rp+vx7wZ2A+ckegMkMRTwItlnzZY61m9vCAJ+exDws696lNq6FkYPK2TCyGIqRw4GYN7cyUyuKOGQIYWYGQB/XV7NKUcc0uM6+nORNkkxBbxI9qhrbgVg4Qub9wv4DduDK+DPu205E0YMZsLIYt57SCkAt10wk9HDisjPzelczi1LNnDslPJ+1aDgzwAKeJHM4O7saWpj485Gqnc1AfCD+16iemcTm3Y2sWlnI63twTUv711VTeXIYqaPG87Hp49hQlkxs696lEcuPW6/ZX7/npeYMLI4oXUq+EVE+mBHw1427Wxk086mMNCD4ZjTrn+STTubMGDM8CLGDi8CYPTQImZUjgjbBjN8cD4TL7ufG/6tKm3boOAfYHTus0h67W5s5c1w2OW3T7zOpp2NnT12gLlXP87YMNjHDBvc2Ru/5qxpjB0+mKFF+Z3Lqpy/kAs+9J7Ub8S7UPAPMANhWMfdaWnroKEluPr26urdNLS0UR/zALj2oTXUtbR1vlfXHLT/y2+epq29g9Z2p62jg7bwq+3RVy6irSNsb+8A4PDvP0iuGWaQm2PkhAeoZl75SNhuYXtQ25m/+gcF+bkU5udSlJ9DUX4uAFc98CqlhXkMKcyjtDD4w1u1cRcjSwZRVlJAYTifCEDj3jbWb2tkfRjwX79zFeu3N/Dmtgb2tnVQWRYcPN3RsJcpFaXMPayCMcOLOO36xay6/JQDlvfDv73MEaOHpnQbDoaCP02S3bPv6AjCdv22BnY3tbKrqZXd4QPgxwtfZlfj/m1zfvpoZ4gbRklh8PH45p9foLQgj+KCXEoK8ykpCEI0J8cYM6yI0sI8iguCx+J125j/4feSl2Pk5+aQl2vk5eRw0rVP8NeLjiEv18jPCdrff8VDPPedk+hwp6PD6XBo73CO+vEj3D1vNh1O2O60dzgnXPMEV3zsCJpbO2hua6d5bztNre3cvfJtSgvz2NPcSvWuJurDHdB3717N9voWtjXsJT/HGFlSAMAXbllKWckgRpYMYmRxASNLBgGwcUcjo4YUMijvnQNokrncndr6FtZuqWft1joAzlnwDG9ua2Bn414mjBzMxLKgtz5j4nA+fdQ4KssGU15SgJlROX8h3z79fenchKRJS/Cb2XqgDmgH2tw9fYNdaXKwPfu9bR1s2N7AuvAHG1fc+xJbdjezZU8zW/c0s62+BYDP3vQcQ4vyGTY4nyFF+Z1fQ8tKCphcUcLQonyGFg3inN89w+0XHB0GeC4FeUG4V85fyAMXH3vA+m9/biNfO6n7HdVRlSO6bR89rOiAtpKC7j+Chw49cF6AI8cPP6Dt0jtXMW/u5P3a7l31Nvf9+xwgCID6lja21+/l+Ksf51NVY9lev5ft9S28taORZW/tBODsBc9QU9fMiOJBjB5W1Fnvn5dtYmJZMe8pK2Z48aBu65L02tmwlzVb61gXBvxZNyxh7dY6DJg6qpSpo4KzYy6aO4mJZcUcOrSI3PBrZOX8hXz6qPHpKj0t0tnjn+vu29K4/oyxo2EvK8Jw+vL/LGNdTT1v7WhkzLAiJoc/2Bg3YjBHVY7gkKEFjBpSSEVpIVO/+wBPfmvuAcu77dm3+OJxkw5oHzdicHI3JE3MjNLC/M4hoFO7Ocd54QsLeWr+CbS1d1BT18Lbu5qo3tXEwhc28491tfxhyXreqG0gNzcIi0vvXMmk8hImlhUzsayYyraOlG5TVO1pbmXd1jrWbq1nzZYg5Kt+9Agtre1MPaSUqaOCv4eLT5zC1FGllJUM6jy//Q/P9P/0x2yjoZ4k6+uQjrt39uK/8adVLN+wk5q6FqaPGwbAh99/KF+tCAJn37h15fyFnD9nYvI2IkLycnM6e/tVwMV3rOT6s48Egv832xv2UvWjR5g5cQRvbGvg7hXVvLmtgbe219NeNp5//n9PMXpoEYcOLeTQYUWMHloIQM2eZspKCsjZd7BCelTf0sa6rXWsq6nv7MUf85NF7GpqZUpFCVNGlXJY2Iu/799n7/cDptuf28jsyWVpqz0TpCv4HXjIzBy4wd0XpKmOpOvNkE5DSxtPv76dR1+t4fE1NZ1fQY8cP4zz50xk6qhScnOCMcePTRudyvIlhplRFh4n6Do00L58Baf++hm+Pe9YNu9uZvOuJjbuaOTZN7YDcPovFrOnqY2KIQXBjmFYsEO4cfEblJcWUFFaGDwPKUjtRqVZXXMrr9XUs25rPetq3gn4nY2tTKooZkpFKVPCXvwfvziLMcOK9tt5/vj+V+IOC0p86Qr+Oe5eHd7H92Eze9Xdn4ydwcwuBC4EGD8++8bftte3sPDFzQDMvHIR08YNZe5hFZw/ZyaTyouZeNn9nDtzQpqrlN7KzTEGdbR1e3yjcv5Cln73ZJpb29m6p5m3dzWzeXcT96x8m+pdTazcuIuauhZqwwfAsT97lIrSQipKCygvDXYGdz6/kfJwumJIASOLM2cnsWZLHW/taAwe4Zk0sQE/taKUyWHA33HhLMYMf2cMHuBnD67J2qHIdEhL8Lt7dfhcY2Z3ATMIbvQSO88CYAFAVVWVp7zIPurNkE7T3nYeenkLd6+oZumGnZz43goAllx2Quf4s2SvwvxcJows7jzv+9I7V3H5Rw+8dXXl/IX84fMzO3cGNXXNADy3fgc1dS3UhAfvdzUGZ2OdcPXjDCl65+D90KLgz/qGJ14Pp4PHkPDAfs2eZgYX5DE4P7ffQ0+1dS3UNbdS19wWPoJarnloDTV7gppr6lrYuifYkc27bTnjRwwOHuH233HhLMYO378H/7MH1zB+pAI+2VIe/GZWDOS4e134+hTgh6muI9HiDem4O8++GVzMdMaVj3Dk+OF84sgx/OpfP0hxQV54KqJCX/ZXWVZMZdk7P9P/wX0vc/Wnpu03T2t7B1O+8wALPlvF7qZW9jS3sifm9Nxt9S28XlsfvNfU1tn+kV/+g8aWNhpb2ykMz96aeeUj5OXkkJ8b/G5i3zVh5l79OC2t7bS0ddAcPgOcev2T4e8m8iktzKM0PPU3N8eYNm4YFeG3klFDCpl55aIDLkPwn397WQGfRuno8Y8C7goPxOQBt7n7g2moI6ka97Zx14pq/v/TG2jrCP5YFn39OCpKC9NcmWSLfeG878yuWN+/5yW+85HDD2ivnL+Q579zEhD8RqKptZ0jLv87d8+bTVu709bhnT++O/0Xi7nxvCoK8nIozM+lIC+Hgrxcpn73AZZ/7+Rulx3vFF8ZWFIe/O7+BjDtXWfMYD/628v8efkmqiaM4HtnHM7sySOZeNn9Cn0ZUHJyjOLwdxTxDpBOKtddVrORfqLYR9c9vJbK+Qs7H9c9vBYIDl7Nu3U5EHzdve8rc7jxvCrmTCnrPM0sWbUA+9XS13YRiRadx99HXcfyX92yh4tuXcZzb+7kgmMnsvDFzVx2ED/zjj1IXDl/4QG3yNvXFns9/ng19qa9u/V1vf3ewbSLyMCj4D8IX/rDMpZu2MmFH5rI1Z+axuBBefzkgVd79W/jBWVPoZ0Midh59NTelx2Zdh4iqaHg76Pauhb+6+9BuFdVDue6T0+naFD8Kz/2NeCzTV92ZH3ZeWhHIdJ/Cv4edD03/0NTylj99h4+ceQYAL5w7LtfZzsqAZ9MiRi6EpF3KPh7EDueP6WiBAfu/OLRTK4o5cZ/vLnfvAqcgaGv3xpEokjB34MdDXv53j2rAfjGqYdxyuGj4p6ho579wKZhJJF3KPjjeHJtLd/68wt8dNqhwDuX8lUgZBcNI0kUKfi7aG5t56cPvsqDq7dw7VnTOGZyGb9b/M6wjnr20dXTt4aHb1sMaIcgmUHBz4EHcadUlPDAxccybLDutiTv7pKTp3JJeWMwMX16Z7u+IUgy9GV4Mh4FP8EfblnJIL53z0tc/alp/MsHx2Bm+sOVg6LjCtkvmT+AhIP/4ealceqOfPC3tXfww7+9zJLXgxtmfPL/jO18T8M6kgw6rpBaAyWE+9ueDJEO/j3Nrcy7dTlmxl8uOoYPXPFQuksS2U/UvzUMtHDOFuY+4O9xQlVVlS9dujShy3xreyPn3/I8syaNZGhhPr987LXO9zL9j0XSYOXK4DlmjH8g6e5GQcABbV1DNbZ9n3i3Ee3ansh1Sv+Y2TJ3rzqgPUrB3/XDdfzUcm7+/IyDXq7IQA/+RIgXzgrtgUvBH7r/xc1cdOtybv7cURx/WEVClikSheCXzBMv+CM1xn/fqrf5wX0vAyj0RSSyInMjlntWVvMff3mBbfXBzZ91IxIRiapI9Pj/smwTP33wVe6eN5upo0rTXY6ISFplffDfuXQj1z60ltsumMnkCoW+iEhWB/8dz73Fzxet49YLZuqm0SIioawM/q6nl9278m2dXiYiEsrKg7vVO5vSXYKIyICVdT3+12vreXxtLX84fwbHTilPdzkiIgNOVvX4t9W38LnfP8+3Tj1MoS8iEkfWBH/T3na+cMtSzpw+mrOOGpfuckREBqysCP72Dudrf1xB5cjBXKqDuCIiPcqKMf6f3P8Ku5ta+cU5R8a9GbqIiAQyPvhPuvYJXqupB+Cw7z7IzIkj+OMXZ6W5KhGRgSujg3/jjkZ2NuzlnnmzmTZuWLrLERHJCBk7xt/a3sFX71jBl4+fpNAXEemDjA3+ax5ay7CifD4/e2K6SxERySgZOdTzxNpa7l5RzcKvziEnRwdzRUT6IuOCv2ZPM9/40yp+cfaRjCwpSHc5IiIZJ6OGejo6nEvuXMk5M8Yza9LIdJcjIpKRMir4f/PE67S2OV89YXK6SxERyVgZM9Tz6RuW8OybOwCY/J0HdL6+iEg/pSX4zew04OdALnCju1/1bv9GIS8ikhgpH+oxs1zg18CHgcOBc8zs8FTXISISVekY458BvObub7j7XuAO4Mw01CEiEknpCP4xwMaY6U1hm4iIpMCAPavHzC40s6VmtrS2tjbd5YiIZI10BH81EHunlLFh237cfYG7V7l7VXm57qYlIpIo6Qj+54EpZjbRzAYBZwP3pqEOEZFISvnpnO7eZmZfAf5OcDrnTe7+UqrrEBGJqrScx+/u9wP3p2PdIiJRN2AP7oqISHIo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGIU/CIiEaPgFxGJGAW/iEjEKPhFRCJGwS8iEjEKfhGRiFHwi4hEjIJfRCRiFPwiIhGj4BcRiRgFv4hIxCj4RUQiRsEvIhIxCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYkYBb+ISMQo+EVEIkbBLyISMQp+EZGIUfCLiESMgl9EJGLM3dNdw7sys1pgQzhZBmxLYzmpEoXtjMI2grYzm2TaNk5w9/KujRkR/LHMbKm7V6W7jmSLwnZGYRtB25lNsmUbNdQjIhIxCn4RkYjJxOBfkO4CUiQK2xmFbQRtZzbJim3MuDF+ERE5OJnY4xcRkYOg4BcRiZiMCX4zO83M1pjZa2Y2P931JIqZ3WRmNWa2OqZthJk9bGbrwufh6awxEcxsnJk9ZmYvm9lLZnZx2J4122pmhWb2nJmtCrfxB2H7RDN7Nvzs/tHMBqW71kQws1wzW2Fmfwuns247zWy9mb1oZivNbGnYlvGf2YwIfjPLBX4NfBg4HDjHzA5Pb1UJczNwWpe2+cAid58CLAqnM10b8HV3Pxw4GpgX/j/Mpm1tAU5w92nAdOA0Mzsa+ClwnbtPBnYC56evxIS6GHglZjpbt3Ouu0+POX8/4z+zGRH8wAzgNXd/w933AncAZ6a5poRw9yeBHV2azwRuCV/fAnw8lTUlg7tvdvfl4es6gsAYQxZtqwfqw8n88OHACcCfw/aM3sZ9zGws8BHgxnDayMLtjCPjP7OZEvxjgI0x05vCtmw1yt03h6+3AKPSWUyimVklcCTwLFm2reHwx0qgBngYeB3Y5e5t4SzZ8tm9HvgW0BFOjyQ7t9OBh8xsmZldGLZl/Gc2L90FSM/c3c0sa865NbMS4C/A19x9T9BRDGTDtrp7OzDdzIYBdwHvTW9FiWdmZwA17r7MzI5PcznJNsfdq82sAnjYzF6NfTNTP7OZ0uOvBsbFTI8N27LVVjM7FCB8rklzPQlhZvkEoX+ru/81bM7KbXX3XcBjwCxgmJnt62Rlw2d3NvAxM1tPMOx6AvBzsm87cffq8LmGYEc+gyz4zGZK8D8PTAnPGhgEnA3cm+aakule4Lzw9XnAPWmsJSHCMeD/Bl5x92tj3sqabTWz8rCnj5kVAScTHMt4DPhkOFtGbyOAu1/m7mPdvZLgb/FRdz+XLNtOMys2s9J9r4FTgNVkwWc2Y365a2anE4wr5gI3ufuP01tRYpjZ7cDxBJd73QpcDtwN3AmMJ7gc9Vnu3vUAcEYxsznAYuBF3hkX/jbBOH9WbKuZfYDgYF8uQafqTnf/oZm9h6BnPAJYAXzG3VvSV2nihEM933D3M7JtO8PtuSuczANuc/cfm9lIMvwzmzHBLyIiiZEpQz0iIpIgCn4RkYhR8IuIRIyCX0QkYhT8IiIRo+AX6QMzu87MvhYz/XczuzFm+hozuzQtxYn0koJfpG+eAo4BMLMcgt9fHBHz/jHA02moS6TXFPwiffM0wWUYIAj81UCdmQ03swLgfcDydBUn0hu6SJtIH7j722bWZmbjCXr3SwiuQjkL2A28GF46XGTAUvCL9N3TBKF/DHAtQfAfQxD8T6WxLpFe0VCPSN/tG+d/P8FQzzMEPX6N70tGUPCL9N3TwBnADndvDy/QNYwg/BX8MuAp+EX67kWCs3me6dK22923packkd7T1TlFRCJGPX4RkYhR8IuIRIyCX0QkYhT8IiIRo+AXEYkYBb+ISMQo+EVEIuZ/AV8y+vDuzj7bAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pyerrors.jackknifing as jn\n",
    "jack1 = jn.generate_jack(c_obs1, max_binsize=120)\n",
    "jack2 = jn.generate_jack(c_obs2, max_binsize=120)\n",
    "jack3 = jn.derived_jack(lambda x: np.sin(x[0] / x[1] - 1), [jack1, jack2])\n",
    "\n",
    "print('Binning analysis:')\n",
    "jack3.print(binsize=25)\n",
    "jack3.print(binsize=50)\n",
    "jack3.print(binsize=100)\n",
    "\n",
    "jack3.plot_tauint()\n",
    "\n",
    "print('Result from the automatic windowing procedure for comparison:')\n",
    "c_obs3.gamma_method(S=1.5)\n",
    "c_obs3.print()\n",
    "c_obs3.gamma_method(S=2)\n",
    "c_obs3.print()\n",
    "c_obs3.gamma_method(S=3)\n",
    "c_obs3.print()\n",
    "\n",
    "c_obs3.gamma_method(S=2)\n",
    "c_obs3.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this specific example the binned Jackknife procedure seems to underestimate the final error, the deduced intergrated autocorrelation time depends strongly on the chosen binsize. The automatic windowing procedure displayed for comparison gives more robust results for this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}