pyerrors/examples/03_pcac_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pyerrors as pe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('./base_style.mplstyle')\n",
    "plt.rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Primary observables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can load data from preprocessed pickle files which contain a list of `pyerror` `Obs`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_obs_names = [r'f_A', r'f_P']\n",
    "\n",
    "p_obs = {}\n",
    "for i, item in enumerate(p_obs_names):\n",
    "    p_obs[item] = pe.load_object('./data/B1k2_' + item + '.p')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the `pyerrors` function `plot_corrs` to have a quick look at the data we just read in "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApcAAAGLCAYAAACSkLnsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAiuUlEQVR4nO3dYYhd12En8P+Rw5pu8MzECkvTdYQ7ap18apHsgL8YdjcyLgRMIZKzRaRsG2pBV/5k1rI/LFXZD468uGBWW5BpS0kcwkYOLIJAjb3dbo0/1ZZIPhm3o5jYrWGJ7JkJ2eLdWmc/vDf20+jNaN5797737pvfDx537j3v3Xvsp9H965x7zim11gAAQBMOzLoCAAAsDuESAIDGCJcAADRGuAQAoDHCJQAAjREuAQBojHAJAEBjPjXNi5VSSpJfSvKzaV4XAOi8O5L8QzVB99ybarhML1i+O+VrAgCL4a4kfz/rSrC7aYfLnyXJO++8k6WlpSlfGgDoos3NzXz+859P9Hx2wrTDZZJkaWlJuAQAWEAG9AAA0BjhEgCAxgiXAAA0RrgEAKAxwiUAAI0RLgEAaIxwCQBAY4RLAAAaM5NJ1Dvn2lpy5dvJ+k+SlUPJka8nBw/PulYAAHNHuLyVKy8klx5LUpLU3va155KHzydHTs64cgAA80W43M21tV6wrNdvLrt0Ojl0vxZMAGBP7n7yB59L8rld3vLe29/8ynvTqk9bPHO5m1efnawcAOATp5K8scvr1Oyq1hwtl7t566XhrZZJ7/hbL023PgBAl11Icqn/8xeTfCfJySRv9o91vtUyES53d89DyQ+/OzxglgO9cgCAPeh3eb+XJHc/+YOtw2++/c2vXJ5FfUopx5KcqrWeaPK8usV388Djk5UDAMyZUsrRUsq5JCeSrDZ9fuFyNwcP90aFlwNJue3G7cPnDeYBAEZ295M/+NUkp/u7p/v7U1NrvVxrPZPk5TbOX2qtbZx3+MVKWUqysbGxkaWlpaldd2LmuQSAmdnc3Mzy8nKSLNdaN2ddn0nc/eQPfifJn6Q3v+FtST5Kb77Db7z9za/8+TTrUko5nuSpWuu9jZ5XuAQA5tmihMt+C+WbGd5zfD3JF97+5lf+blr1aStc6hYHAJiO300vRA5Tk3xjinVpzeKMFtd1DQDMt7uze8Pe3dOpRrsWI1x2eYlGoRgA9ou302u53Clgvj21mrSo++Gyy0s0djkUAwCj+rMkT+xQVpL86RTr0pruP3PZ1SUaB0Nx/ejG7aXTvXIAYGG8/c2v/G16z1VeT2+UePrb6+mNFp/aYJ42jRQuSykXSyl1yOuNtip4S11dorGroRgAGFt/uqEvJPlW/9C30hsl/uczqM6dbZx01JbL95M8mOTegdczSX6v4Xrt3T0P9SY1H2ael2jsaigGACbSb6E83989P+0Wy4EVes4kOVpKuVBKebSp84/6zOXLtdZXBiq3kuRarXUma2Im6S3B+MPv7l4+j6xbDgDMQD+3XU4vXDZupJbLWuuL2w49VWt9psH6jK6rSzRatxwAWEBjjxYvpRxN8je3eM/tSW4fOHTHuNfb1ZGTvVHhXZrSZysUXzqdG0aLp853KAYAxnL3kz/4XJLP9Xe/uLW9+8kfbL3lvbe/+ZX3pl6xhk0yFdFTtdYTt3pPkj+Y4Bp7d/BwcuzsVC7VmC6GYgBgXKdycy76zsDPf5jk7NRq05KxwmUpZXWPb306yR8N7N+R5N1xrrmwuhiKAYBxXEhyaZfyzrdaJuO3XB5PcvVWb6q1fpjkw639UsqYlwMA6LZ+l/dCBMjdjDuJ+teSmOUbAIAbjBsuV9Ob8xIAAD42brhcSbLeXDUAAFgE44bL9Wi5BABgm7EG9NRaP9N0RQAAFtrZ5cF5Lod5L2c3Oj/gZ5J5LgEA2Lth81wO2r/zXAIAMLLBeS6/mN4E6ieTvNk/1vlWy0S4BACYjl6Xdy9Anl3eOvpmzm5cnlWV2iBcAgDsQ6WUJ/o/finJ1VrrmSbOK1wCAOwzpZRzg2GylHKxlHKx1npi0nOPOxURAADjOLv8q0lO9/dO9/enppSykuRYf7vl6STHSymrk55fuAQAmJazy7+T3gCe3+4f+e0kb+bs8r+bck1W+68tVweOT0S4BACYhl4L5Z+kl79u6x+9rb//pzm7/CvTqEatdb3W+pla6+BAoq1QeXXYZ0bhmctFdW0tufLtZP0nycqh5MjXk4OHZ10rANjPfjfJ9Qxv3KtJvpHkqanW6BOnkrxSaxUuGeLKC8mlx5KU9P6sluS155KHzydHTs64cgCwb92d3XuN755ONW5USjma5FiSe5s4n27xRXNtrRcs6/WkfnTj9tLpXjkAMAtvp9dyuVv5LJxLcm+tdb2JkwmXi+bVZycrBwDa8mfZOXuVJH86xbr0LlrKhSSnmgqWiXC5eN56qddKOUy93isHAKbv7Mbfpvdc5fUkH/WPftTf/0bObvzdNKtTSnk0ybmt5yxLKav9LvKJCJeL5p6HkrLD11oO9MoBgNk4u/HnSb6Q5Fv9I99K8oX+8akppRxPspJktZRyrL9/Jg2MFi+11knPsfeLlbKUZGNjYyNLS0tTu+6+cm0tOX/f8NbLciA5/bpR4wB0yubmZpaXl5Nkuda6Oev6NOLs8tEkbyS5d9pri/cnT/9gWFmttUx6fi2Xi+bg4d6o8HIgKbfduH34vGAJAPtcf57LMuzVxPlNRbSIjpxMDt1vnksAmCdnlz+X5HP9vS9+vD27vPWO93J2472p16thwuWiOng4OXZ21rUAAD5xKskfbDv2nYGf/zDJ2anVpiXCJQDAdFxIcmmX8s63WibCJQDAdPS6vBciQO7GgB4AABojXAIA0BjhEgCAxgiXAAA0ZqwBPaWUJ5Ks93ffr7W+2FiNAADorJHDZSnl5SSnaq1X+4ubv5GkkRndAQDotpG6xUspjya5XGu9miS11stJ7m2jYgAAdM+oz1yeS/Ly4IF+wAQAgL13i5dSVpKsJFnpt2AmyeFa65kW6gUAQAeN8szlan97Z631+SQppRwrpVystZ4Y9oFSyu1Jbh84dMd41QQAoAtG6Ra/s799fetArfWVJMdLKavDP5KnkmwMvN4dp5IAAHTDKOHy6rbtlvUkR3f4zNNJlgded41SOQAAumXP3eL9qYeSXvf44CCelV0+82GSD7f2+58HAGBBjTpa/HI+6R7ffhwAgH1u1HB5JsnHg3f6o8Zf3Jr3EgCA/W2kFXpqra+UUlZLKecGjg0dKQ4AwP4z8vKPW9MQAQDAdqN2iwMAwI6ESwAAGjNytzjk2lpy5dvJ+k+SlUPJka8nBw/PulYAwBwQLhnNlReSS48lKUlqb/vac8nD55MjJ2dcOQBg1nSLs3fX1nrBsl5P6kc3bi+d7pUDAPuacMnevfrsZOUAwMITLtm7t17qtVIOU6/3ygGAfU24ZO/ueSgpO/yRKQd65QDAviZcsncPPD5ZOQCw8IRL9u7g4d6o8HIgKbfduH34vOmIAABTETGiIyeTQ/eb5xIAGKrUWqd3sVKWkmxsbGxkaWlpatcFALprc3Mzy8vLSbJca92cdX3YnW5xAAAaI1wCANAY4RIAgMYIlwAANEa4BACgMcIlAACNES4BAGiMcAkAQGOESwAAGiNcAgDQGOESAIDGCJcAADRm7HBZSllpsB4AACyAkcJlKeVYKaWWUmqSD0opa6WU1ZbqBgBAx3xqxPevJLm3//N6rfVqs9UBAKDLRg2XSXK11rredEUAAOg+A3oAAGjMOC2Xj5RS3u///KVa65md3lhKuT3J7QOH7hjjegAAdMSoLZdXk7xea32x1vpikrVSyoVd3v9Uko2B17vjVRMAgC4otdbxP9wbKb6W5DPDnsPcoeXy3Y2NjSwtLY19XQBg/9jc3Mzy8nKSLNdaN2ddH3Y3Trf4x2qtV0spSbKa5PKQ8g+TfLi1338vAAALas/d4qWUlVLKB4PzWppIHQCAQaO2XL6+bW7L1SSptd7UagljubaWXPl2sv6TZOVQcuTrycHDs64VALBHew6Xtdb1UsrL2w4/lWTH0eI
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pe.plot_corrs([p_obs['f_A'], p_obs['f_P']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Secondary observables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way of generating secondary observables is to write the desired math operations as for standard floats. `pyerrors` currently supports the basic arithmetic operations as well as numpy's basic trigonometric functions.\n",
    "\n",
    "We start by looking at the unimproved pcac mass $am=\\tilde{\\partial}_0 f_\\mathrm{A}/2 f_\\mathrm{P}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "uimpr_mass = []\n",
    "for i in range(1, len(p_obs['f_A']) - 1):\n",
    "    uimpr_mass.append((p_obs['f_A'][i + 1] - p_obs['f_A'][i - 1]) / 2 / (2 * p_obs['f_P'][i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more complicated secondary obsevables or secondary observables we use over and over again it is often useful to define a dedicated function for it. Here is an example for the improved pcac mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pcac_mass(data, ca=0, **kwargs):\n",
    "    return ((data[1] - data[0]) / 2. + ca * (data[2] - 2 * data[3] + data[4])) / 2. / data[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can construct the derived observable `pcac_mass` from the primary ones. Note the additional argument `ca` with which we can provide a value for the $\\mathrm{O}(a)$ improvement coefficient of the axial vector current."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "impr_mass = []\n",
    "for i in range(1, len(p_obs['f_A']) - 1):\n",
    "    impr_mass.append(pcac_mass([p_obs['f_A'][i - 1], p_obs['f_A'][i + 1], p_obs['f_P'][i - 1],\n",
    "                             p_obs['f_P'][i], p_obs['f_P'][i + 1]], ca=-0.03888694628624465))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the error of an observable we use the `gamma_method`. Let us have a look at the docstring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "?pe.Obs.gamma_method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can apply the `gamma_method` to the pcac mass on every time slice for both the unimproved and the improved mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "masses = [uimpr_mass, impr_mass]\n",
    "for i, item in enumerate(masses):\n",
    "    [o.gamma_method() for o in item]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now have a look at the result by plotting the two lists of `Obs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA10AAAGLCAYAAAA8g/hEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABDwUlEQVR4nO3df3Ac533n+c9DOSFVFIEhSMO2RDI0SEkuWyuRICQn2rgUHoHbqJRVshZ/2MVSTlqfwKsLtZLLKYL8I2Vmry4UWJHX2vCqjlCyUsVWOfwh7kYX1VaOUiSWf+QsQxCtcqociQQZkZZVFH8MQLEkulb83h/9NNgYTg8GwHT3TPf7VTXEdD/zdD+NAcD+zvM838eZmQAAAAAAyZiTdQMAAAAAIM8IugAAAAAgQQRdAAAAAJAggi4AAAAASBBBFwAAAAAkiKALAAAAABJE0AUAAAAACfpE1g1oRs45J+lGSRezbgsAAGgpCyS9ayyECiCCoKu6GyWdzroRAACgJS2R9IusGwGgeRB0VXdRkk6dOqW2tras2wIAAFrA+Pi4li5dKjFSBkAFgq4a2traCLoAAAAAzAqJNAAAAAAgQQRdAAAAAJAggi4AAAAASBBBFwAAAAAkiKALAAAAABJE0AUAAAAACSLoAgAAAIAEEXQBAAAAQIIIugAAAAAgQQRdAAAAAJCgT2TdAAAAADSH5dtf/Iykz9R4yS9PPnHfL9NqD5AXBF0AAAAIbZH0zRrlfyppZzpNAfKDoAsAAAChvZJe8M8/J+k5SZsl/dzvo5crQc65QUldknaZ2UjW7UHjEHQBAABAkuSHDv5SkpZvfzHc/fOTT9yXagDgnOtW0OvWL6ksaUjSXjMbTbMdGdgl6YSC4Bc5QtAFAACApuJ7ebY457okjZrZQNZtSoOZlZ1zeQ8sC4nshQAAAJhk+fYXb5a01W9u9dtIx/msG4DGI+gCAADAhOXbX3xYwRyuP/S7/lDSz5dvf/GhzBoFtDiGFwIAAEDSRA/XX2ryB/PX+a9/tXz7iz84+cR9x9JvWcA51ytp0G8+oiDpRIekNWa2xTnXr6CnaJMiySicc+sl7fD1dvk6JUmLwqGL/jWDkkYUzKnqk6RI+TZJ4dC/LjPb7c834NuxwcwO+vloL/vXbjCzUZ8g4yf+dWUzG4pc0zb/NDx2Rx3fhymvp+L4Zf99kZkd9Pu7JHX7l90p6bCZvVRP3Zg2zei9qactkbpS5HsYtz+ujVki6AIAAEDo30u6ouqjoUzS13T1Zj91ZvaSc25AQVDUEQkgjjvnBiMBkiQ9LWmNr3fQ7zsgaSRMyOGcG3TO7TWzLZHX7FBwI79PQYAg59xeSQfCQMA51+WcO2xmfX4O1t6wLWY24pzbZWa7/WsPSxqM1D3gnBv11zIo6VzktSXf7qm+D1NeT6Tdx8NAxJ+7y59vr4LgZrekg/57uMbMynXUbdh748W2xQeYlcHi+rj9U33vskLQBQAAgNBy1Z5+sjydZtR0XkFPU7RXpjL5xIiCnp+osoKkHNHX7pJ0wQcFo/41XZFemBF/M78xDGQkyfdedTjnen2w0eGc647UK0sTWRh7Ktp6WEGSkGFJ28zMRY47nUQaNa/Hb/dHj6+gVy7sGdqiyfPHRiX1Kgh6uqaoG2em701sW/z2BufcfjMr++/9S/K9i1X2NyXmdAEAACB0UkFPV63yZlCusn28Yl+1YXqT6vlenbKuDm2Trg0SeqvsC1/X558PKQgcwmF2+/3+HknnnXPrw4eCgOMn/riV1zFdk+pXXE93lfLRsCfLB2sdzrl+364OXf2e1aw7nTapjvemVlsiQxovOOded85tM7ORuP11tC8T9HQBAAAg9F8kbYspc5L+KsW21NJsGf72SnpdQeBV2dNTdS5UODwuK36+1p1mtsFvb2rQoaf93kzVFj+Ms0tBoDrgnFtkZgNx+xtwDQ1HTxcAAAAkSSefuO9tBfO2rkj62O/+2G9/LcskGg1Sim74OVQlBUPe4gwrGMpWqVtBj1XYUzNakdghtq4PFKoNs5uuSfUrrqfq8f18tJKCeWYbKo/lh0TG1p1le68xVVt871fJ97QNKZgL1hu3v9HtaxSCLgDFcPE96d2j8Y+L72XYOABoHiefuO9ZSbdK+mu/668l3er3N4vKoYOlKbZDYcAR2iFpqGJe1CR+yNpL0Z4pH5iUK3qw9ioIHg5W1D3og7GwbklStz/nUEVZmMUvrv11X0/k+BM9l/61YVZBVdQNv6ddU9SdynTfm5pt8a/vn1xFozX2NyVnZlm3oek459okjY2NjamtrS3r5gBohFd2SUeeiC+/Z7u0NrOEXAByYHx8XO3t7ZLUbmbjWbdntpZvf7FbwZC5NSefuC/VuTI+qNmi4Ka6rGDO1F4FN9qDCno0hhQkd9ihYEjkiIJEEqOVr/EJKnr9MQYUmccVyRwYpjzv9vUGo8GYT04Rzk1aUSU1e0nSjmrD28Ishf6856OBWSQtuxT0ku3w1zlQR4r22OuZ6tz+vCsUJPYo++9bmKVxqFbdmPZ0a+bvTWxbIqeIpoYfkrSx2v46552ljqCrCoIuIIcuvne1N+vsW9KhR6QvPy0tviXYt+DTwQMAZoigq7mFQZWZrZnyxS0gb9eTdyTSAJC+aABUTRIBULVjLr5FunFVY88DAGhmpawb0GClrBuA+hB0AUjf8DMM9QOAJrR8+4ufkfQZv/m58Ovy7S+GL/nlySfu+2XqDZsl3ys0oGAO1GC1IYCtJG/XUwQEXQDS1/OwdOu9wfO4oX4AgCxskfTNin3PRZ7/qaSdqbWmQXwK96ZdOHe68nY9RUDQBSB9DPUDgGa1V9ILNcpbrpcLaAapBF3RVJiSSpVZVWZSp55jVmaDqZVxBQAAoOj80EECK6DBEg+6fHBUiqTjXD/V2NOp6tRzTOfcYUlbzGzUp7B8XcFK6gCARssiOQoAAC0ijZ6uAUl94YaZHXTOPe33z7ROzXIflI2EayuY2YhzjnSaAJAUkqMAABAr0aArXLm6yirfJedct1+le1p1FCyWNtUxByVtiBZWOxcAJKpIvT8kRwEAIFbSPV1dMfvLvqxaIDRVnThlBWkzRxWsWVCKzPu6ZtXwKOfcXElzI7sW1DgPANSnSL0/JEcBACBW0kFXR8z+8zXKpqpTnqI8DMw6zGxICtYycM4dMLMNMXV36Nr0qADy6Nxx6bWh4PlrQ9KXviEtWpHMuej9AdBqdrZH1+mq5pfaOUaiDWCa8pgyPgzahsMdZvaSc+6wc67asERJ2iXpW5HtBZJOJ9hGoDkUafibJL3xXemFR69u//R7weP+PdLqzY0/H70/ySrazy+QjmrrdEW15DpdQNaSDrrOx+zvqFE2VZ2pysOgqjK4KkvqrrJfZnZZ0uVw2zmSHKIgijT87dzxIOCyK1f3hc9f2Cot+83keryQjCL9/ALpia7T9TkFCyNvlvRzv49eLmAGkg66RqUgOYaZlSP7S6oS/NRZp2a5TxEvXTtnrDSjKwDyrEjD3974juJXjXBBee/OFBuEWSvSzy+QlmDoYBBY7WwP9/5cO8dISNaCnHODCu6Jd5FULluJBl1mVvaJLa6ZixX3xtdTp45jjqj63DB+2ICorIe/pTm/qvyOJIspNF+OlpLVzy/DGoHEOee6FAx13Kbgfm9IQeBQ9uV7JfVLeknSoJm9VOdxBxVkwY6b5583uySdUNCDiQylMadrUNJ6SeFCxv2KrNHlf6l6w6QX9dSpo3xAQcr4lyLlB2PmcwHIQtrzq0rLVLOnq7Ss8edEPjGsEUicv2cbCJcLqsxCbWZbnHPlWtmpY/xE0rlGtbPZRTozkLHEgy4zG3LObfOBT0nSoopfkF4FQdJQvXXqKH/JOdflP80I9xXlEw2g+WUxv2r1g9IPn4optKA8L+iJSRbDGlEEO9tvlrTVb23VzvZd2jn2dpZNqmLawZOZHUyiIU0uLh8CUpRK9kIz212jbEiRgKueOnWWX3NMAE0ii/lVi1YEvWgvbA3OYR9Lbk5Qdv+efCXRoCcmWVkMaySQRpp2tj8s6S91dUz2H0r
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pe.plot_corrs([impr_mass, uimpr_mass], xrange=[0.5, 18.5], label=['Improved pcac mass', 'Unimproved pcac mass'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tertiary observables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now construct a plateau as (tertiary) derived observable from the masses. At this point the distinction between primary and secondary observables becomes blurred. We can again and again resample objects into new observables which allows us to modulize the analysis. Note that `np.mean` and similar functions can be applied to the `Obs` as if they were real numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 4.79208242e-03 +/- 2.09091228e-04 +/- 1.90500140e-05 (4.363%)\n",
      "  t_int\t 1.09826949e+00 +/- 1.84087104e-01 S = 2.00\n"
     ]
    }
   ],
   "source": [
    "pcac_plateau = np.mean(impr_mass[6:15])\n",
    "pcac_plateau.gamma_method()\n",
    "pcac_plateau.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use a weighted average with given `plateau_range` (passed to the function as kwarg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def weighted_plateau(data, **kwargs):\n",
    "    if 'plateau_range' in kwargs:\n",
    "        plateau_range = kwargs.get('plateau_range')\n",
    "    else:\n",
    "        raise Exception('No range given.')\n",
    "    \n",
    "    num = 0\n",
    "    den = 0\n",
    "    for i in range(plateau_range[0], plateau_range[1]):\n",
    "        if data[i].dvalue == 0.0:\n",
    "            raise Exception('Run gamma_method for input first')\n",
    "        num += 1 / data[i].dvalue * data[i]\n",
    "        den += 1 / data[i].dvalue\n",
    "    return num / den"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 4.78698515e-03 +/- 2.04149923e-04 +/- 1.85998184e-05 (4.265%)\n",
      "  t_int\t 1.06605715e+00 +/- 1.79069383e-01 S = 2.00\n"
     ]
    }
   ],
   "source": [
    "w_pcac_plateau = weighted_plateau(impr_mass, plateau_range=[6, 15])\n",
    "w_pcac_plateau.gamma_method()\n",
    "w_pcac_plateau.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case the two variants of the plateau are almost identical\n",
    "\n",
    "We can now plot the data with the two plateaus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA10AAAGLCAYAAAA8g/hEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABUM0lEQVR4nO3de3gc13nn+d8BKd4BNEESliiShkDqMrZGFxCSM76MrAic2CtHyWPxYo/WM9Z4BDw7kUf2kyxBzvPMhsk+GwrcKGNNmNkl5Iy0sTUOL9ImijUeh1QkxbITyxBEaeVYkUiQIqCLKRJogKR4EYizf9RpsNDo6i4AXV3d1d/P84BA1elTdaobaNbb55z3GGutAAAAAADRqIm7AQAAAACQZARdAAAAABAhgi4AAAAAiBBBFwAAAABEiKALAAAAACJE0AUAAAAAESLoAgAAAIAIzY67AeXIGGMkLZd0Ku62AACAilIr6R3LQqgAfAi6clsuaSDuRgAAgIq0QtLbcTcCQPkg6MrtlCT19/errq4u7rYAAIAKMDIyopUrV0qMlAGQhaArj7q6OoIuAAAAADNCIg0AAAAAiBBBFwAAAABEiKALAAAAACJE0AUAAAAAESLoAgAAAIAIEXQBAAAAQIQIugAAAAAgQgRdAAAAABAhgi4AAAAAiBBBFwAAAABEiKALAAAAACJE0AUAAAAAESLoAgAAQEUzxqSMMc1xtwMIQtAFAACAsmKMaTbG7DLGWGPMkDGmy33tMsbsNca0+R7bIukZSfvja/HUGWPajDF7424HSmN23A0AAAAA/Ky1fZI6XO9Vn7W2M1NmjElJeskY02Wt7bbW9hpj7pM05QDGGNNure0uWsPDnbNF0iZJKUn0zlUJeroAAABQMay1aUm7JHX5dqenebh1M23PVFlre10QWVE9c5gZgi4AAABUmrSklOv1mjI3B6xLRe5pyp5bltmebjuRHAwvBAAAQKVZLW/YYTroAS74aXGbt0jab6094Lbb5AVczS74kqTtmeO5fT9zj0lnhiDmO6YbNpgJ5Fa7x2yVtFlSh6SSDmNEeSHoAgAASKimLU+3KdwQumNHH7zzT7Lq/pak+yT9sEDd/UcfvDMTzKhpy9PzJP1eUPlMucCnXdIdBR66S15QtEPSPmPMYWPMWmtt2lq7zxiTltTsny/mjr9fUpcvmNprjOlz2/mO2euCtV2ZY1lrO10whipH0AUAAJBc8+UlbChkMMe+RZJ+GqL+/Bz7UgXKp6LVGNPujrnE7bsqXy+X06GJ19Unr4drX1AFFyC1+nrEJG/uVYekAyGOmet5LNROVAGCLgAAgOQ6q3A3/acD9oWpezbHvnSB8qnomU6GQWttn5tPtVFeMNTgvvJplTRojFnv25eSN9RwuscECLoAAACSyg3rm9bQvuzhhlOod05SZ8EHRswYs1nSLdbaDW57U4HHjyfAsNbm7A2b6jGdVNg2I7kIugAAAJAoLltgl7XW+HanXFmLtbY3R7UWST3KkdHQBWSDIY6Z1uSeL3rCQNAFoEqces/7ClJ7ufcFACg3YYKWVK46xpiUb+5X5jjNknrlzcdqdo/LLMLca4zZ51802QVwLa5OoWMO+tvi6ma3bSrXhYQw1tq421B2jDF1koaHh4dVV1cXd3MAFMOz26XnHwwuv22LdPvW0rUHQOKMjIyovr5ekuqttSNxt6eSuSCoQ1669bS8dOvbcyXP8KVqb5O0I5ON0A0FXC0vEUZaXpC1S9JeX0CVSRd/2D9vzO0/6eoNZoYbhjxmJulHn3tMh7ygrNNae8C1d5Ok9W5/t6SXpjNvDZWDoCsHgi4ggfw9XSfekJ68T/riI9LSa7x99HQBmCGCLgBBGF4IoPTiGOqX65hLr5GW31Tc8wAAAGQh6AJQej2PMtQPAABUDYIuAKXXeq907ee9n4OG+gEAACQEQReA0mOoHwAAqCIlCbpcFpeMlLV2x0zrhDmmyzCTdpuDQQvdAQAAAEBUIg+6MmkzM0GRMWa9MaYrk85zOnXCHNMYs19Sh7W2z6XmfEmSmXQyAMDMsQ4aAACBStHT1SlpXWbDWrvPGPOI2z/dOnnLXVDWa63tc+W9xpi1RboeAEA2kqMAABAo0qDLrcLdnAl+fFLGmBZrbe9U68itHl7gmF2SNvgLc50LACJVTb0/JEcBACBQ1D1dzQH7064sVyBUqE6QtKRmY0yfvFXAU755X6sLDGecK2mub1dtnvMAQDjV1PtDchQAAAJFHXQ1BOwfzFNWqE66QHkmMGuw1nZLkjGmzRiz11q7IaDuVkm/G1AGIElOHpZe7PZ+frFb+sxvS0tWR3Muen8AAICSmTI+E7T1ZHZYaw8YY/YbY3INS5Sk7ZL+yLddK2kgwjYC5aGahr9J0svflZ76+qXtV77nfd21U7r5nuKfj96faFXb7y8AoGJFHXQNBuxvyFNWqE6h8kxQlR1cpSW15Ngva+15Secz28aQ5BBVopqGv5087AVcduzSvszPT90vrfqV6Hq8EI1q+v0FkJfLCdAQ8OE6ELuog64+yftDsNamfftTyhH8hKyTt9yliJcmzxlLTesKgCSrpuFvL39HwatGGK+8bVsJG4QZq6bfX6DKGGOa5WWlbpf3wbkbF66UvA/ad1lrD7jHtkh6xJVVxKdnbi1ZSbpF3v1rvqzeSIBIgy5rbdoltpg0Fysom2CYOiGO2avcc8PIYAj4xT38rZTzq9LHJNmAQuvKUVHi+v1lWCMQOddj1eGCr76stVhTkl5ya7R2u6WB7pO0d6rnMca0Z3IAlEqOtWX3Fsg9gASoKcE5uiStz2y4jIL+X7RmX5bBUHVClHfKlzLele+jyxkoIy9/V9rZ6s2pkrzvO1ullx+P5nypVcrb05VaFc15kTw9j0rdtwV/9TwadwuBRHMjnXbJux/MSE/zcOsKP6R4XMDY5r5nbJe03gWYSKjIE2lYa7uNMZtd4JOStCSrC7VNXpDUHbZOiPIDLpjr8u3j0wOgXMQxv+rmr0g/fjig0HrlSUFPTLQY1giUg7S85YGyp5uE4oKercq/HNGUZc8ty2xLGvS1s1kTp8H0+fbTQZBQJcleaK3dkaesW76AK0ydkOUl7SoGMAVxzK9astrLUvjU/d457EXJuM7+u3YmK4kGCSaiFcewRgJpINtqecMO00EPcD1HLW7zFkn7M/PA5H3o3yxvjdfMh/TbM8dz+37mHpP2LUMUeEw3t6zL1cn8p7JV0mZJHZK63fEXZzU1E/gRcCVYElPGAyh3cc2vuvkerxftRw9JBx+XbvxytPPI4kJPTPIQSAPjXODTLumOAg/dJS8o2iFpnzHmsDFmrbU2ba3dZ4xJS2rOTmJhjNkvqcsXTO01xvS57XzH7HXB2q7Msay1nS4Yy6dD0gGmwSQbQReA0otzftWS1dKt7V7QdWt7aQKuUiYMkeJPkILiI5DGdG2rb1O4eUvHtG34T7Lq/pak+yT9sEDd/do2fGB8a1v9PEm/F1g+da3+KSVu31UhhhV2aOJSQ33yerj2BVVwAVKrr0dMkva7Yx0IccxcSxsFttOdr03S2qDHIBkIugCUXjXNryr1gsxxK3WAWS0IpDF98xVu2ZxcwcIiST8NUX9+jn2pAuVT0TOdaSNuGaFmY8xGedfXoNzZrf1aJQ0aY9b79qXkDTWc7jHz6ZK0djrz0lBZCLoAlF61zK+qtgWZqy3ABCrDWYXL7Hc6YF+Yumdz7EsXKI+cWwvrlkwyNWPMpgKPH0+qYa3N2Rs21WM6qYBj7ZLUQcBVHQi6AMSjGuZXVdOCzHEGmPSuAcG8YX3TG9qXPdwwfL1zmriUT8m5rIFd1lr/m3DKlbUErBfbIqlHOTIauoBsMMQx05rc8zWpJ8wNl+zyZTlslpQKWscWla8U63QBKHfZN60nD5fmvJn5VVLp5leVUjUtyBwmwIzkvCVe780vrr8boPqEGb6XylUnaz2szHH82QKb3eMyizD3ykuQMb6GrDtGS8hjDvrb4h47oW1u6GJKXubENrfdKbIXJho9XUC1Y0hYdKppQeY4Asw4e9f4uwEi5YK
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pe.plot_corrs([impr_mass, uimpr_mass], plateau=[pcac_plateau, w_pcac_plateau], xrange=[0.5, 18.5],\n",
    "              label=['Improved pcac mass', 'Unimproved pcac mass'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Refined error analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two way of adjusting the value of S. One can either change the class variable `Obs.S_global`. The set value is then used for all following applications of the `gamma_method`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 4.79208242e-03 +/- 2.02509166e-04 +/- 2.05063968e-05 (4.226%)\n",
      "  t_int\t 1.03021214e+00 +/- 1.94552148e-01 S = 3.00\n"
     ]
    }
   ],
   "source": [
    "pe.Obs.S_global = 3.0\n",
    "pcac_plateau.gamma_method()\n",
    "pcac_plateau.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively one can call the gamma_method with the keyword argument S. This value overwrites the global value only for the current application of the `gamma_method`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 4.79208242e-03 +/- 2.04669865e-04 +/- 1.97135904e-05 (4.271%)\n",
      "  t_int\t 1.05231340e+00 +/- 1.88061498e-01 S = 2.50\n"
     ]
    }
   ],
   "source": [
    "pcac_plateau.gamma_method(S=2.5)\n",
    "pcac_plateau.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can have a look at the respective normalized autocorrelation function (rho) and the integrated autocorrelation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAGfCAYAAAD1WR7GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA00ElEQVR4nO3dfWxV953n8c/PBmwDti8m4SEmNBhIG0IoMXloZzrN7BaqSFM0mlkTtEpXWqmq2VQ7aip1cNM/hslqFdbMdLdVq82E7KorNVkNgd2qk1mJGejDjDp9SIJLiUvaEkxKoA40Mdc2D+bB/u4f51znPvte+957zrnn/ZKuju/53XPv79rXxx//no4zMwEAACB4DUFXAAAAAB6CGQAAQEgQzAAAAEKCYAYAABASBDMAAICQIJgBAACEBMEMAAAgJAhmAAAAIUEwAwAACIl5QVcAQHw55xKSnpfULalL0oCk19Ie0iHpgJkdyjquW9JTkhKSnssuDwPnXJekPknJtN3PSdoq6aiZDVXodbol7ZT3veiSNGBmfSUee0zSXr8+Sedcj6SdZrajEnUDUD7HJZkABM05t1XSEUnbzOxoVtkRSR1mtiVrf5ek05J2mdn+mlW2BH7dDkr6hJkl0/Y/J+kxSVsqEcz8ULbVzPal7cv7/SpwfPYfgCF5P4OKhEYA5aMrE0DY7ZDU7YeaaX54CGuA6Je0Nz2USZKZ7arw6+xKD2WpffK+X70lHH/If3yfvEC2llAGBItgBiDU0sLN1iDrUaZuZXZhpqtk695jzrn+9B1pwWpbCccPmdl+M9uX3VIJIBgEMwCh5ncLSt74s6hIymuFyudABV9nSN7YsnwK7QcQYgz+BxBa/uSAg/ICyGeLPK5H3kSBhKQHJX02uxvRD3i75I1Lk6S1pQ6Sn4W9kg46507LC2hHU/Uxs4oFzHzjyNKC7JESniLhnNut91v3qvk9AVACghmAMNnmh7EOeV1x3ZL6Zxjcv01SX6oLzzl3UN4Yr+nxXP4g+efTg4xzrssPTluyQ9xcmdkh51yfX4+D/usNyZtBmjEmzB8715X7LEWfv1g35S55QauULtOu9HFvzrle59yxUiYOAKgOZmUCCFyhWZl+a86DhZZv8IPVQHq5f8wuM1ub9bicgOeHuJEqDMpPPX9C3ti4nf424de3KsHHby07Jm826Kxa5vyZmjvCuAQJEAeMMQMQWn7rUs6MzCyvZt1Ppt/xW8u6JOUb3H5E3vIVVWFmSTM7ZGY7zGyJvG7NUmdMzsZBeaFqLt2lSZU2cQBAFRDMAITdkKS5BJlUN+FIgfKE37JVMf6Ytxx+0DwqqeItZn7r395SZ1c6554rUM8Rldm1CqByGGMGIOySktctOMuxYKnlIzqUu4RFQspYkqNStslbIyyfI/ImKEiqzBgzv/s24woJzrnuGVrOeuV9P7Lr2aH8rYsAaoBgBiDsUi1dqUs2lcXMBpxzSXljvLIHxD+owgFqLh5zzvUVCHxrlTZjcq7j2/xWr2SeMWE7Vfz7tS97BqY/Ri0hf8ICgNqjKxNAmCTy7EuNL9sqeRMF0paE6CjxOXZI6kvvsvTHnnWryDIcaY9NOOdOlzk27Hn/NdKfp1vSA5W6hJT/fLv8r3v9225/0dn3Zqj/gTzvp1/SfhabBYJDixmAQPnhINVq1O+c60pfUsJv8dohaZff8pWQNOKPqUpIeso5t9TM+vxA0iOpy+8ifM7MBszsqP8c/f4MTUlaqvKWykgt4VFKqHrJzHb5QWmnvy/hv59Kji/7rv+8+a6KkD2TNSPE+t9XpV05oEvSkbBddxSIG5bLAIASOed6CS4AqomuTAAoXSLoCgCobwQzACiBP54rStfrBBBBBDMAKM0DDIoHUG2MMQMAAAgJWswAAABCguUyyuCcc5LukDQedF0AAECktEr6rc3QVUkwK88dks4FXQkAABBJqySdL/YAgll5xiXp7bffVltbW9B1AVAvxsel48elzZul1tagawOgwsbGxnTnnXdKJfS4Ecxmoa2tjWAGoHKckxYtktraCGZAzDH4HwCCNn++tHKltwUQa7SYAUDQmpulD34w6FoACAFazAAgaFNT0pUr3hZArBHMACBoV65Ir77qbQHEGsEMAAAgJAhmAAAAIRH5wf/Oua2SdpnZjhIf35t2N2Fm+6pTMwAAgPJENpg557ol7ZSUkNRV4jG9Sgtjzrke51y/mfVVraKSJqdMr5wZ0cXxCS1rbdZDazrU2OCq+ZIAosQ5qaHB2wKINTfDJZtCzznXI+kpM9tSwmNPS9pmZkNp+y6Z2ZISX6tN0ujo6GjJC8weHhzW0y+f1PDoxPS+le3N2rN9gx7duLKk5wAAANE1Njam9vZ2SWo3s7Fij43NGDPnXEJSV3oo8yX81reKOzw4rCdeGMgIZZL0zuiEnnhhQIcHh6vxsgAAIKJiE8xUuLszWaRs1ianTE+/fFL52iNT+55++aQmp6LdYgmgAq5elV57zdsCiLU4BbOOAvtHCpU555qcc22pm6SSL2L3ypmRnJaydCZpeHRCr5wZKfUpAdSryUnp8mVvCyDW4hTMZuMpSaNpt3OlHnhxvHAom83jAABA/YtTMCvUNNVRpGyvpPa026pSX2xZa3NFHwcAAOpfnILZkDQ9CSBdIlWWzcyum9lY6iZpvNQXe2hNh1a2N6vQ5Hcnb3bmQ2sK9bACAIC4iU0wM7OkvACWk4TMbKDSr9fY4LRn+wZJyglnqft7tm9gPTMAUkuLdO+93hZArNVDMCs0cL8ra5V/SeqX1JP2mF5JVVtc9tGNK/Xsp7u1oj2zu3JFe7Oe/XQ365gB8MybJ91+u7cFEGuRXWA2beX/HnnLXeyXdMzM9vvlvZL6zGxt1nG75S2RkZC0tJxV/2ezwKzEyv8AZnDjhnThgrR8ubRgQdC1AVBh5SwwG9lgFoTZBjMAKGp8XDp2TNqyRWoteVUeABHByv8AAAARRDADAAAICYIZAABASBDMACBo8+ZJS5cyKxOAOAsAQNBaWqT77gu6FgBCgBYzAAiamXTzprcFEGsEMwAI2uXL0r/8i7cFEGsEMwAAgJAgmAEAAIQEwQwAACAkCGYAAAAhwXIZABC0xYulj31MamwMuiYAAkYwA4CgOcfisgAk0ZUJAMG7dk06ccLbAog1ghkABO3WLWlkxNsCiDWCGQAAQEgQzAAAAEKCYAYAABASBDMACFpzs7R+vbcFEGvMzwaAoM2fL3V2Bl0LACFAixkABO3mTenCBW8LINYIZgAQtIkJ6Y03vC2AWCOYAQAAhATBDAAAICQIZgAAACFBMAOAoDU2Sm1t3hZArLFcBgAEbeFCqbs76FoACAFazAAAAEKCYAYAQRsfl37wA28LINYIZgAAACFBMAMAAAgJghkAAEBIEMwAAABCguUyACBoixZJDz8sNTUFXRMAASOYAUDQGhqklpagawEgBOjKBICgTUxIb7zhbQHEGsEMAIJ286Z04YK3BRBrBDMAAICQIJgBAACEROQH/zvnetPuJsxsXxnHJCQtlbTXzJKVrx0AAEDpIh3M/IA1Hcaccz3OuX4z6ytyzG5J+1NBzDmXkNQvaVf1awwAeTQ1SXfdxXIZACLfldkn6VDqjpkdktRb+OGSpG3prWP+113VqBwAlGTBAi+YLVgQdE0ABCyywcxv6eoys6GsooRzrrvIoR1+qxkAhMOtW9LIiLcFEGuRDWYq3MqVLFImea1s/c65I865hHOObkwAwbp2TTpxwtsCiLUojzHrKLB/pEiZzOyoc26bpCOSLknakafVTZLknGuSlD7oo3WWdQUAAJhRlFvMZsU51yWpW9ISSfslHcya2ZnuKUmjabdzNakkAACIpSgHs5EC+zuKlElSv5ntM7Okme2StE3Sc35gy7ZXUnvabdVcKgwAAFBMlIPZkDQ9CSBdIlWWzZ8UkFFmZkcl7ZO0NfvxZnbdzMZSN0njc682AGRJXcS8IcqnZACVENkxZmaWdM4NyWshS2aVDZT5dKdVIMwBQNUtWiQ9/HDQtQAQAlH/96xfUk/qjj9WrC/tflf6+DE
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGJCAYAAAC5Lib1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuZUlEQVR4nO3df2wc533n8c+Xv/WDuysqps3YVuOVHSSKmjqSbF/Q9A6NyZ6L1C1ykOwCDg4oUEv1AbkkwEE85x/B/0SgD2gSBEUqNUCBq42DLF4D1DlEd5JzaJMi5x9icjnFaROLjhUldOmI3F1KMn8/98fO0svlLrmkZvaZ3Xm/gAW5M7Oc73K5w88+zzPPmHNOAAAACFeb7wIAAABaESELAAAgAoQsAACACBCyAAAAIkDIAgAAiAAhCwAAIAKELAAAgAh0+C7AFzMzSe+XNOO7FgAA0FR6Jf3KbTDZaGJDlooB66rvIgAAQFO6S9Iv19sgySFrRpJ+8YtfKJVK+a4FQIt4pzCrF167qscO3aXbUj2+ywEQskKhoLvvvluqoycsySFLkpRKpQhZAEIzqy717Nip3lRKKUIWkGgMfAeAEHV3tmv/nWl1d7b7LgWAZ4lvyQKAMKW3dWpo3+2+ywAQA7RkAUCIFpaW9evrc1pYWvZdCgDPCFkAEKLpG/P6m++/pekb875LAeAZIQsAACAChCwAAIAIELIAAAAiQMgCgJC1t5nvEgDEgG1w2Z2WZWYpSfl8Ps9kpAAAoC6FQkHpdFqS0s65wnrb0pIFAAAQASYjBYAQXbs+p3M/fluPfOQO7d7Z7bsctIjJwqwmZ+bWLO/v7VY/l2+KLUIWAIRoadlpsjCnpeVkDsVANJ5/+Yq++tLP1iz/3MP36QtDH/RQEepByAIAIOaeeGiPhvbdrjcmr+vzZ36orzx+v+7t36n+XlpL44yQBQBAzPWnelZ1C97bv1P770x7rKj55HI5ZTKZhu6TkAUAAEL17LPPSpKuXbsmSXr66ad1+vRpHT9+PJSf++qrryqbzWpkZGTd7S9cuKChoaGV+9lsVufPn1c2m72lOupFyAKAEKW2depTHx1Qalun71LQYpaWnX50NSdJ+tHVnD48kIrlnGzHjh3TkSNHNDg4uLLsyJEjt/xzh4eHV4WqI0eO6MiRIzp79mzNx+RyOV28eFGSlMlkGhauSpgni3myAAAxd+7ShJ558XVN5GdXlg2ke3Ti0X16ZP+Ax8rW2rVrl6anp1cty+VyevLJJ9cNROvJ5XJ6+OGH9dJLL610+Y2NjengwYO6fPlyzfA0OjqqwcHBULsJmScLADy5Mbeoi29N68bcou9S0CLOXZrQU8+NrQpYkvR2flZPPTemc5cmPFVWXS6XUy6XW7Usk8ms6rbbivHxcY2Pj6/cLwWr8mVxQ8gCgBDdmFvUP/z0HUIWQrG07PTMi6+rWp9TadkzL74eqylDBgcHdfDgQV24cGHV8qNHj275Z2YyGU1PT+vAgQMry0rhaqMuwBdeeEGjo6MaHR3V8PDwlmvYCkIWAAAx9cqbU2tasMo5SRP5Wb3y5lTjitrA2bNnV1quzExDQ0NrAlcYTp06pcHBwXVDVjab1aFDh3T48GEdPnxYe/fu1bFjx0KvpRYGvgMAEFOTM7UD1la2a4RMJqOLFy9qbGxMZ86cWTnD7+zZszp8+LCk+gfCDw0NVW0BGxsb04ULF1YGtddS3vIlFVvZjh07ppGRkYZM50DIAgAgpvp767tkTr3bNdKBAwdWQs7w8LCefPLJlZC11QHwJcPDw7p48eKmg1L5OK7KABYFugsBIERdHW3K3rZDXR0cXnHrHrynTwPpHtWaqMFUPMvwwXv6GllWTblcTqOjo2uWj4yMVB0QvxXHjh3TqVOnNgxYuVxOu3btWjUwPoz9b0ZsWrLMLCPpseDuXkkZScPOudwGjytvR8w4556Noj4AqEdme5f+6P47fZeBFtHeZjrx6D499dyYTFo1AL4UvE48ui9W82WdP39+pcWqXDabXQlGW+0uPH36tIaHh1e1SOVyuZqtUocOHVo1ZqsUuBrRiiXFKGRJGpF0yjk3JklmdkrSWUk1z/kMAtZKsDKzw2Y24pxr7OkDABBYWnaaW1xSd0d7rP7xoXk9sn9AX//MgTXzZN3R4HmySmOsNppl/YUXXtDQ0NCqoDU6Orrq/la6C0dHR5XL5Vamcsjlcjp//vxKPZX1VZs24uTJkxvWH6bYTEZqZuclnS8LTMcljTjnah6lzOyypCHn3HjZsmnn3K469sdkpABCN1mY1fMvX9ETD+1Zda054FYtLTudefWKvvjNS/rSp/fr8Qf2NDTIj46O6sknn9TFixdrntGXy+X0wgsv6NChQzpz5ox2796ta9euaffu3bd0SZ1S1181pRxTmqLh8uXLq9aXLsVz+fJlHTx48JamkpA2NxlpbEJWpaAlK+ucq9qSFXQvTleGMDNzkg6WWsTW+fmELAChI2QhSpd+mdcffO17+tZnP+HlAtFjY2NeLk8TJ00/47uZZSUNSlpvMotar3BunXUAADSdycKsLv0yrzcmr0uS3pi8rku/zGuy0NipG8bHxxMdsDYrTmOyJK2Mszom6Vh5N2AVtU6lmKq2zsy6JXWXLerdcpEAADTQ8y9f0Vdf+tnK/c+f+aEk6XMP36cvDH2wITXkcjkC1ibFLmQ5505LOm1m583sQIhnCz4t6URIPwsAgIZ54qE9Gtp3+5rl/b3dVbaORiaTadhZea0izmOyBiWdl7S3WouWmR2QdLHKmKxpSU8650YrlldrybrKmCwAYVpedlpYXlZnW5vaOLsQaDlNNybLzDJmdjYYzF5SClaDNR42XnpsxfJM2WNXOOfmnHOF0k3SzC0VDQBVtLWZujvaCVgA4hGyVByoPqjVY6kywdeq47KCSUrHVWX81UZnFgJAVKZvzOtvx65q+sa871IAeBaLkBWEotMV3YKPSxpzzl2QimccVszuLhUnMF2Z3SxYz0SkALxZWFrWW9duamFp2XcpADyL08D3k2ZWPg1rRtLDZfcHVQxQp0sLnHOnzex4aeZ3SbuZ7R0AAMRBbEJW0P1XMyCVzjqsspxrFQIAgNiJRXchAABAqyFkAUCIdvZ06Hc/1K+dPbHpKADgCUcBAAjR9q4O3X93xncZAGKAliwACNHswpJ+MlHQ7MKS71IAeEbIAoAQFd5d0LlLb6vw7oLvUgB4RsgCAACIACELAAAgAoQsAACACBCyACBEHe1tGkj3qKOdwyuQdEzhAAAh6tvRpT9+cI/vMgDEAB+1AAAAIkDIAoAQTRZm9eXzP9VkYdZ3KQA8I2QBAABEgJAFAAAQAUIWAABABAhZAAAAEWAKBwAIUd+OLv3Jb39AO7s5vAJJx1EAAELU0d6mzPYu32UAiAG6CwEgRPmbCzp3aUL5mwu+SwHgGSELAEI0t7ikn0zMaG5xyXcpADwjZAEAAESAkAUAABABQhYAAEAECFkAEKLt3R36V9nd2s4UDkDicRQAgBDt7O7Qx/fu9l0GgBigJQsAQjS3uKSf//oGZxcCIGQBQJjyNxf0zR/8knmyABCyAAAAokDIAgAAiAAhCwAAIAKELAAIUVubKbO9U21t5rsUAJ4xhQMAhOh9O7v1J799j+8yAMQALVkAAAARIGQBQIjemZnTX/79Zb0zM+e7FACeEbIAIETOOb07vyTnnO9SAHjWMiHLzDK+awAAACiJ1cB3MzsefPuApHHn3PAG2w9KOl92f1zSkHNuPLoqAQAANhabkGVmI+WhyszOmtlZ59yRdR6WkXQw+D5HuAIAAHERi5AVdPUNmlnGOZcLFp+UdNHMshuEp/GyxwCAV5ntXXr8gbuV2d7luxQAnsVpTFY2uJWMly0HgKbQ1dGm92e2qasjTodXAD7EoiUraInaVbG4FK426gJ8zMymgu8f2GgcFwBEaWZ2QWNXcjqwJ6Penk7f5QDwKBYhq4Zjki5s1FWoYnfhmCSZWZ+ZnXLOHavc0My6JXWXLeoNtVoAkPTu/JLG3prWh+/oJWQBCRfL9mwzOyBpUNJ6g97
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pcac_plateau.plot_rho()\n",
    "pcac_plateau.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Critical slowing down"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pyerrors` also supports the critical slowing down analysis of arXiv:1009.5228"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 4.79208242e-03 +/- 2.28649024e-04 +/- 1.67571716e-05 (4.771%)\n",
      "  t_int\t 1.31333644e+00 +/- 5.19554793e-01 tau_exp = 10.00,  N_sigma = 1\n"
     ]
    }
   ],
   "source": [
    "pcac_plateau.gamma_method(tau_exp=10, N_sigma=1)\n",
    "pcac_plateau.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The attached tail, which takes into account long range autocorrelations, is shown in the plots for rho and tauint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAGfCAYAAAD1WR7GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6bklEQVR4nO3dfXAc933n+c9v8EziYQhSfBAfJIKibFGULIF6TGzrLia9upxVcbKglT37rlLlE7lKdstWVUJavqsoutqVFsr5Epe9K4vaPded5UoocsvRarORI/opJz/ogTBNwZRtEaAtkaYIieAMABIAQeB7f3QP3JgnzAxmMN2D96sK1Zz+TTd+mAGnP/g9tTMzAQAAoPpi1a4AAAAAPAQzAACAkCCYAQAAhATBDAAAICQIZgAAACFBMAMAAAgJghkAAEBIEMwAAABCor7aFQAQXs65uKSnJXVL6pLUJ+m1wFM6JR00s8Npx3VLelhSXNJT6eWoHf7vSK+Z7c1XLmlAUkJS3MyeWKz6AVHjWPkfwHycczslvShpl5kdSSt7UVKnme1I298l72K818wOLFpli+Cc6zKzwWrXI6r8AH5I0qCZ7crxnAF5vzeD/uMeSfeb2e7FqykQHXRlAlio3ZK6nXNPBXf6F+LQhh4/VHRXux5R5JyLO+cOSbpfXitYruftk9QXDL9+6+lO//UHkIZgBmBBzCzh/3NnNetRgvurXYGoMrOEme02s/3KH77vl/Rqlv2vKXq/L8CiIJgBWBC/y1Lyxp9FgnNuj6R9ZThPl3OuO/AapFqTumgRkuS1SCay7E9Iytr1CSx1DP4HUDJ/YPchea0mD+R5Xo+8iQJxSbdLeiDQ0pZ6TpekvfLGpUnSFr9Fptx17pGUGg+31zmXCgj7g3UK1Fn+8w8Fx9f59d0vaY+/TQ1of9jfNyxpS4F12qnfhJgtks6nBsj7IXJvoHyHmQ0653rlhctBSU/JC8a98iZpPOBvJWmlvNd9f/prXkn+74bkvQ7phkU3MpAVwQxAMXb5F9xOeS0e3fJm5OUb3L9LXihIDf4+JC9AzM7i81uXng5OIPBbnQbkBZFEuX4Af4zTYT/wZJ0xmm2AunNuwDm328z6/PMMygt2e9LOv985dz748+Xjf6+9wcHzzrmnnHO9Zrbff20P+K9bcLLCi/JmOAa/zw7nnPnfe3fqdfO/xynn3ObAvqf0m/BWkFwD/HPoXGA5sCQRzAAU48VAq9EBf3D3Lkl5Z12mzXx8VZmhJRXW5hzjnEu1AhUUcsosfQxUn1+P9LpkG2OVKOL7PK3M1sZD8oJXsMXwAXnhap/fmrYrxxIVCXmte7N1MLPDzrmH/e+1299Xjdc0hVAG5EAwA1AyM3vCb0l6Ks+FPn3wdyL4wG8t65J0RJleVBWCWaBVrVteQEuoyNalQvjnj0vq8lu1UjolzWnJM7OEc+4jko4657ZobmhLl6378IjKMK6uCKk65Aph2eoILHkEMwALNShvTFWp4SkVeHJdqOPOuXgVxkd9S9KR1Dg359wOlb+lJ/WzH86ynlpGK6SZ9TnnDkjaWcLrMSB5P9tivJZ+kJS84JkurhAvpQJUE8EMwEIlpAVd8FMX6E5ldgHGpTlLclSMc25PYKzcIUmv5Zp8UMDCtPECv23qHF0qIKj4gXFA0m3ztFJms0VSYhHHmEle9+/KLPs7JR0s8lzAkkAwA7BQqZau1C2biuK3AiXkdRmmtxLdrrQuvTJKpD2OB/69U/5YrICuwDHBuqafRypwNqb/sw/KG6eXfkeF1FiyoIf9yQWHJQ04517McburbC17OyU9G/jei9E9fFDZl8XoVp5ZvMBSxjpmAIoRz7IvteL/Tslb+iGwrle2gJDtHLsl7Q8ssRBcmX/eC7i/dthA+gzJebwmL/ilvlcwVCZSZX55qkUrVb/gz/CaAkHM/xlukzdurJAWqb2S9qSthRYMgal9h+S3MqVmhEp6OviaBcwJlf7rElf+cWmVcEBpr4Nfl8Op2a0A5uJemQDySltHa1DeEhNPpD2nx3/OIXkB4Ii89bx65AWMA35LT6+/r0veRfup1AXaD0fBdcxWSnq8kG5MP5yckjcmrKB7MPphIbX+1/ngz5RWNiBp2J/ZeNT/2Z4KLP8Rlzfb8VX/Zx2WF0h7/X/vlTeJIYOZucD36w2cQ6lu1bTXbH9gfbND+s3r+2yqBcw5d0G/CWDD/nFbyt1C5tery6+D5LVsDirtPfNfn4flvY5xSSsrsT4dUCsIZgBqRto4sSXJD2YP5OjiBBBydGUCqCXxalcAABaCYAagJmQZJ7ZUxcUCrkBkEcwA1IrbgveyXGr8m6kf8h+mxvMBiBjGmAEAAIQELWYAAAAhQTADAAAICVb+L4Lzbvx2taTRatcFAABESpukX9s8Y8gIZsW5WtLpalcCAABE0gZJZ/I9gWBWnFFJevvtt9Xe3l7tumA+o6PSsWPSLbdIbW3Vrg0AYIkaGRnRxo0bpQJ63AhmJWhvbyeYRYFz0vLlUns7wQwAEAkM/kftamiQ1q3ztgAARAAtZqhdzc3S+95X7VoAAFAwWsxQu2ZmpIsXvS0AABFAMEPtunhRevVVbwsAQAQQzAAAAEKCYAYAABASkR/875zbKWmvme0u8Pl7Ag/jZvZEZWoGAABQnMgGM+dct6T7JcUldRV4zB4Fwphzrsc512tm+8tVr+kZ0yunhjU0OqHVbc26Y3On6mKuXKdHMZyTYjFvCwBABLh5btkUes65HkkPm9mOAp47IGmXmQ0G9l0wsxUFfq92SclkMpl1gdkX+s/q0edP6GxyYnbfuo5mPXLfNt27fV0h3wIAANSYkZERdXR0SFKHmY3ke+6SGWPmnItL6gqGMl/cb31bkBf6z+rBZ/rmhDJJeic5oQef6dML/WcX+i0AAECNWzLBTLm7OxN5ygoyPWN69PkTytb2mNr36PMnND0T7dbJyLl0SXrtNW8LAEAELKVg1plj/3CuMudck3OuPfUlKesNF185NZzRUhZkks4mJ/TKqeEiq4wFmZ6Wxsa8LQAAEbCUglkpHpaUDHydzvakodHcoayU5wEAgKVpKQWzXM1VnXnKHpfUEfjakO1Jq9uaC6pAoc8DAABL01IKZoPS7CSAoHiqLJ2ZTZrZSOpL0mi2592xuVPrOpqVa1EGJ2925h2bc/WmAgAALKFgZmYJeQEsIx2ZWd9Czl0Xc3rkvm2SlBHOUo8fuW8b65kttpYW6cYbvS0AABFQC8Es18D9rrRV/iWpV1JP4Dl7JJVlcdl7t6/Tk5/q1tqOud2Vazua9eSnulnHrBrq66WrrvK2AABEQGQXmA2s/N8jb7mLA5KOmtkBv3yPpP1mtiXtuH3ylsiIS1pZzKr/8y0wK7Hyf6hcviydOyetWSM1Nla7NgCAJaqYBWYjG8yqoZBghhAZHZWOHpV27JDasq50AgBAxbHyPwAAQAQRzAAAAEKCYAYAABASBDPUrvp6aeVKZmUCACKDKxZqV0uLdNNN1a4FAAAFo8UMtctMmprytgAARADBDLVrbEz6/ve9LQAAEUAwAwAACAmCGQAAQEgQzAAAAEKCYAYAABASLJeB2tXaKn3wg1JdXbVrAgBAQQhmqF3OsbgsACBS6MpE7Rofl44f97YAAEQAwQy168oVaXjY2wIAEAEEMwAAgJAgmAEAAIQEwQwAACAkCGaoXc3N0tat3hYAgAhgLQHUroYGaf36atcCAICC0WKG2jU1JZ07520BAIgAghlq18SE9MYb3hYAgAggmAEAAIQEwQwAACAkCGYAAAAhQTBD7aqrk9rbvS0AABHAchmoXcuWSd3d1a4FAAAFo8UMAAAgJAhmqF2jo9J3v+ttAQCIAIIZAABASBDMAAAAQoJgBgAAEBIEMwAAgJBguQzUruXLpTvvlJqaql0TAAAKQjBD7YrFpJaWatcCAICC0ZWJ2jUxIb3xhrcFACACCGaoXVNT0rlz3hYAgAggmAEAAIQEwQwAACAkIj/43zm3J/AwbmZPFHFMXNJKSY+bWaL8tQMAAChcpIOZH7Bmw5hzrsc512tm+/M
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGJCAYAAAC5Lib1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9W0lEQVR4nO3df3xc9X3v+fd3RtLot0ayLdtgGywbA8b8ksCEhKQJyA1pQrNpbdNeSHv7I3bIvSkk26KQTZdLtxtWvre3SbMJsdvuvS2wrUEt3ZLcuLUhTZqblB8WlCoQwJaDsWMssH5LntFo9N0/ZkaMRjPWyDpH5xzN6/l4zEOac2bOfDQzmvOe7/d7vsdYawUAAABnhbwuAAAAYCkiZAEAALiAkAUAAOACQhYAAIALCFkAAAAuIGQBAAC4gJAFAADggjKvC/CKMcZIukDSiNe1AACAQKmT9DM7x2SjJRuylApYJ7wuAgAABNIaSSfPdYNSDlkjkvTmm2+qvr7e61owh7eHY3rs+RPaed0araiv9LocAECJGh4e1tq1a6UiesJKOWRJkurr6wlZARBThSpralVXX696QhYAIAAY+I5AiJSHteXCBkXKw16XAgBAUUq+JQvB0FBVrm2bV3pdBgAARaMlC4GQSE7pndG4Eskpr0sBAKAohCwEwsDYhB7+0RsaGJvwuhQAAIpCyAIAAHABIQsAAMAFhCwAAAAXELIQGOGQ8boEAACKZuY47c6SZYyplzQ0NDTEZKQAAKAow8PDamhokKQGa+3wuW5LSxYAAIALmIwUgXBmNK4DP35Lt16xSstqI16XAwAL0jccU99IfNby5rqImjl12JJByEIgJKes+objSk6VZvc2gKXl0WeO66tPvT5r+d23XKLPbdvkQUVwAyELAIBFdscN67Rt80od6RvVPftf1Fduv0Ybm2vVXEdL/VJCyAIAYJE111fO6Bbc2FyrLRc2eFgR3EDIAgAAgXXo0CHt3btXjz/++Kx1+/btm/59cHBQ995772KWxtGFCIb6qnJ99KrVqq8q97oUAHBEcsrqpRODkqSXTgwy5nSeuru71dHRoccff1y9vb2z1u/bt0+Dg4PatWuXdu3apZaWFnV0dCxqjcyTxTxZAIBFdqDnlB548mWdGopNL1vdUKn7b9usW7es9rCy4Onq6tKDDz6ow4cPz1i+YcMGHTx4UC0tLdPLGhsbNTAwsKDHY54sLDlj8UkdfmNAY/FJr0sBgAU50HNKdz3SPSNgSdJbQzHd9Ui3DvSc8qiypWNwcFC9vb0zAlZmeXd396LVQchCIIzFJ/X9194mZAEItOSU1QNPvqx8fUiZZQ88+TJdhwuUr/tQkqLRaMF1biBkAQCwSJ491j+rBSublXRqKKZnj/UvXlFzOHTokDZs2KANGzaosbFx+vcNGzZocHDQ6/Ly6u/P//w1NTUVXOcGji4EAGCR9I0UDljnczu39fb2qre3V0ePHpUk7dmzx5Uj9Hbs2FHU7bZt26Zdu3Y5/vhuIWQBALBImuuKO2VOsbdzW1NT03So6e3tVTQadeVx8k2/sBBNTU15l/f39xdc5wa6CxEIFWUhtayoUUUZb1kAwbV1fZNWN1TKFFhvlDrKcOv6xQsC55Idqrq6unTdddd5V8w8ZAa853ZnDg4OzhoM7ybftGQZY6KSdqavbpAUldRhrR2c437Z7YZRa+0eN+qDt6LVFfr4NRd6XQYALEg4ZHT/bZt11yPdMtKMAfCZ4HX/bZsVDhWKYd7Zv3+/a5N5Ot1dGI1G1dLSov7+/lmtb62tredT4nnxTciS1Clpr7W2W5KMMXslPS5pW6E7pAPWdLAyxmw3xnRaaxd3tjG4LjllFZ9MKlIW9uWHDwAU69Ytq/XQna2z5slatcjzZHV3d2v//v3q7Oyc87aFBrh3dHRMD4Dfvn27ent7tXv3bnV0dGjnzp1qa2tTZ2enotGoduzYoc7OTjU1NengwYPq6OiYblVaSHdhoYHsHR0d6urqmg6G+/btK+pvdZJvJiM1xhyUdDArMN0rqdNaW3CPaow5KmmbtbY3a9mAtbaxiMdjMtIA6RuO6dFnjuuOG9bNON8XAARVcspq/3PH9cUnevTlT2zR7devW9QvkV1dXfrUpz6lw4cPz9mF1tXVpeeee25GSNm3b5+OHj06vWzbtm06ePDgdNDavXu32tvbp1uS2tra9NRTTykajWpwcFBtbW3TA+rPRyYkdnV1qbe3V7t27VJbW9uMlq49e/ZMP96ZM2ccCVnzmYzUNy1Z1trcFqsNkg4Vun26e7ElO2ClRY0xrZkWMQAA/CgcMrpqTVSSdNWa6KK30m/fvr3o8Unbt2/X9u3bZyw7fPiwotGoDh1K7aoz3XAtLS3atm2b9u/fP+s+mcAVjUbV39+vwcHB8x5M39raqtbW1nMGp8U+V2Eu34SsbMaYFkntOkdXoaRC74zB9DpCFgDAl/qGY+obietI36gkTf9srossamt9b2/vrCBUrLa2Ng0ODqq9vV3SzMHmLS0tOnPmjLq6umZsPxOqBgcH1dTU5NrRin7hu5CVHme1W9LuPK1U2QodetGfb50xJiIpkrWo7ryLBABgAR595ri++tTr09fv2f+iJOnuWy7R57ZtWpQaFnqk3a5du7Rnzx7t27dPTU1NamlpUVdXlzo6OnT48GG1trZqw4YN2rt373QX3mOPPaampiY999xzOnjwoFN/im/5ZkxWrtwxWnnWt6fXm5zlR5Uay7UvZ/l/knR/7nYYkxUMjMkCsJRkWrJyLXZL1mJqa2ubdRLnIJrPmCw/h6x2SQclbcjXomWMaZV0OE/IGpD0KWttV87yfC1ZJwhZwTA1ZZWYmlJ5KKQQRxcCQKB0d3frlltuUWdnZ6BmbM8ncAPf04PY/1SpcDSYXpwJVu2S9uW5W2/mvjlzaUWz7jvNWhuXNP21wRh21EESChlFQmGvywAAnIfW1lYNDAx4Xcai88v02ZmB7tljqaLpn3nHZaWDVa/yjL/iyMKlZ2BsQn/bfUIDYxNelwIAQFF8EbLSoWhfTrfg7ZK6rbWHpNQRhzmzu0upCUynD1tIr2ci0iUokZzSG2fGlUhOeV0KAABF8UV3YdqDxpjsyS6ikm7Jut6uVICa7jq01u4zxtybmfld0jJmewcAAH7gm5CV7v4rGJDSRwvOGpvFuQoBAIAf+aK7EAAAYKkhZCEQaivL9KHLmlVb6ZvGVwAAzok9FgKhuqJM16yNel0GAABFoyULgRBLJPXKqWHFEkmvSwEAoCiELATC8NmEDvS8peGzCa9LAQCgKIQsAAAAFxCyAAAAXEDIAgAAcAEhC4FQFg5pdUOlysK8ZQEAwcAUDgiEppoK/crWdV6XAQBA0WgWAAAAcAEhC4HQNxzTHx98TX3DMa9LAQCgKIQsAAAAFxCyAAAAXEDIAgAAcAEhCwAAwAVM4YBAaKqp0G+872LVRnjLAgCCgT0WAqEsHFK0usLrMgAAKBrdhQiEofGEDvSc0tB4wutSAAAoCiELgRCfTOqVUyOKTya9LgUAgKIQsgAAAFxAyAIAAHABIQsAAMAFhCwEQnWkTO9pWaZqpnAAAAQEeywEQm2kTDduWOZ1GQAAFI2WLARCfDKpn74zxtGFAIDAIGQhEIbGE3rihZPMkwUACAxCFgAAgAsIWQAAAC4gZAEAALiAkIVACIWMotXlCoWM16UAAFAUpnBAICyvjeg33rfe6zIAACgaLVkAAAAuIGQhEN4eieub3zuqt0fiXpcCAEBRCFkIBGutzk4kZa31uhQAAIqyZEKWMSbqdQ0AAAAZvhr4boy5N/3r9ZJ6rbUdc9y+XdLBrOu9krZZa3vdqxIAAGBuvglZxpjO7FBljHncGPO4tXbHOe4WldSW/n2QcAUAAPzCFyEr3dXXboyJWmsH04sflHTYGNMyR3jqzboPlqhodYVuv36totUVXpcCAEBR/DQmqyV9yejNWo4SV1EW0gXRKlWU+ektCwBAYb5oyUq3RDXmLM6Eq7m6AHcaY/rTv18/1zguBNNILKHu44NqXRdVXWW51+UAADAnX4S
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pcac_plateau.plot_rho()\n",
    "pcac_plateau.plot_tauint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additional information on the ensembles and replicas can be printed with print level 2 (In this case there is only one ensemble with one replicum.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result\t 4.79208242e-03 +/- 2.28649024e-04 +/- 1.67571716e-05 (4.771%)\n",
      "  t_int\t 1.31333644e+00 +/- 5.19554793e-01 tau_exp = 10.00,  N_sigma = 1\n",
      "1024 samples in 1 ensembles:\n",
      " : ['B1k2r2']\n"
     ]
    }
   ],
   "source": [
    "pcac_plateau.print(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Monte Carlo history of the observable can be accessed with `plot_history` to identify possible outliers or have a look at the shape of the distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAFyCAYAAAC5qt3eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABj/UlEQVR4nO29f5Bdx3Xf+W0OOaOBOIPBwCYgGATMoRmoNlpIS5ohy7AoxQJTMrmKGWMobqVSW17vGjQjW1RKDhBGa6tcyUYC7JXNpWmasGuz5U25zAiwqFCWVEuYikjBRcUSTY2pXaFEDoMfxg9ZGAwGJAYz4LD3j3n90NPT3bf73r6/3vt+qqaA9+7re/v2z9PnnD4tpJQghBBCCCH1cE3dGSCEEEII6WcojBFCCCGE1AiFMUIIIYSQGqEwRgghhBBSIxTGCCGEEEJqhMIYIYQQQkiNUBgjhBBCCKmRa+vOQBMRQggAmwBcrDsvhBBCCGkVIwBOyYhArhTG7GwCcLLuTBBCCCGklWwG8LehP65EGBNC7NY+jkkp9xdNI4TY0/nv7QCmpZR7iz5T4yIAnDhxAqOjoxHJCCGEENKvzM3N4cYbbwQiLWui7OOQOkJRVxgSQkwCuN0UnmLSCCH26emFEJ8HACnl/XmfaTx/FMCFCxcuUBgjhBBCSBBzc3NYu3YtAKyVUs6FpqtCGHsNwN1Symntu/NSynV50gghxgD8BYAPSSlnO9duBfBtADdLKafzPNN4PoUxQgghhESRVxgrdTdlR3Ca0IWiDmMdASpvmonOn0L9diLPMwkhhBBC6qLs0BYTju9nPde8aaSUs1LKdVLKlyxppvM8UwgxJIQYVX9Y3glBCCGEEFI6ZQtj447vZzzX8qR5EMDhjjYsT/pHAFzQ/riTkhBCCCGV0Pqgrx3T404A9xe4zWcArNX+NifIGiGEEEJIJmWHtphxfD/uuRabZh+A25Qzf55nSikXACyoz8sxXwkhhBBCyqdszdg00HXK1xnDVaf73GmEEE8CeFATxPI+kxBCCCGkFkoVxjpCktWPy3DAj07TiSW2T+2aFEJMCCFuzfNMQgghhJC6qMJnbB+ASfWhI0TpAVsnjGj5IWkmsazpmhBC7Ox83ourmi9vekIIIYSQplB60Fege3TRLJYFqPVG9PzdAPZKKW8OSdMxP563PUdK2XX28j0zIL8M+koIIYSQKBobgb+NUBgjhDSd+cUlHD17Eds2jGB4cKDu7BBCkF8Yq+SgcEIIIemYX1zCPY8+j9fPXcJN69fgyw/fRYGMkBbT+jhjhBDSbxw9exGvn7sEAHj93CUcPXux5hwRQopAYYwQ0ljmF5fw8olZzC8u1Z2VRrFtwwhuWr8GAHDT+jXYtoEnuBHSZugzZoE+Y4TUD01xfugzRkjzyOszRs0YIaSR0BTnZ3hwAO+7cYyCGCE9AIUxQkgjoSmOENIv0ExpgWZKQpoBTXGEFIN9qFoY2oIQ0nMoUxwhJB76XbYHmikJIYQ0Eu6mLQb9LtsDNWOEEEIaB7U6xVF+l6oM6XfZXCiMEUIIaRw2rQ5N1nEMDw7gyw/fRZ+xFkAzJSGk8dBc1X9wN20aGAKlHXA3pQXupiSkOdBc1b9wJyBpGwz6SgjpSeiE3L9Qq0P6BQpjhJBG00ZzFc2qhJAYaKa0QDMlIc2iTeYqmlVJk2hT3+kFaKYkhPQsbTJXVWFWTa15oyavN1ELg/seP4J7Hn2e9dtgGNqCEEISUnZsp9SaN2ryeheGB2kPFMYIISQhZcd2sk2w2zaM5H4eJ+zehUFf2wOFMUIISUyZZ2qaE+yWdWsKabY4YfcuDPraHujAb4EO/ISQJqM7ZR89exH3PX6ke+3pj+2IFgTp5E1IGujATwghfYK+oSFF6I82bZAgpBehZswCNWOEkDZBzRYhzYCaMQdCiLG680AIIWVCzRYh7aYSB34hxG7t45iUcn+KNEKInQAelFLeb/n+We3zNIC7pZTT0ZknhBBCCCmR0oWxjlDVFaaEEJNCiH1Syr150wghbgXwAIAxABOWW4wBuK3z/1kKYYQQQghpKqX7jAkhXoOhlRJCnJdSriuaRggxCeARKeVtlu8PSylnc+aZPmOEkFZAfzFCmkMjfcY6/loTFs3UWEe7lSQNIYT0I67jbvIeb8RjkQiph7LNlDYTIgDMdq69lCiNjY8KIWY6/7/dZxYlhJA24orGnycILI9FIqQ+yt5NOe74fsZzLU8ak2kA35JSHpRSHgTwmhDiSdePhRBDQohR9QeAIagJIY3HFmMs70HlVRxwTgix05PHIUkpTe3ZYQBPCiH2OvzIHgHw6dIzRgjpC6ry47Idd5P3eCMei0RIfZQtjM04vh/3XMuTxouUcloIAbjNnJ8B8Dnt8wiAk3meRQjpb6o295nnYOY9j5DnGBJSH2WbKacBa+DVMXUtUZouQogxIcR5IcSE/p0vjZRyQUo5p/4AUD9PCMlFE8x9eYPAMnhsdXCzBNEpVTMmpZztBFwdx7IDvn7N6oifJ42Fbxm7MSci0xNCSC5o7iNZcLMEManiOKR9ACbVh05A173a5wkj2n5mGo1VDv0dn7Bnja8fcaQnhJCkKHPf0x/bwUmWWGmC9pQ0i0oOChdC7MGylmsMwHo9zIQStKSUN0ekURH4J7Gs9ToA4NtSygNGegC42bwWkN/GBH1lQEfSr7Dtswx6FWrGepe8QV8rEcbaRlOEMXZY0q+w7TezDOoQDntVIO3V9+p3GhmBnxSDqmzSr7DtN68MXNH+e+2ZVcHNEkSHwliDsQV0JKQfYNtvXhnUIRxOnZxtlEBKSFn0ZNDXXoFxf0i/wrbfvDKoepfo/OIS9hz8Tvfz1gYIpISUBX3GLDTFZ4wQQppElX5OL5+YxX2PH+l+fmr3nbhjYn2pzySkKPQZI62DQQ8JCesHTekrKf2cst7JNNNu3zxW+JmENBVqxixQM1Y+TdwpRkjVhPSDXuwroe/EHYekbVAzRioh1Qq9aTvFSLNpimYoNSH9oBf7Sug79cOOw15t2/3I/OISpk7O5kpLB34STMoVuu4MvGl0CFvWrUmcW9IruNpdL2hNQpzie/F4pV58pzz0otazqZQ9Xqi6fO3UD3OlpzBGgrGtZt9341iuew0PDuDQQztw72Mv4NSFy9j1xBEORMSKrd1t2zDSE5NYyI7Jpu2qTEEvvFOKyT3lmErcVCH06nWZB5opSTCp4x4dP38Jpy9cBtA75heSTaxZxtbuesl0F2KK60VzXZF3qtu0pwej/fDvfh3fnD63Ki8hedyybg3etfYdAJoRS65XqWK80MepPFAzRoJJvZqlqcJOL5jfXORZodraHdtO/9IE054+uR+bmccDB15cZUIP2Zix64kjOH3hMjaNDuHQQzt6rr83hSrGCzVOffvVv8X7fzc+PYUxEoVazaa6V9tNFYpUAlQTJpoyyWuWMdtdL7UdEkcTTHv65K7Q8xKSR/03p+YWcPz8JYxfP5g0nyHjUi8v/hR5xwuzbLLKanhwIHcIFgpjpFZSCnd1kVKAasJEUyYpV6i90HZIPEXbUArhQ03uUydnsefQFI4ZedHzuHV8GAtXljC/uLTieWVra/o1bIqL2PHCLJtDD+3ArieOlFZWFMYIKUhKAarXzW9labT6YXVfNm0pwyJtKKXwMTw4gDsm1uOrlryYwtoDB17EptEhfOnjd3W1X2Vrd2O1c724+FPkadtm2Rz+3tlSy4rCGCEFSa3t6XXzW2qNVj+t7oFyhKa2lWHeNlSG8OHKy/DgAIauG8AxzRR572Mv4LlPfnCF0FaW8NOvYVNM9LZtCsQ+zLLZ+e4NpZYVhTFCCpJagFIO6k0RyJquMemX1T1QntDUL2VYtfCxbcMINo0O4dTcAgDg9IXL3bItu1/1a9gUE9M3zxSIXdjKpsyyojBGSAJSrnDr0lLYJoc2aEz6YXWvKEto6pcyrFr4GB4cwJc+fhfufewFnL5wuVu2VfWrkHGp130vfQJxFraNQ2UtlCmMVUDsjgzS39ShpXBNDm3QmPTD6h5YrqOFK0vYun7NKofxovRLGQL
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pcac_plateau.plot_history()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If everything is satisfactory, dump the `Obs` in a pickle file for future use. The `Obs` `pcac_plateau` conatains all relevant information for any follow up analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcac_plateau.dump('B1k2_pcac_plateau')"
   ]
  },
  {
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}