{
 "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": [
    "In this example we look at the analysis of the current quark mass (PCAC mass) on a test gauge field ensemble with fixed Schrödinger functional boundary conditions in the temporal direction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can load data from preprocessed files which contains lists of `pyerror` `Obs` and convert them to `Corr` objects as explained in the previous example. We use the parameter `padding` to keep track of the fixed boundary conditions at both temporal ends of the lattice. This allows us to specify absolut temporal positions without having to keep track of any shifts in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been written using pyerrors 2.0.0.\n",
      "Format version 0.1\n",
      "Written by fjosw on 2022-01-06 11:27:27 +0100 on host XPS139305, Linux-5.11.0-44-generic-x86_64-with-glibc2.29\n",
      "\n",
      "Description:  SF correlation function f_A on a test ensemble\n",
      "Data has been written using pyerrors 2.0.0.\n",
      "Format version 0.1\n",
      "Written by fjosw on 2022-01-06 11:27:34 +0100 on host XPS139305, Linux-5.11.0-44-generic-x86_64-with-glibc2.29\n",
      "\n",
      "Description:  SF correlation function f_P on a test ensemble\n"
     ]
    }
   ],
   "source": [
    "p_obs_names = [r'f_A', r'f_P']\n",
    "\n",
    "p_obs = {}\n",
    "for i, item in enumerate(p_obs_names):\n",
    "    tmp_data = pe.input.json.load_json(\"./data/\" + item)\n",
    "    p_obs[item] = pe.Corr(tmp_data, padding=[1, 1])\n",
    "    p_obs[item].tag = item"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the method `Corr.show` to have a quick look at the data we just read in "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 16164 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_obs['f_A'].show(comp=p_obs['f_P'], y_range=[-0.8, 8])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constructing the PCAC mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PCAC mass is defined as\n",
    "\\begin{align*}\n",
    "am(x_0)=\\frac{a\\tilde{\\partial}_0 f_\\mathrm{A}(x_0)+a^2c_\\mathrm{A}\\partial_0^{\\ast}\\partial_0^{}f_\\mathrm{P}(x_0)}{2f_\\mathrm{P}(x_0)}+\\mathrm{O}(a^2)\\,.\n",
    "\\end{align*}\n",
    "\n",
    "We now need to obtain the first derivative of f_A and the second derivative of f_P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_deriv_fA = p_obs['f_A'].deriv()\n",
    "first_deriv_fA.tag = r\"First derivative of f_A\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "second_deriv_fP = p_obs['f_P'].second_deriv()\n",
    "second_deriv_fP.tag = r\"Second derivative of f_P\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use these to obtain the unimproved PCAC mass:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "am_pcac = first_deriv_fA / 2 / p_obs['f_P']\n",
    "am_pcac.tag = \"Unimproved PCAC mass\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And with the inclusion of the improvement coefficient $c_\\mathrm{A}$ also the $\\mathrm{O}(a)$ improved PCAC mass:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "cA = -0.03888694628624465\n",
    "am_pcac_impr = (first_deriv_fA + cA * second_deriv_fP) / 2 / p_obs['f_P']\n",
    "am_pcac_impr.tag = \"Improved PCAC mass\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take a look at the time dependence of the PCAC mass with the method `Corr.show`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 26096 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "am_pcac_impr.show(comp=am_pcac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plateau values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now construct a plateau as a derived observable from the masses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit with 1 parameters\n",
      "Method: Levenberg-Marquardt\n",
      "`ftol` termination condition is satisfied.\n",
      "chisquare/d.o.f.: 0.2704765091136813\n",
      "Result\t 5.03431904e-03 +/- 5.38835422e-04 +/- 8.24919899e-05 (10.703%)\n",
      " t_int\t 5.15384615e-01 +/- 1.25000000e-01 S = 2.00\n",
      "64 samples in 1 ensemble:\n",
      "  · Ensemble 'test_ensemble' : 64 configurations (from 1 to 64)\n"
     ]
    }
   ],
   "source": [
    "pcac_plateau = am_pcac_impr.plateau([7, 16])  # We manually specify the plateau range here\n",
    "pcac_plateau.gamma_method()\n",
    "pcac_plateau.details()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the data with the two plateaus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 31548 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "am_pcac_impr.show(comp=am_pcac, plateau=pcac_plateau)"
   ]
  },
  {
   "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": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 15068 characters]\n",
      "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 we can save the `Obs` in a file for future use. The `Obs` `pcac_plateau` conatains all relevant information for any follow up analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcac_plateau.tag = \"O(a) improved PCAC mass extracted on the test ensemble\"\n",
    "pcac_plateau.dump(\"pcac_plateau_test_ensemble\", datatype=\"json.gz\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
