{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "adc5eb9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pyerrors as pe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a921a326",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('./base_style.mplstyle')\n",
    "plt.rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd0fac2c",
   "metadata": {},
   "source": [
    "A correlator can not only contain a list of Obs, but also a list of matrices of obs. \n",
    "This is useful, if there are multiple sources and sinks used. In our example, the sources have a different Gaussian smearing applied. \n",
    "\n",
    "We will load such a correlator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c49f7c4a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been written using pyerrors 2.0.0+dev.\n",
      "Format version 0.1\n",
      "Written by jan on 2022-01-27 10:49:51 +0100 on host endwings, Linux-5.13.0-27-generic-x86_64-with-glibc2.10\n"
     ]
    }
   ],
   "source": [
    "matrix_V1V1 = pe.input.json.load_json(\"./data/matrix_correlator_V1V1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a56dca5f",
   "metadata": {},
   "source": [
    "We can print the content at timeslice 0. As we can see, it is a matrix of Obs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dd28fe5b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Obs[0.95214(67)] Obs[0.01240(11)] Obs[0.005965(72)] Obs[0.002719(40)]]\n",
      " [Obs[0.01241(12)] Obs[0.004389(60)] Obs[0.002672(41)] Obs[0.001432(25)]]\n",
      " [Obs[0.005975(74)] Obs[0.002672(41)] Obs[0.001741(29)] Obs[0.000990(18)]]\n",
      " [Obs[0.002729(41)] Obs[0.001433(25)] Obs[0.000990(18)] Obs[0.000596(12)]]]\n"
     ]
    }
   ],
   "source": [
    "print(matrix_V1V1.content[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b965e080",
   "metadata": {},
   "source": [
    "Let us try to get the effective mass. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bd58d1f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Something is wrong\n"
     ]
    }
   ],
   "source": [
    "try:\n",
    "    matrix_V1V1.m_eff() #This does not work! \n",
    "except:\n",
    "    print(\"Something is wrong\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "044503dc",
   "metadata": {},
   "source": [
    "Many methods we could use for regular correlators do not work with matrix-correlators. \n",
    "In order to get the effective mass, we need to convert to a regular correlator first. \n",
    "\n",
    "One way to do it, is to pick an element out of the matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "90e54c2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "single_smearing = matrix_V1V1.item(0,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c25a23c",
   "metadata": {},
   "source": [
    "**Corr.index(i,j)** picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice. \n",
    "But there is a more usefull way to retrieve a single value per timeslice. \n",
    "We might want a linear combination of different sources and sinks. \n",
    "We can formalize this as\n",
    "\n",
    "$$C_{\\textrm{projected}}(t)=v_1^T \\underline{C}(t) v_2$$\n",
    "\n",
    "If we choose the vectors to be $v_1=v_2=(0,1,0,0)$, we should get the same correlator as in the cell above. \n",
    "\n",
    "Thinking about it this way is usefull in the Context of the generalized eigenvalue problem (GEVP), used to find the source-sink combination, which best describes a certain energy eigenstate.\n",
    "A good introduction is found in https://arxiv.org/abs/0902.1265."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1d0e8a06",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  -0.71920537  162.3903938  -492.21321233  714.06804297]\n"
     ]
    }
   ],
   "source": [
    "vec = matrix_V1V1.GEVP(t0=3,ts=6,state=0)\n",
    "assert len(vec) == matrix_V1V1.N\n",
    "print(vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d7aeb89",
   "metadata": {},
   "source": [
    "As we see, the eigenvector is of length **matrix_V1V1.N** and contains regular floats. \n",
    "We can use it to project the correlator. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fb0b7b57",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 19248 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "matrix_V1V1.projected(vec).m_eff().show(comp=single_smearing.m_eff())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f26cdaf2",
   "metadata": {},
   "source": [
    "There is a lot going on in this line of code. \n",
    "We start with our matrix correlator and we project it, using the vector we got from the GEVP routine. \n",
    "\n",
    "This gives us a new correlator with one Obs per timeslice. We then calculate its effective mass and plot it. \n",
    "We tell the **Corr.show** method to show another correlator as a comparison. \n",
    "\n",
    "We can see, that the projected correlator (*blue*) converges to a mass plateau much quicker than the single smearing level."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06c366ed",
   "metadata": {},
   "source": [
    "## Example Analysis\n",
    "\n",
    "We can use what we learned so far to perform an actually usefull analysis. \n",
    "The correlator **matrix_V1V1** we looked at corresponds to vector-charmonium. \n",
    "\n",
    "We might be interested in the mass of the $J/\\Psi$ state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "def98bc0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been written using pyerrors 2.0.0+dev.\n",
      "Format version 0.1\n",
      "Written by jan on 2022-01-27 10:49:51 +0100 on host endwings, Linux-5.13.0-27-generic-x86_64-with-glibc2.10\n",
      "Data has been written using pyerrors 2.0.0+dev.\n",
      "Format version 0.1\n",
      "Written by jan on 2022-01-27 10:49:51 +0100 on host endwings, Linux-5.13.0-27-generic-x86_64-with-glibc2.10\n",
      "Fit with 1 parameter\n",
      "Method: Levenberg-Marquardt\n",
      "`ftol` termination condition is satisfied.\n",
      "chisquare/d.o.f.: 0.692947327942094\n",
      "--- The mass was calculated to be 3079.8(7.9) MeV ---\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 18360 characters]\n",
      "text/plain": [
       "<Figure size 640x395.55 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# We do not just have V1V1, but also the two other spacial directions. We can average over them for better statistics. \n",
    "matrix_V2V2= pe.input.json.load_json(\"./data/matrix_correlator_V2V2\")\n",
    "matrix_V3V3= pe.input.json.load_json(\"./data/matrix_correlator_V3V3\")\n",
    "matrix_VnVn=(matrix_V1V1 + matrix_V2V2 + matrix_V3V3) / 3. \n",
    "\n",
    "#We then solve the GEVP to get eigenvectors corresponding to the ground state. \n",
    "\n",
    "vec_ground = matrix_VnVn.GEVP(t0=3, ts=6, state=0)\n",
    "\n",
    "#Now we project the matrix-correlators to get new correlators belonging to the ground state.\n",
    "\n",
    "corr_ground = matrix_VnVn.projected(vec_ground)\n",
    "\n",
    "# We get the effective mass using the periodic cosh method. \n",
    "m_eff_Jpsi = corr_ground.m_eff(variant=\"cosh\")\n",
    "\n",
    "# We can now pick a plateau range and get a single value for the mass. \n",
    "m_eff_Jpsi.gamma_method()\n",
    "m_Jpsi = m_eff_Jpsi.plateau([5, 24])\n",
    "\n",
    "# We can now visually compare our plateau value to the data\n",
    "\n",
    "m_eff_Jpsi.show([10, 25], plateau=m_Jpsi)\n",
    "\n",
    "# Since the lattice spacing is known, we can multiply with hbar*c/a to see it in physical units \n",
    "\n",
    "m_Jpsi = m_Jpsi * 197 / 0.0653\n",
    "\n",
    "#As a last step we call the gamma method to get the error \n",
    "\n",
    "m_Jpsi.gamma_method()\n",
    "\n",
    "print(\"--- The mass was calculated to be\" , m_Jpsi, \"MeV ---\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c55c6066",
   "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": 5
}
