{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6a4ccad5-2d9d-44ae-802f-455a9cc17c05",
   "metadata": {},
   "source": [
    "# Evolving a PDF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e0dcd0f",
   "metadata": {},
   "source": [
    "## Method 1: Using apply_pdf\n",
    "\n",
    "In this first part we will compute the eko and subsequently apply the initial PDF \"manually\" calling a dedicated function. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23fcf9a1-f18a-4178-8b67-062df92abbdc",
   "metadata": {},
   "source": [
    "Next, we need to load the PDF that we want to evolve. EKO uses the same interface as lhapdf to query for the actual values of PDFs. However, for the scope of this tutorial we want to avoid the complication of dealing with an external dependency. Therefore we will use the toy PDFs as they were established by the Les Houches benchmark setting. They are provided in the `banana-hep` package available from PyPI, so first run in your shell:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f24a7ec-1829-41ef-a20b-343b50a9393e",
   "metadata": {},
   "source": [
    "`$ pip install banana-hep`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "134094a0-8df4-4e5b-b90a-951483fbe8af",
   "metadata": {},
   "source": [
    "and then in your python interpreter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8bd8824b-de64-491f-b279-3ffaed97a6d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "import eko\n",
    "from banana import toy\n",
    "pdf = toy.mkPDF(\"\", 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5de940ea-1dbe-455c-91ae-274b26de62b1",
   "metadata": {},
   "source": [
    "Now, we have all ingredients at hand to evolve the PDF set with the operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e110fe95-7f57-4b09-a224-ab66da5d7d64",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ekobox.apply import apply_pdf\n",
    "with eko.EKO.read(\"./myeko.tar\") as evolution_operator:\n",
    "    evolved_pdfs = apply_pdf(evolution_operator, pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c2e752c-4a7c-4e2a-bd97-94f54907184a",
   "metadata": {},
   "source": [
    "The returned object `evolved_pdfs` is a dictionary which maps the requested final scales onto further dictionaries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "74ddb375-82a6-449f-984d-171ad315ef3e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys([10000.0])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evolved_pdfs.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ddf91bb-bc2c-4833-b624-32c5ac2d717a",
   "metadata": {},
   "source": [
    "Each final scale contains a dictionary with a key `pdfs`, where all PDF values are hold, and a key `errors` , where the integration errors are hold.\n",
    "Finally, inside each of those there is a mapping of Monte Carlo particle identifiers onto a the PDF value at the requested interpolation points.\n",
    "\n",
    "E.g. to access the gluon PDF at $Q^2 = 10000\\,\\text{GeV}^2$ you can run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6cab9fd1-f4c0-4fc4-aad8-6b6a6b91340c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.58966998e+04,  4.74768306e+02,  3.57939013e+01, -1.03946292e+01,\n",
       "        0.00000000e+00])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evolved_pdfs[10000.0][\"pdfs\"][21]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ead0923f-d7b7-4ce2-afba-3300c72c467a",
   "metadata": {},
   "source": [
    "Note that we return the actual PDF and not the momentum fraction times the PDF."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e925d2c9",
   "metadata": {},
   "source": [
    "## Method 2: Using evolve_pdfs\n",
    "\n",
    "In this second part we illustrate how to get (and install) directly a LHAPDF set evolved with eko. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6799ff5f",
   "metadata": {},
   "source": [
    "First, we define our initial PDF. Here, we will use the same toy PDF as in the previous example, but any LHAPDF-like object will do the job!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ca2ed7a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from banana import toy\n",
    "pdf = toy.mkPDF(\"\", 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbbae6d6",
   "metadata": {},
   "source": [
    "Now, we set the theory inputs: in this example we will evolve our toy PDF at LO and create a new LHAPDF object with\n",
    "a size two `mu2grid`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4a68494f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ekobox.cards import example\n",
    "\n",
    "th_card = example.theory()\n",
    "op_card = example.operator()\n",
    "# here we replace the grid with a very minimal one, to speed up the example\n",
    "op_card.xgrid = eko.interpolation.XGrid([1e-3, 1e-2, 1e-1, 5e-1, 1.])\n",
    "op_card.mugrid = [(10.,5), (100.,5)]\n",
    "# set QCD LO evolution\n",
    "th_card.orders = (1,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f1323db",
   "metadata": {},
   "source": [
    "Finally, we are ready to run eko and install the new PDF set.\n",
    "Note, that if the evolved PDF already exist the code will overwrite it.\n",
    "\n",
    "Additionally, you can set `path` to load a precomputed EKO, while setting `store_path` you can save the produced EKO and reuse it later.\n",
    "You can also iterate on the given PDF objects (e.g. replicas)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8c18c327",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Overwriting old PDF installation\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "install_pdf Evolved_PDF\n"
     ]
    }
   ],
   "source": [
    "from ekobox.evol_pdf import evolve_pdfs\n",
    "path = pathlib.Path(\"./myeko2.tar\")\n",
    "path.unlink(missing_ok=True)\n",
    "evolve_pdfs(\n",
    "    [pdf],\n",
    "    th_card,\n",
    "    op_card,\n",
    "    install=True,\n",
    "    name=\"Evolved_PDF\",\n",
    "    store_path=path\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "888e659d",
   "metadata": {},
   "source": [
    "Now, you can access the evolved PDF as all the other PDF sets (note that this requires the Python interface of lhapdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cf532a93-7a56-4e45-abd0-50e8652e9a59",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LHAPDF 6.4.0 loading /home/felix/local/share/LHAPDF/Evolved_PDF/Evolved_PDF_0000.dat\n",
      "Evolved_PDF PDF set, member #0, version 1\n"
     ]
    }
   ],
   "source": [
    "import lhapdf\n",
    "evolved_pdf = lhapdf.mkPDF(\"Evolved_PDF\", 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4622dc19",
   "metadata": {},
   "source": [
    "To obtain the value of the gluon PDF at a given scale you can simply do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9a85c78d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "has gluon? True\n",
      "xg(x=0.01, Q2=89.1) = 4.399576390180355\n"
     ]
    }
   ],
   "source": [
    "pid = 21 # gluon pid\n",
    "Q2 = 89.10 #  Q^2 in Gev^2\n",
    "x = 0.01 # momentum fraction \n",
    "\n",
    "# check that the particle is present\n",
    "print(\"has gluon?\",evolved_pdf.hasFlavor(pid))\n",
    "# now do the lookup\n",
    "xg = evolved_pdf.xfxQ2(pid, x, Q2)\n",
    "print(f\"xg(x={x}, Q2={Q2}) = {xg}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58b26d42",
   "metadata": {},
   "source": [
    "## A more Realistic Example: Benchmark to CT14llo\n",
    "\n",
    "In this part of the tutorial we do an eko benchmark showing how PDFs evolved with eko can reproduce the values from the original LHAPDF grids."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5afc630",
   "metadata": {},
   "source": [
    "First, we need to set up the theory and operator runcards to match the settings used to produce the chosen PDF, here we will use `CT14llo`.\n",
    "\n",
    "We have to use LO evolution and we choose to dump our PDF into grids with 5 values of `Q2` and 60 points in x-space logarithmically spaced between 1e-7 and 0.1 and linearly spaced from 0.1 to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "00f4f8a4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LHAPDF 6.4.0 loading /home/felix/local/share/LHAPDF/CT14llo/CT14llo_0000.dat\n",
      "CT14llo PDF set, member #0, version 1; LHAPDF ID = 13205\n"
     ]
    }
   ],
   "source": [
    "from math import nan\n",
    "import numpy as np\n",
    "import lhapdf\n",
    "from ekobox.cards import example\n",
    "from eko.interpolation import make_grid\n",
    "from eko.quantities.heavy_quarks import QuarkMassRef,HeavyQuarks\n",
    "\n",
    "# get the PDF object\n",
    "ct14llo = lhapdf.mkPDF(\"CT14llo\")\n",
    "\n",
    "# setup the operator card\n",
    "op_card = example.operator()\n",
    "op_card.xgrid = eko.interpolation.XGrid(make_grid(30, 30)) # x grid\n",
    "op_card.mugrid = [(float(q),5) for q in np.geomspace(5., 100, 5)] # Q2 grid\n",
    "op_card.mu0 = 1.295000 # starting point for the evolution \n",
    "\n",
    "# setup the theory card - this can be mostly inferred from the PDF's .info file\n",
    "\n",
    "th_card = example.theory()\n",
    "th_card.orders = (1,0) # QCD LO\n",
    "th_card.heavy.masses = HeavyQuarks([QuarkMassRef([1.3,nan]), QuarkMassRef([4.75,nan]), QuarkMassRef([172.,nan])]) # quark mass\n",
    "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n",
    "th_card.couplings.scale = 91.1876 # the reference scale at which alpha_s is provided\n",
    "th_card.couplings.num_flavs_ref = 5 # the number of flavors active at the alpha_s reference scale\n",
    "th_card.couplings.max_num_flavs = 5  # the maximum number of flavors active in the alpha_s evolution\n",
    "th_card.couplings.num_flavs_init = 3 # the number of flavors active at the reference scale\n",
    "th_card.num_flavs_max_pdf = 5  # the maximum number of flavors active in the pdf evolution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0a51d8b",
   "metadata": {},
   "source": [
    "Next, we run the evolution using method 2 and save the new PDF. Due to the extended x grid and Q2 grid this might take a minute so please be patient ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "05b81dca",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Overwriting old PDF installation\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "install_pdf my_ct14llo\n"
     ]
    }
   ],
   "source": [
    "from ekobox.evol_pdf import evolve_pdfs\n",
    "path = pathlib.Path(\"./myeko_ct14llo.tar\")\n",
    "path.unlink(missing_ok=True)\n",
    "evolve_pdfs(\n",
    "    [ct14llo],\n",
    "    th_card,\n",
    "    op_card,\n",
    "    install=True,\n",
    "    name=\"my_ct14llo\",\n",
    "    store_path=path\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df0e805f",
   "metadata": {},
   "source": [
    "Now, we can compare the values given by the original PDF set and the one evolved with eko, both at different `x` and `Q2` scales, for a chosen parton,\n",
    "here we look at the gluon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2b02aab1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LHAPDF 6.4.0 loading /home/felix/local/share/LHAPDF/my_ct14llo/my_ct14llo_0000.dat\n",
      "           x            Q2       ct14llo    my_ct14llo  relative_diff\n",
      "0   0.000010     25.000000  7.635785e+01  7.630461e+01       0.000697\n",
      "1   0.000173     25.000000  3.194273e+01  3.192092e+01       0.000683\n",
      "2   0.003000     25.000000  1.081843e+01  1.081086e+01       0.000701\n",
      "3   0.051962     25.000000  1.958956e+00  1.958632e+00       0.000166\n",
      "4   0.900000     25.000000  1.922415e-05  1.955026e-05      -0.016963\n",
      "5   0.000010    111.803399  1.333957e+02  1.332985e+02       0.000729\n",
      "6   0.000173    111.803399  4.777286e+01  4.773664e+01       0.000758\n",
      "7   0.003000    111.803399  1.341028e+01  1.339967e+01       0.000791\n",
      "8   0.051962    111.803399  1.978216e+00  1.978130e+00       0.000044\n",
      "9   0.900000    111.803399  6.644805e-06  6.753652e-06      -0.016381\n",
      "10  0.000010    500.000000  1.967032e+02  1.965456e+02       0.000801\n",
      "11  0.000173    500.000000  6.291393e+01  6.286095e+01       0.000842\n",
      "12  0.003000    500.000000  1.542347e+01  1.540996e+01       0.000876\n",
      "13  0.051962    500.000000  1.947465e+00  1.947391e+00       0.000038\n",
      "14  0.900000    500.000000  2.929060e-06  2.977306e-06      -0.016471\n",
      "15  0.000010   2236.067977  2.633266e+02  2.631109e+02       0.000819\n",
      "16  0.000173   2236.067977  7.708540e+01  7.701938e+01       0.000856\n",
      "17  0.003000   2236.067977  1.700410e+01  1.698928e+01       0.000872\n",
      "18  0.051962   2236.067977  1.893923e+00  1.893971e+00      -0.000025\n",
      "19  0.900000   2236.067977  1.544450e-06  1.570997e-06      -0.017189\n",
      "20  0.000010  10000.000000  3.314097e+02  3.311351e+02       0.000829\n",
      "21  0.000173  10000.000000  9.023010e+01  9.015279e+01       0.000857\n",
      "22  0.003000  10000.000000  1.825934e+01  1.824402e+01       0.000839\n",
      "23  0.051962  10000.000000  1.830992e+00  1.831183e+00      -0.000104\n",
      "24  0.900000  10000.000000  9.288458e-07  9.447927e-07      -0.017169\n",
      "my_ct14llo PDF set, member #0, version 1\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# load evolved pdf\n",
    "my_ct14llo = lhapdf.mkPDF(\"my_ct14llo\", 0)\n",
    "\n",
    "pid = 21 # gluon pid\n",
    "\n",
    "# collect data\n",
    "log = {\"x\": [], \"Q2\" : [],  \"ct14llo\": [], \"my_ct14llo\": [], \"relative_diff\": []} \n",
    "for q in np.geomspace(5., 100, 5):\n",
    "    q2 = q**2.\n",
    "    for x in np.geomspace(1e-5, 0.9, 5):\n",
    "        value = ct14llo.xfxQ2(pid, x, q2)\n",
    "        my_value =  my_ct14llo.xfxQ2(pid, x, q2)\n",
    "        log[\"x\"].append(x)\n",
    "        log[\"Q2\"].append(q2)\n",
    "        log[\"ct14llo\"].append(value)\n",
    "        log[\"my_ct14llo\"].append(my_value)\n",
    "        log[\"relative_diff\"].append((value - my_value) / value)\n",
    "\n",
    "print(pd.DataFrame(log))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6805c26d",
   "metadata": {},
   "source": [
    "As you can see EKO is able to reproduce the numbers from the original LHAPDF grid mostly below the permille level.\n",
    "\n",
    "The accuracy is mainly limited by the number of points in the `x` and `Q2` grids that can be finer to achieve higher precision.\n",
    "\n",
    "You can also notice that at large-x the gluon pdf vanishes so the worst accuracy of our benchmark is not worrying at all. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "eko-KkPVjVhh-py3.10",
   "language": "python",
   "name": "eko-kkpvjvhh-py3.10"
  },
  "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.10.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "7e472fcaae86e99be842e6a461f25738bfe926fddef76aeba76e9314d132dd8f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}