diff --git a/examples/python/dyaa.py b/examples/python/dyaa.py deleted file mode 100755 index a9620b57..00000000 --- a/examples/python/dyaa.py +++ /dev/null @@ -1,184 +0,0 @@ -#!/usr/bin/env python -import argparse - -import numpy as np -import pineappl - -def int_photo(s, t, u): - """ - Photon-photon matrix element. - - Parameters - ---------- - s : float - Mandelstam s - t : float - Mandelstam t - u : float - Mandelstam u - - Returns - ------- - float : - matrix element - """ - alpha0 = 1.0 / 137.03599911 - return alpha0 * alpha0 / 2.0 / s * (t / u + u / t) - - -def hadronic_pspgen(mmin, mmax): - r""" - Hadronic phase space generator. - - Parameters - ---------- - mmin : float - minimal energy :math:`\sqrt{s_{min}}` - mmax : float - maximal energy :math:`\sqrt{s_{max}}` - - Returns - ------- - s : float - Mandelstam s - t : float - Mandelstam t - u : float - Mandelstam u - x1 : float - first momentum fraction - x2 : float - second momentum fraction - jacobian : float - jacobian from the uniform generation - """ - smin = mmin * mmin - smax = mmax * mmax - - r1 = np.random.uniform() - r2 = np.random.uniform() - r3 = np.random.uniform() - - tau0 = smin / smax - tau = pow(tau0, r1) - y = pow(tau, 1.0 - r2) - x1 = y - x2 = tau / y - s = tau * smax - - jacobian = tau * np.log(tau0) * np.log(tau0) * r1 - - # theta integration (in the CMS) - cos_theta = 2.0 * r3 - 1.0 - jacobian *= 2.0 - - t = -0.5 * s * (1.0 - cos_theta) - u = -0.5 * s * (1.0 + cos_theta) - - # phi integration - jacobian *= 2.0 * np.math.acos(-1.0) - - return [s, t, u, x1, x2, jacobian] - - -def fill_grid(grid, calls): - """ - Fill grid with points. - - Parameters - ---------- - grid : pineappl.Grid - grid to fill - calls : int - number of events - """ - - # in GeV^2 pbarn - hbarc2 = 389379372.1 - - for _ in range(calls): - # compute phase space - s, t, u, x1, x2, jacobian = hadronic_pspgen(10.0, 7000.0) - # build observables - ptl = np.sqrt((t * u / s)) - mll = np.sqrt(s) - yll = 0.5 * np.log(x1 / x2) - ylp = np.abs(yll + np.math.acosh(0.5 * mll / ptl)) - ylm = np.abs(yll - np.math.acosh(0.5 * mll / ptl)) - - jacobian *= hbarc2 / calls - - # cuts for LO for the invariant-mass slice containing the - # Z-peak from CMSDY2D11 - if ( - ptl < 14.0 - or np.abs(yll) > 2.4 - or ylp > 2.4 - or ylm > 2.4 - or mll < 60.0 - or mll > 120.0 - ): - continue - - # build event - weight = jacobian * int_photo(s, u, t) - q2 = 90.0 * 90.0 - # put - grid.fill(x1, x2, q2, 0, np.abs(yll), 0, weight) - - -def main(calls, pdfname, filename): - """ - Generate the grid. - - Parameters - ---------- - calls : int - number of events - pdfname : str - if given, write the predictions - filename : str - if given, write the grid to this path - """ - # create a new luminosity function for the $\gamma\gamma$ initial state - lumi_entries = [pineappl.lumi.LumiEntry([(22, 22, 1.0)])] - # only LO $\alpha_\mathrm{s}^0 \alpha^2 \log^0(\xi_\mathrm{R}) \log^0(\xi_\mathrm{F})$ - orders = [pineappl.grid.Order(0, 2, 0, 0)] - bins = np.arange(0, 2.4, 0.1) - params = pineappl.subgrid.SubgridParams() - grid = pineappl.grid.Grid.create(lumi_entries, orders, bins, params) - - # fill the grid with phase-space points - print(f"Generating {calls} events, please wait...") - fill_grid(grid, calls) - - # load pdf for testing - if pdfname: - import lhapdf # pylint: disable=import-outside-toplevel - - pdf = lhapdf.mkPDF(pdfname, 0) - pdg_id = int(pdf.set().get_entry('Particle')) - # perform convolution - dxsec = grid.convolve_with_one(pdg_id, pdf.xfxQ2, pdf.alphasQ2) - for i in range(len(dxsec)): - print(f"{bins[i]:.1f} {bins[i + 1]:.1f} {dxsec[i]:.3e}") - - # write the grid to disk - if filename: - print(f"Writing PineAPPL grid to disk: {filename}") - grid.write(filename) - - -if __name__ == "__main__": - # expose the arguments as script - ap = argparse.ArgumentParser() - ap.add_argument("--calls", default=100000, type=int, help="Number of events") - ap.add_argument( - "--pdf", - default="NNPDF31_nlo_as_0118_luxqed", - type=str, - help="plot with PDF set", - ) - ap.add_argument("--filename", type=str, help="Output path") - args = ap.parse_args() - main(args.calls, args.pdf, args.filename) diff --git a/examples/python/positivity.py b/examples/python/positivity.py index 426a3600..26e05819 100755 --- a/examples/python/positivity.py +++ b/examples/python/positivity.py @@ -11,14 +11,15 @@ def main(filename, Q2): pid = 4 # init pineappl objects - lumi_entries = [pineappl.lumi.LumiEntry([(pid, lepton_pid, 1.0)])] + lumi_entries = [pineappl.boc.Channel([(pid, lepton_pid, 1.0)])] orders = [pineappl.grid.Order(0, 0, 0, 0)] bins = len(xgrid) - bin_limits = list(map(float, range(0, bins + 1))) + # NOTE: `bin_limits` have to be `np.ndarray` + bin_limits = np.array([float(i) for i in range(0, bins + 1)]) # subgrid params - default is just sufficient params = pineappl.subgrid.SubgridParams() # inti grid - grid = pineappl.grid.Grid.create(lumi_entries, orders, bin_limits, params) + grid = pineappl.grid.Grid(lumi_entries, orders, bin_limits, params) limits = [] # add each point as a bin for bin_, x in enumerate(xgrid): @@ -31,13 +32,15 @@ def main(filename, Q2): # create and set subgrid = pineappl.import_only_subgrid.ImportOnlySubgridV1( array[np.newaxis, :, np.newaxis], - [Q2], - xgrid, - [1.0], + np.array([Q2]), # `q2_grid` has to be `np.ndarrary` + np.array(xgrid), # `x_grid` has to be `np.ndarrary` + np.array([1.0]), # `x_grid` has to be `np.ndarrary` ) - grid.set_subgrid(0, bin_, 0, subgrid) + grid.set_subgrid(0, bin_, 0, subgrid.into()) # set the correct observables - normalizations = [1.0] * bins + normalizations = np.array( + [1.0] * bins + ) # `normalizations` has to be `np.ndarray` remapper = pineappl.bin.BinRemapper(normalizations, limits) grid.set_remapper(remapper) @@ -55,4 +58,4 @@ def main(filename, Q2): if __name__ == "__main__": - main("charm-positivity.pineappl", 1.5 ** 2) + main("charm-positivity.pineappl", 1.5**2) diff --git a/pineappl_py/.gitignore b/pineappl_py/.gitignore index 83eb351c..be08e5e8 100644 --- a/pineappl_py/.gitignore +++ b/pineappl_py/.gitignore @@ -1,4 +1,5 @@ - +.ipynb_checkpoints +docs/source/*.pineappl.lz4 ### Python ### # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/pineappl_py/docs/source/advanced.ipynb b/pineappl_py/docs/source/advanced.ipynb new file mode 100644 index 00000000..b499b128 --- /dev/null +++ b/pineappl_py/docs/source/advanced.ipynb @@ -0,0 +1,447 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8f2d1864-a5c6-44e0-b226-b615814cb29f", + "metadata": {}, + "source": [ + "# A (slightly) Advanced Tutorial\n", + "\n", + "In the following tutorial, we will calculate **leading-order (LO)**\n", + "cross section of the production of a same-flavour opposite-sign (SFOS) \n", + "lepton-pair (also known as Drell-Yan lepton-pair production) at the \n", + "LHC: $\\mathrm{p}\\mathrm{p} \\to \\ell\\bar{\\ell}$ @ 7 TeV. In particular,\n", + "we are going to look at the differential distribution in the rapidity\n", + "of the lepton pair, in the setup given by .\n", + "\n", + "That is, we'll write a Monte Carlo integrator that calculates a part of \n", + "this process and produces an interpolation grid. The following steps will\n", + "provide a tangible illustration on how to fill a PineAPPL grid with the\n", + "various ingredients. These steps involve the computation of the matrix\n", + "element and the generation of the phase space." + ] + }, + { + "cell_type": "markdown", + "id": "2feab5c9-ddda-40a7-be1d-6111bdc566e5", + "metadata": {}, + "source": [ + "## Computing an observable\n", + "\n", + "A physical observable that involves two hadrons in the initial states is computed as:\n", + "$$ \\left\\langle \\frac{d\\sigma^{hh \\to F}}{\\mathrm{d} \\mathcal{O}} \\right\\rangle = \\sum_{a,b} \\int \\mathrm{d} x_1 \\mathrm{d} x_2 f_a^h (x_1) f_b^h (x_2) \\frac{d\\sigma_{ab \\to F} (x_1, x_2)}{\\mathrm{d} \\mathcal{O}} $$\n", + "where the partonic cross section is a perturbative series in the two couplings (strong coupling: $\\alpha_s (M_\\mathrm{Z}^2) = 0.118$ and electromagnetic coupling $\\alpha(0) \\approx 1/137 \\approx 0.0073$):\n", + "$$ \\frac{d\\sigma_{ab} (x_1, x_2)}{\\mathrm{d} \\mathcal{O}} = \\sum_{n,m} \\alpha_\\mathrm{s}^n \\alpha^m \\frac{d\\sigma_{ab}^{n,m} (x_1, x_2)}{\\mathrm{d} \\mathcal{O}} $$\n", + "Since $\\alpha_s \\gg \\alpha$ we usually look only at the lowest order $m$ and calculate corrections in $n$: this is what we refer as QCD corrections. However, this isn't always reliable, sometimes electroweak (EW) corrections are needed.\n", + "\n", + "Inserting the perturbative expansion into the main formula:\n", + "$$ \\left\\langle \\frac{d\\sigma^{hh \\to F}}{\\mathrm{d} \\mathcal{O}} \\right\\rangle = \\sum_{a,b} \\sum_{n,m} \\int \\mathrm{d} x_1 \\mathrm{d} x_2 f_a^h (x_1) f_b^h (x_2) \\alpha_\\mathrm{s}^n \\alpha^m \\frac{d\\sigma_{ab \\to F}^{n,m} (x_1, x_2)}{\\mathrm{d} \\mathcal{O}} $$\n", + "\n", + "We call $d\\sigma_{ab \\to F}^{n,m} (x_1, x_2) / \\mathrm{d} \\mathcal{O}$ the **interpolation grid**. It has the advantage that one can very quickly (less than a second) perform the integrals above with any PDF set, which is very important for PDF extraction." + ] + }, + { + "cell_type": "markdown", + "id": "54d8e0b0-d664-4f49-9bf3-a167c82dd5a1", + "metadata": {}, + "source": [ + "## Compute matrix elements\n", + "\n", + "The first step in computing theory predictions is the computation of the (partonic) matrix elements (amplitudes). The next step is to sum all the amplitudes and take the modulus squared. It is common practice to also account for the flux factor and the spin and color sums together with their eventual average. Recall to average on the input and to sum on the output. In our example we find:\n", + "$$ \\frac {1}{2 s} |\\mathcal M_t + \\mathcal M_u |^2 = \\frac{\\alpha^2}{2s} \\left(\\frac t u + \\frac u t\\right) $$" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "6bc0bf80-787e-427c-aec4-f835309eaf76", + "metadata": {}, + "outputs": [], + "source": [ + "def photon_photon_matrix_element(s: float, t: float, u: float) -> float:\n", + " alpha0 = 1.0 / 137.03599911\n", + " return alpha0 * alpha0 / 2.0 / s * (t / u + u / t)" + ] + }, + { + "cell_type": "markdown", + "id": "cb8f5606-83ca-45e0-bc1a-87600b0cc1de", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Determine phase space decomposition\n", + "\n", + "Given the initial states with momenta $k_1$ and $k_2$ we need to integrate the squared matrix elements over all possible momenta, that is all momenta which fulfill momentum conservation and which are on-shell: $ p_i^2 = m_i^2 $. In general the Lorentz invariant phase-space (LIPS) for $n$ particles is\n", + "\n", + "$$ \\int \\mathrm{d} \\mathrm{LIPS} = \\int \\left( \\prod_{i=1}^n \\mathrm{d}^4 p_i \\right) \\, \\delta^{(4)} \\left( k_1 + k_2 - \\sum_{i=1}^n p_i \\right) \\prod_{i = 1}^n \\delta \\left( p_i^2 - m_i^2 \\right) $$\n", + "\n", + "and has $4n$ integration dimensions, reduced to $3n - 4$ through the momentum conservation ($-4$) and on-shell conditions ($-n$).\n", + "\n", + "In our example we have two massless final state particles ($n = 2$ and $m_1 = m_2 = 0$), so effectively we integrate over $3n - 4 = 2$ dimensions. We choose to integrate over these two variables:\n", + "1. $\\cos \\theta$, where $\\theta$ measures the angle of one of the leptons w.r.t. the beam axis and\n", + "2. the angle $\\phi$, which is another angle transverse to the beam axis.\n", + "\n", + "Matrix elements do not depend on the angle $\\phi$, since the collision is symmetric around the beam axis." + ] + }, + { + "cell_type": "markdown", + "id": "f84f85c0-be69-4007-9871-9ce4797f4fe3", + "metadata": {}, + "source": [ + "## Compute phase space integrals\n", + "\n", + "Our master formula,\n", + "\n", + "$$ \\sigma^{\\mathrm{p}\\mathrm{p} \\to \\ell\\bar{\\ell}} = \\sum_{a,b} \\int \\mathrm{d} x_1 \\mathrm{d} x_2 f_a^\\mathrm{p} (x_1) f_b^\\mathrm{p} (x_2) \\sigma_{ab \\to \\ell\\bar{\\ell}} (x_1, x_2) $$\n", + "\n", + "requires us to integrate over all possible momentum fractions $x_1$ and $x_2$ of the two PDFs. We do this by rewriting the integral into $\\tau$, relative centre-of-mass energy squared and $y$, the rapidity relating the hadronic and partonic centre-of-mass frames:\n", + "\n", + "$$ \\int_0^1 \\mathrm{d} x_1 \\int_0^1 \\mathrm{d} x_2 = \\int \\mathrm{d} \\tau \\int \\mathrm{d} y $$\n", + "\n", + "The exact form of this transformation isn't really important, but it is chosen such that the jacobian contains the inverse flux factor, cancelling the flux factor multiplied to the squared matrix elements above.\n", + "\n", + "We approximate the integrals numerically by using a Monte Carlo integration, which computes the average of the integrand evaluated using uniformly chosen random numbers $r_1, r_2, r_3$:\n", + "\n", + "$$ \\int_0^1 \\mathrm{d} r_1 \\int_0^1 \\mathrm{d} r_2 \\int_0^1 \\mathrm{d} r_3 f(r_1, r_2, r_3) \\approx \\frac{1}{N} \\sum_{i=1}^N f(r_1^i, r_2^i, r_3^i) $$\n", + "\n", + "Translated to Python code this reads:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "357d81ff-5c6b-4514-82f7-e5575b80215a", + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "import numpy as np\n", + "from typing import Tuple\n", + "\n", + "np.random.seed(1234567890)\n", + "\n", + "def hadronic_ps_gen(\n", + " mmin: float, mmax: float\n", + ") -> Tuple[float, float, float, float, float, float]:\n", + " r\"\"\"Hadronic phase space generator.\n", + "\n", + " Parameters\n", + " ----------\n", + " mmin :\n", + " minimal partonic centre-of-mass energy :math:`\\sqrt{s_{min}}`\n", + " mmax :\n", + " maximal partonic centre-of-mass energy :math:`\\sqrt{s_{max}}`\n", + "\n", + " Returns\n", + " -------\n", + " s :\n", + " Mandelstam s\n", + " t :\n", + " Mandelstam t\n", + " u :\n", + " Mandelstam u\n", + " x1 :\n", + " first momentum fraction\n", + " x2 :\n", + " second momentum fraction\n", + " jacobian :\n", + " jacobian from the uniform generation\n", + "\n", + " \"\"\"\n", + " smin = mmin * mmin\n", + " smax = mmax * mmax\n", + "\n", + " r1 = np.random.uniform()\n", + " r2 = np.random.uniform()\n", + " r3 = np.random.uniform()\n", + " \n", + " # generate partonic x1 and x2\n", + " tau0 = smin / smax\n", + " tau = pow(tau0, r1)\n", + " y = pow(tau, 1.0 - r2)\n", + " x1 = y\n", + " x2 = tau / y\n", + " s = tau * smax\n", + "\n", + " jacobian = tau * np.log(tau0) * np.log(tau0) * r1\n", + "\n", + " # theta integration (in the CMS)\n", + " cos_theta = 2.0 * r3 - 1.0\n", + " jacobian *= 2.0\n", + "\n", + " # reconstruct invariants (in the CMS)\n", + " t = -0.5 * s * (1.0 - cos_theta)\n", + " u = -0.5 * s * (1.0 + cos_theta)\n", + "\n", + " # phi integration\n", + " jacobian *= 2.0 * math.acos(-1.0)\n", + "\n", + " return [s, t, u, x1, x2, jacobian]" + ] + }, + { + "cell_type": "markdown", + "id": "65f003f3", + "metadata": {}, + "source": [ + "Now we can test the integration by generating a phase-space point between $s_\\text{min} = (10~\\text{GeV})^2$ and $s_\\text{max} = (7000~\\text{GeV})^2$ (our hadronic centre-of-mass energy):" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "5f1be2c6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Values of the Mandelstam variables:\n", + "s = 1.476157e+04\n", + "t = -1.643205e+03\n", + "u = -1.311836e+04\n", + "\n", + "Values of the partonic variables x1 and x2:\n", + "x1 = 3.648218e-02\n", + "x2 = 8.257633e-03\n", + "\n", + "Check the sum of the Mandelstam variables:\n", + "s+t+u = 0.000000e+00\n" + ] + } + ], + "source": [ + "[s, t, u, x1, x2, jacobian] = hadronic_ps_gen(10.0, 7000.0)\n", + "\n", + "# print(\n", + "# \"s = {:.6e}\\nt = {:.6e}\\nu = {:.6e}\\n\\nx1 = {:.6e}\\nx2 = {:.6e}\\n\\ns + t + u = {:.6e}\"\n", + "# .format(s, t, u, x1, x2, s + t + u)\n", + "# )\n", + "print(\"Values of the Mandelstam variables:\")\n", + "print(f\"s = {s:.6e}\\nt = {t:.6e}\\nu = {u:.6e}\")\n", + "\n", + "print(\"\\nValues of the partonic variables x1 and x2:\")\n", + "print(f\"x1 = {x1:.6e}\\nx2 = {x2:.6e}\")\n", + "\n", + "print(\"\\nCheck the sum of the Mandelstam variables:\")\n", + "print(f\"s+t+u = {s+t+u:.6e}\")" + ] + }, + { + "cell_type": "markdown", + "id": "8a135784-1ce2-49f2-8c21-b1447c8eb83a", + "metadata": {}, + "source": [ + "## Join phase space integration and matrix elements\n", + "\n", + "Finally, we have to\n", + "- put the integral together with the squared matrix elements,\n", + "- transform the phase-space variables into the well-known LAB quantities, and\n", + "- we want to simulate the setup from CMS DY 7 TeV, see: \n", + "\n", + "This means, we need to add phase-space cuts:\n", + " $$ p_\\mathrm{T}^\\ell > 14~\\text{GeV}, \\quad |y^\\ell| < 2.4, \\quad 60~\\text{GeV} < M_{\\ell\\bar{\\ell}} < 120~\\text{GeV}$$\n", + "and we want the differential cross section w.r.t. $|y_{\\ell\\bar{\\ell}}|$, with bin limits $0 < |y_{\\ell\\bar{\\ell}}| < 2.4$, in steps of $0.1$." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "5bb89769-55ad-43b8-8844-41538ae2e9a1", + "metadata": {}, + "outputs": [], + "source": [ + "import pineappl\n", + "\n", + "def fill_grid(grid: pineappl.grid.Grid, calls: int):\n", + " \"\"\"Fill grid with points.\"\"\"\n", + "\n", + " # in GeV^2 pbarn\n", + " hbarc2 = 389379372.1\n", + "\n", + " # perform Monte Carlo sum\n", + " for _ in range(calls):\n", + " # compute phase space\n", + " s, t, u, x1, x2, jacobian = hadronic_ps_gen(10.0, 7000.0)\n", + "\n", + " # build observables\n", + " ptl = np.sqrt((t * u / s))\n", + " mll = np.sqrt(s)\n", + " yll = 0.5 * np.log(x1 / x2)\n", + " ylp = np.abs(yll + math.acosh(0.5 * mll / ptl))\n", + " ylm = np.abs(yll - math.acosh(0.5 * mll / ptl))\n", + "\n", + " # apply conversion factor\n", + " jacobian *= hbarc2 / calls\n", + "\n", + " # cuts for LO for the invariant-mass slice containing the Z-peak from CMS (7 TeV): https://arxiv.org/abs/1310.7291\n", + " if (\n", + " ptl < 14.0\n", + " or np.abs(yll) > 2.4\n", + " or ylp > 2.4\n", + " or ylm > 2.4\n", + " or mll < 60.0\n", + " or mll > 120.0\n", + " ):\n", + " # continuing means this we don't call fill below that means this event counts as zero or it 'cut away'\n", + " continue\n", + "\n", + " # build event\n", + " weight = jacobian * photon_photon_matrix_element(s, u, t)\n", + " # set factorization and renormalization scale to (roughly) the Z-boson mass\n", + " q2 = 90.0 * 90.0\n", + " \n", + " # fill the interpolation grid\n", + " grid.fill(x1, x2, q2, 0, np.abs(yll), 0, weight)" + ] + }, + { + "cell_type": "markdown", + "id": "f65fcfc8-2fd9-45e2-9bff-df12e0759392", + "metadata": {}, + "source": [ + "We want our results stored in an interpolation grid, which is independent of PDFs and the strong coupling. To create a `Grid`, we need to give it a few bits of information. We have to tell it that\n", + "\n", + "- our initial state is photon-photon, or in PDG Monte Carlo IDs `(22, 22)`\n", + "- the perturbative order in $\\alpha^2$\n", + "- as per CMS's setup we bin the observable from $0$ to $2.4$ in steps of $0.1$." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "26191d31-ea97-4b57-880c-2a20713cff6b", + "metadata": {}, + "outputs": [], + "source": [ + "def generate_grid(calls: int) -> pineappl.grid.Grid:\n", + " \"\"\"Generate the grid.\"\"\"\n", + " # create a new luminosity function for the $\\gamma\\gamma$ initial state\n", + " lumi_entries = [pineappl.boc.Channel([(22, 22, 1.0)])]\n", + " # only LO $\\alpha_\\mathrm{s}^0 \\alpha^2 \\log^0(\\xi_\\mathrm{R}) \\log^0(\\xi_\\mathrm{F})$\n", + " orders = [pineappl.grid.Order(0, 2, 0, 0)]\n", + " bins = np.arange(0, 2.4, 0.1)\n", + " params = pineappl.subgrid.SubgridParams()\n", + " grid = pineappl.grid.Grid(lumi_entries, orders, bins, params)\n", + "\n", + " # fill the grid with phase-space points\n", + " print(f\"Generating {calls} events, please wait...\")\n", + " fill_grid(grid, calls)\n", + " print(\"Done.\")\n", + "\n", + " return grid" + ] + }, + { + "cell_type": "markdown", + "id": "6823abe6", + "metadata": {}, + "source": [ + "We have to play a bit with the Monte Carlo statistics, to produce smooth results. To generate the full theory predictions, we must also use our master formula and convolute the interpolation grid with the two photon PDFs. Finally, let's plot the result:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "94775264-b3e1-4d3f-b8a6-596ed52da98f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Generating 1000000 events, please wait...\n", + "Done.\n", + "LHAPDF 6.5.0 loading /Users/tanjona/miniconda3/envs/nnpdf/share/LHAPDF/NNPDF31_nnlo_as_0118_luxqed/NNPDF31_nnlo_as_0118_luxqed_0000.dat\n", + "NNPDF31_nnlo_as_0118_luxqed PDF set, member #0, version 2; LHAPDF ID = 325100\n" + ] + } + ], + "source": [ + "import lhapdf\n", + "\n", + "# generate interpolation grid: increase this number!\n", + "grid = generate_grid(1000000)\n", + "\n", + "# perform convolution with PDFs: this performs the x1 and x2 integrals\n", + "# of the partonic cross sections with the PDFs as given by our master \n", + "# formula\n", + "pdf = lhapdf.mkPDF(\"NNPDF31_nnlo_as_0118_luxqed\", 0)\n", + "bins = grid.convolve_with_one(2212, pdf.xfxQ2, pdf.alphasQ2)" + ] + }, + { + "cell_type": "markdown", + "id": "8d665de1-ac23-40df-b219-a30e3e6f8475", + "metadata": {}, + "source": [ + "**NOTE:** If you do not have `NNPDF31_nnlo_as_0118_luxqed` installed, you can do so with the following command:\n", + "\n", + "```\n", + " !lhapdf install NNPDF31_nnlo_as_0118_luxqed\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "fe14ee4b-8309-4d78-9682-8bd37b91992e", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAFvCAYAAABAYhLAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsVUlEQVR4nO3df3RU9Z3/8dfMQDIECAYjQSA1QiwISiKhCdBScBdlXStq9btgXaHZSHfVVDGrrmm3sKg1UiDiEQSrRlpdhfagrrt2URuI9QcWl8BBUCjBQqCSEMqPISFMwtz7/QMyNU2AzMydmZu5z8c5OWfuzed+8p7PSfI69zP3fq7LNE1TAAAg4bnjXQAAAIgNQh8AAIcg9AEAcAhCHwAAhyD0AQBwCEIfAACHIPQBAHCIHvEuIN4Mw9CXX36pvn37yuVyxbscAABCYpqmjh8/rkGDBsntPve5vOND/8svv1RmZma8ywAAICL79u3TkCFDztnG8aHft29fSdKePXuUlpYW52qcJRAIaPfu3Ro2bJg8Hk+8y3EMxj0+GPf4SfSx9/l8yszMDObZuTg+9Num9FNTU5WamhrnapwlEAioT58+Sk1NTcg/RLti3OODcY8fp4x9Vz6i5kI+AAAcgtAHAMAhCH0AAByC0AcAwCEIfQAAHILQBwDAIQh9AAAcgtAHAMAhCH0AAByC0AcAwCEcvwxvUEuT1JIUeT89UySe1gcAsCFC/wzPkpFSsgVhPSRfKnqH4AcA2A7T+1bbv1FqPhLvKgAA6IAz/TMCt7wopfUPv4NTJ6VXZ5x+bRjWFAUAgIUI/TbJqVKvC8I/vrXZslIAAIgGQt8JTFNqPWFdf1ysCADdEqGf6ExTqpgq7fu9dX1ysSIAdEtcyJfoWk9YG/gSFysCQDfFmb6T/L9fSt7U8I/nYkUA6NYIfSfxcrEiADgZ0/sAADgEoQ8AgEMQ+gAAOAShDwCAQxD6AAA4BKEPAIBDcMueXVm1dG6LhcvvAgC6NULfjqKxdC4AwPGY3rejaCydm/51qWcva/sEAHQrnOnbXaRL57bp2UvqkRx5PwCAbovQt7tIl84FAOAMpvcBAHAIQh8AAIcg9AEAcAg+00d4Wk9ILU2R9REInL49EQAQE4R+NEQaiN1hQZ2nroy4C4+kr6WPli5bH3k9AIDzIvSjwYJAtKUeXiljlFS/3bIuUw5tVaD5qNQ33bI+AQCdI/StEoVAtN2COi6XNG2p5PdFPi1/6qT06ozTrw0j8toAAOdF6FvFykBsY8cFdVwuydsv8n5amyPvAwAQEkLfSlYFIgAAUcAtewAAOARn+og/K27/65lyeqYFAHBWtgz9ZcuWaeHChaqrq1NOTo6efvpp5efnd9p25cqVKiwsbLcvOTlZJ0+ejEWpsIBnaU7knQzJl4reIfgB4BxsN72/evVqlZSUaN68eaqurlZOTo6mTp2qgwcPnvWY1NRUHThwIPi1d+/eGFaMsPTwyhwwyrr+9m+Umo9Y1x8AJCDbnemXl5dr9uzZwbP3FStW6K233lJFRYUefvjhTo9xuVwaOHBgLMtEpFwuGTc8rZraA8oekCKPO8wzdG79A4Aus1Xot7S0aNOmTSotLQ3uc7vdmjJlijZs2HDW4xobG3XJJZfIMAyNGTNGjz/+uEaN6vws0u/3y+/3B7d9Pp8kKWCYChgsCRtLAVMK9ExRILmfFG7ou5PkaesvYJxe2hfnFAgEZBiGAoxVTDHu8ZPoYx/K+7JV6B86dEiBQEAZGRnt9mdkZGjHjh2dHjN8+HBVVFRo9OjROnbsmBYtWqQJEyZo+/btGjJkSIf2ZWVlmj9/fof9Xxw6odSWY9a8EXSJYUqHj/tVI1/Yme86dVLDz7yu2bNHZvJhy+pLVIZh6PDhw6qpqZHbbbtP+BIW4x4/iT72jY2NXW5rq9APx/jx4zV+/Pjg9oQJE3T55Zfr2Wef1aOPPtqhfWlpqUpKSoLbPp9PmZmZGpqeorQLucc+lgKGqRqZyh6YGv70fmtS8GV2VpbU50JriktggUBANTU1ys7OlsfjOf8BsATjHj+JPvZtM9ZdYavQT09Pl8fjUX19fbv99fX1Xf7MvmfPnrrqqqtUU1PT6feTk5OVnNxxlTuP2xV+8CBsbpcrsrH/ynEej1tKwD/oaHC73fJ4PAn5D9DOGPf4SeSxD+U92WqeIykpSXl5eaqsrAzuMwxDlZWV7c7mzyUQCOjTTz/VxRdfHK0yAQDolmx1pi9JJSUlmjVrlsaOHav8/HwtWbJETU1Nwav5Z86cqcGDB6usrEyS9Mgjj2jcuHHKzs7W0aNHtXDhQu3du1d33nlnPN8GAAC2Y7vQnz59uhoaGjR37lzV1dUpNzdXa9euDV7cV1tb2+5CjCNHjmj27Nmqq6tTWlqa8vLy9NFHH2nkyJHxegsAANiS7UJfkoqLi1VcXNzp96qqqtptP/nkk3ryySdjUBUAAN2brT7TBwAA0UPoAwDgELac3gfCYsXT+iSe2AcgYRH6SBxPXWlNPzyxD0CCYnof3VsPr5Rh4dP6JJ7YByBhcaaP7s3lkqYtlfw+yYzwgUk8sQ9AgiP00f25XJLXgucmtDZH3gcA2BjT+wAAOAShDwCAQxD6AAA4BKEPAIBDEPoAADgEoQ8AgEMQ+gAAOAShDwCAQxD6AAA4BCvyAZ3hiX0AEhChD3SGJ/YBSEBM7wNteGIfgATHmT7Qhif2AUhwhD7wVdF4Yh/XBwCwCUIfiDauDwBgE3ymD0QD1wcAsCHO9IFo4PoAADZE6APREo3rAwAgAkzvAwDgEIQ+AAAOQegDAOAQhD4AAA5B6AMA4BCEPgAADsEte0B3YsWSvoFA5GsHAOiWCH2gO7FgSV+PpK+lj5YuWx95PQC6Fab3AbuLwpK+KYe2Ss1HLe0TgP1xpg/YHUv6ArAIoQ90ByzpC8ACTO8DAOAQhD4AAA5B6AMA4BCEPgAADmHL0F+2bJmysrLk9XpVUFCgjRs3dum4VatWyeVy6aabbopugQAAdEO2C/3Vq1erpKRE8+bNU3V1tXJycjR16lQdPHjwnMft2bNHDzzwgCZOnBijSgEA6F5sd8teeXm5Zs+ercLCQknSihUr9NZbb6miokIPP/xwp8cEAgHdfvvtmj9/vt5//30dPXr0rP37/X75/f7gts/nO92HYSpgsDRpLAUMU4bJuMeUYcpz5mXAME4vyYuYCAQCMgxDAcY85hJ97EN5X7YK/ZaWFm3atEmlpaXBfW63W1OmTNGGDRvOetwjjzyiAQMGqKioSO+///45f0ZZWZnmz5/fYf8Xh04oteVY+MUjZIYpHT7uV418crviXY0zuE6d1PAzr3fv2StXr6PxLMdRDMPQ4cOHVVNTI7fbdpOsCS3Rx76xsbHLbW0V+ocOHVIgEFBGRka7/RkZGdqxY0enx3zwwQd64YUXtGXLli79jNLSUpWUlAS3fT6fMjMzNTQ9RWkXWrD4CbosYJiqkansganykPqx0ZoUfDks6xJ5Ui+KYzHOEggEVFNTo+zsbHk8nvMfAMsk+ti3zVh3ha1CP1THjx/XHXfcoeeee07p6eldOiY5OVnJyckd9nvcLoInDtwuF2MfS18ZZ4/bnZD/AO3MfWbMGffYS+SxD+U92Sr009PT5fF4VF9f325/fX29Bg4c2KH97t27tWfPHt1www3BfcaZ9cR79OihnTt3atiwYdEtGgCAbsJWH24kJSUpLy9PlZWVwX2GYaiyslLjx4/v0H7EiBH69NNPtWXLluDXtGnTdPXVV2vLli3KzMyMZfkAANiarc70JamkpESzZs3S2LFjlZ+fryVLlqipqSl4Nf/MmTM1ePBglZWVyev16oorrmh3/AUXXCBJHfYDAOB0tgv96dOnq6GhQXPnzlVdXZ1yc3O1du3a4MV9tbW1CXn1JQAA0Wa70Jek4uJiFRcXd/q9qqqqcx67cuVK6wsCACABcMoMAIBDhHSm/+abb4b8A6655hr16tUr5OMAAIC1Qgr9UB9k43K5tGvXLg0dOjSk4wAAgPVCnt6vq6uTYRhd+kpJSYlGzQAAIAwhhf6sWbNCmqr/x3/8R6WmpoZcFAAAsF5I0/svvvhiSJ0vX748pPYAACB6LLl63zRNmSaPRwUAwM4iCv0XXnhBV1xxhbxeb3B1vOeff96q2gAAgIXCXpxn7ty5Ki8v1w9/+MPguvgbNmzQ/fffr9raWj3yyCOWFQkAACIXdugvX75czz33nG677bbgvmnTpmn06NH64Q9/SOgDAGAzYU/vt7a2auzYsR325+Xl6dSpUxEVBQAArBd26N9xxx2dXp3/85//XLfffntERQGIgdYTUktT5F9cxAt0GyFN75eUlARfu1wuPf/883rnnXc0btw4SdLvf/971dbWaubMmdZWCcBynqU51nQ0JF8qekdyuazpD0DUhBT6mzdvbredl5cnSdq9e7ckKT09Xenp6dq+fbtF5QGwVA+vzAGj5Dpo4d/o/o1S8xEppb91fQKIipBCf/369dGqA0AsuFwybnhaNbUHlD0gRR53BGfnp05Kr844/dowrKkPQFSFffX+V7UtzONieg+wP5dLZlJvqVc/KZLQb222riYAMcHiPAAAOASL8wAA4BAszgMAgEOwOA8AAA7B4jwAADhERFfvv/DCC2ddnOerC/mUl5dHViUAAIhY2KG/bds2jRkzRlLHxXm2bdsWbMdtfAAA2EPYoc9CPQAAdC8hfaa/detWGSGsvLV9+3Yu6gMAwCZCCv2rrrpKf/7zn7vcfvz48aqtrQ25KAAAYL2QpvdN09RPfvITpaSkdKl9S0tLWEUBAADrhRT63/72t7Vz584utx8/frx69eoVclEAAMB6IYV+VVVVlMoAAADRZslT9gA4XOsJqaUp8n56pkjc5gtEDaEPIHJPXWlNP0PypaJ3CH4gSiJ6tC4AB+vhlTJGWdvn/o1S8xFr+wQQxJk+gPC4XNK0pZLfJ5lmZH2dOim9OuP06xDWAgEQGkIfQPhcLsnbL/J+Wpsj7wPAeTG9DwCAQ4R0pn/ppZeG9QCdOXPm6N577w35OAAAYJ2QQn/lypVh/ZCsrKywjgMAANYJKfQnTZoUrToAAECU8Zk+AAAOEXHob9u2TRUVFaqurm63v7GxMewH7ixbtkxZWVnyer0qKCjQxo0bz9r2tdde09ixY3XBBReod+/eys3N1UsvvRTWzwUAIJFFHPrLly/Xhg0b9L//+7+6/fbbVV5erubmZvn9fv3TP/1TyP2tXr1aJSUlmjdvnqqrq5WTk6OpU6fq4MGDnbbv37+/fvzjH2vDhg3aunWrCgsLVVhYqLfffjvStwYAQEKJOPQXLFigr3/96/rggw909OhRvfnmmxoxYoQef/zxsM70y8vLNXv2bBUWFmrkyJFasWKFUlJSVFFR0Wn7yZMn6+abb9bll1+uYcOG6b777tPo0aP1wQcfRPrWAABIKBEvztOnTx89+OCDevDBB+X3+7Vr1y41NDTowIEDId/e19LSok2bNqm0tDS4z+12a8qUKdqwYcN5jzdNU+vWrdPOnTu1YMGCTtv4/X75/f7gts/nkyQFDFMBI8JVxRCSgGHKMBn3WLPluBumPGdeBgKGFAjEtZxoCAQCMgxDgQR8b3aX6GMfyvuKOPS3bdumjRs3Kjc3V2PGjNEVV1wh6fRn+jfffHNIfR06dEiBQEAZGRnt9mdkZGjHjh1nPe7YsWMaPHiw/H6/PB6PnnnmGV1zzTWdti0rK9P8+fM77P/i0AmlthwLqV5ExjClw8f9qpFPbp6vEjN2HHfXqZMafuZ1zZ49MpMPx7WeaDAMQ4cPH1ZNTY3cbq6hjqVEH/vGxsYut4049JcvX66WlhYdOHBAixcvVl5enu666y75/X7dd999evnllyP9EefVt29fbdmyRY2NjaqsrFRJSYmGDh2qyZMnd2hbWlqqkpKS4LbP51NmZqaGpqco7UILlhNFlwUMUzUylT0wVR67pI8D2HLcW5OCL7OzsqQ+F8avligJBAKqqalRdna2PB7P+Q+AZRJ97NtmrLsi4tBfsGCBli9frnXr1sntduvNN9/UU089pVtvvTXkz/TT09Pl8XhUX1/fbn99fb0GDhx41uPcbreys7MlSbm5ufr8889VVlbWaegnJycrOTm5w36P22Wff4AO4na5GPs4sN24f6UOj3FSCpyMvM+eKbZ7RK/b7ZbH40nI4LG7RB77UN5T2KG/b98+ZWZmWvqZflJSkvLy8lRZWambbrpJ0ulpmcrKShUXF3e5H8Mw2n1uD6AbeepKa/oZki8VvWO74AfiKezQHzFihP71X/9VDz/8sFJSUiSdPotu+0xfkmbMmBFyvyUlJZo1a5bGjh2r/Px8LVmyRE1NTSosLJQkzZw5U4MHD1ZZWZmk05/Rjx07VsOGDZPf79dvfvMbvfTSS1q+fHm4bw1ArPXwShmjpPrt1vW5f6PUfERK6W9dn0A3F3bov/vuu7r//vv1wgsv6Kc//am+//3vd2gTzgUT06dPV0NDg+bOnau6ujrl5uZq7dq1wYv7amtr2/Xb1NSku+++W/v371evXr00YsQIvfzyy5o+fXq4bw1ArLlc0rSlkt8nmRHeVXDqpPTqmRMOw4i8NiCBuEwzsr+wX/7yl/rxj3+sAQMGaMmSJZo4caJVtcWEz+dTv379dHjrO0q7MD3e5ThKwDC168AxXXZxP/t8tuwACT/urc3Si9edfv3AbqmPPf6uA4GAdu3apcsuuywhP1e2s0Qf+7YcO3bsmFJTU8/ZNuJ7F2bOnKmdO3fq+uuv13XXXadbb71Vf/zjHyPtFgAAWMyyGxavvfZa3XnnnXr99dc1cuRIPfTQQyHdOwgAAKIr7M/0V6xYoU8++USffPKJPv/8c7ndbl1xxRX6l3/5F+Xk5GjVqlUaOXJk8IE4AAAgvsIO/Z/+9KcqKCjQzJkzNW7cOOXl5alXr17B7//gBz/Q448/ru9///vatm2bJcUCAIDwRXSf/vkUFRXpJz/5Sbg/AgAAWCiqixBnZGRo3bp10fwRAACgi0I607/00ktDXmVPkubMmaN777035OMAAIB1Qgr9lStXhvVDsrKywjoOAABYJ6TQnzRpUrTqAAAAUZZ4DxYGAACdCulM/6vPoT+f8vLykIsBAADRE1Lob968ud12dXW1Tp06peHDh0uS/vCHP8jj8SgvL8+6CgEAgCVCCv3169cHX5eXl6tv3776xS9+obS0NEnSkSNHVFhY2O0eugMAgBOE/Zn+4sWLVVZWFgx8SUpLS9Njjz2mxYsXW1IcAACwTtih7/P51NDQ0GF/Q0ODjh8/HlFRAADAemGH/s0336zCwkK99tpr2r9/v/bv3681a9aoqKhI3/3ud62sEQAAWCCip+w98MAD+t73vqfW1tbTnfXooaKiIi1cuNCyAgEgbK0npJamyPromSKFsRIpYEdhh35KSoqeeeYZLVy4ULt375YkDRs2TL1797asOACIyFNXRt7HkHyp6B2CHwkh5On9uXPnatOmTcHt3r17a/To0Ro9ejSBDyD+eniljFHW9bd/o9R8xLr+gDgK+Ux///79uu6665SUlKQbbrhB06ZN09/+7d8qKSkpGvUBQGhcLmnaUsnvk0wz/H5OnZRenXH6tWFYUxsQZyGHfkVFhQzD0Icffqj//u//1pw5c3TgwAFdc801uvHGG/Wd73xH/fv3j0atANA1Lpfk7RdZH63N1tQC2EhYV++73W5NnDhRP/vZz7Rz5079/ve/V0FBgZ599lkNGjRI3/72t7Vo0SL96U9/srpeAAAQprBv2WtsbAy+vvzyy/XQQw/pww8/1L59+zRr1iy9//77evXVVy0pEgAARC7sq/f79eunX/3qV7rlllva7b/oootUVFSkoqKiiIsDAADWCftM3zRNPfvss/rmN7+pb33rW5ozZ44++eQTK2sDAAAWCjv0pdNP3RszZoy+9a1vafv27Zo4caIeeOABq2oDAAAWCnt6X5JeeeUVXXPNNcHtrVu36sYbb9TgwYN1//33R1wcAACwTthn+v3791dmZma7faNHj9bSpUu1fPnyiAsDAADWCjv0c3Nz9eKLL3bYn52drdra2oiKAgAA1gt7ev+xxx7T1VdfrS+//FJ33323Ro8eraamJj3++OO69NJLrawRAABYIOzQHzdunD7++GPde++9mjhxoswzy116vV79+te/tqxAAABgjYgu5MvJydF7772ngwcPatOmTTIMQwUFBUpPT7eqPgAAYJGQQr+kpOS8bSorKyVJ5eXl4VUEAACiIqTQ37x5c7vt6upqnTp1SsOHD5ck/eEPf5DH41FeXp51FQIAAEuEFPrr168Pvi4vL1ffvn31i1/8QmlpaZKkI0eOqLCwUBMnTrS2SgAAELGwb9lbvHixysrKgoEvSWlpaXrssce0ePFiS4oDAADWCTv0fT6fGhoaOuxvaGjQ8ePHIyoKAABYL+zQv/nmm1VYWKjXXntN+/fv1/79+7VmzRoVFRXpu9/9rpU1AgAAC4R9y96KFSv0wAMP6Hvf+55aW1tPd9ajh4qKirRw4ULLCgQAANYI+0w/JSVFzzzzjP785z9r8+bN2rx5sw4fPqxnnnlGvXv3jqioZcuWKSsrS16vVwUFBdq4ceNZ2z733HOaOHGi0tLSlJaWpilTppyzPQCErPWE1NIU+deZRcyAeIlocR5J6t27t0aPHm1FLZKk1atXq6SkRCtWrFBBQYGWLFmiqVOnaufOnRowYECH9lVVVbrttts0YcIEeb1eLViwQNdee622b9+uwYMHW1YXAAd76sqIu/BI+lr6aOmy9edtC0RL2Gf60VJeXq7Zs2ersLBQI0eO1IoVK5SSkqKKiopO2//nf/6n7r77buXm5mrEiBF6/vnnZRhGcJEgAAhLD6+UMcrSLlMObZWaj1raJxCKiM/0rdTS0qJNmzaptLQ0uM/tdmvKlCnasGFDl/o4ceKEWltb1b9//06/7/f75ff7g9s+n0+SFDBMBQym3mIpYJgyTMY91hj3EHznacl/XDKNyPo5dVKe1bdJkgKnTkmBgAXFoasCgYAMw1AgQcc9lPdlq9A/dOiQAoGAMjIy2u3PyMjQjh07utTHv/3bv2nQoEGaMmVKp98vKyvT/PnzO+z/4tAJpbYcC71ohM0wpcPH/aqRT25XvKtxDsY9HJENlOuUS8PPvN69Z69cvY5GXBG6zjAMHT58WDU1NXK7bTfBHbHGxsYut7VV6EfqiSee0KpVq1RVVSWv19tpm9LS0nbPEPD5fMrMzNTQ9BSlXdgvVqVCp884a2Qqe2CqPKRPzDDucdCaFHw5LOsSeVIvimMxzhMIBFRTU6Ps7Gx5PJ54l2O5thnrrrBV6Kenp8vj8ai+vr7d/vr6eg0cOPCcxy5atEhPPPGEfvvb357zwsLk5GQlJyd32O9xu/gHGAdul4uxjwPGPca+Ms4etzshg8fu3GfGPRHHPpT3ZKt5jqSkJOXl5bW7CK/torzx48ef9bif/exnevTRR7V27VqNHTs2FqUCANDt2OpMXzr9+N5Zs2Zp7Nixys/P15IlS9TU1KTCwkJJ0syZMzV48GCVlZVJkhYsWKC5c+fqlVdeUVZWlurq6iRJffr0UZ8+feL2PgAAsBvbhf706dPV0NCguXPnqq6uTrm5uVq7dm3w4r7a2tp2F2IsX75cLS0tuvXWW9v1M2/ePP3Hf/xHLEsHAMDWbBf6klRcXKzi4uJOv1dVVdVue8+ePdEvCACs0ra6X6R6pkgurslAaGwZ+gCQqDxLc6zpaEi+VPQOwY+Q2OpCPgBISD28MgdYu7qf9m+Umo9Y2ycSHmf6ABBtLpeMG55WTe0BZQ9IiexWyVMnpVdnnH5tRLhSIByH0AeAWHC5ZCb1lnr1U0RLIbY2W1cTHIfpfQAAHILQBwDAIQh9AAAcgtAHAMAhCH0AAByC0AcAwCEIfQAAHILQBwDAIQh9AAAcgtAHAMAhCH0AAByC0AcAwCEIfQAAHILQBwDAIQh9AAAcgtAHAMAhesS7AABAmFpPSC1NkffTM0VyuSLvB7ZH6ANAd/XUldb0MyRfKnqH4HcApvcBoDvp4ZUyRlnb5/6NUvMRa/uELXGmDwDdicslTVsq+X2SaUbW16mT0qszTr82jMhrg+0R+gDQ3bhckrdf5P20NkfeB7oVpvcBAHAIQh8AAIcg9AEAcAhCHwAAhyD0AQBwCEIfAACH4JY9AABL+joEoQ8AYElfh2B6HwCciiV9HYczfQBwKpb0dRxCHwCcjCV9HYXpfQAAHILQBwDAIQh9AAAcwnahv2zZMmVlZcnr9aqgoEAbN248a9vt27frlltuUVZWllwul5YsWRK7QgEA6GZsFfqrV69WSUmJ5s2bp+rqauXk5Gjq1Kk6ePBgp+1PnDihoUOH6oknntDAgQNjXC0AAN2LrUK/vLxcs2fPVmFhoUaOHKkVK1YoJSVFFRUVnbb/xje+oYULF2rGjBlKTk6OcbUAAHQvtrllr6WlRZs2bVJpaWlwn9vt1pQpU7RhwwbLfo7f75ff7w9u+3w+SVLAMBUwIrxPFSEJGKYMk3GPNcY9PhJ+3A1TnjMvAwFDCgTiWs5XBQIBGYahgI1qslIo78s2oX/o0CEFAgFlZGS025+RkaEdO3ZY9nPKyso0f/78Dvu/OHRCqS3HLPs5OD/DlA4f96tGPrlZsTNmGPf4SPRxd506qeFnXtfs2SMz+XBc6/kqwzB0+PBh1dTUyO221QS3JRobG7vc1jahHyulpaUqKSkJbvt8PmVmZmpoeorSLrRggQp0WcAwVSNT2QNT5UnE/4I2xbjHR8KPe2tS8GV2VpbU58L41fJXAoGAampqlJ2dLY/Hc/4Dupm2GeuusE3op6eny+PxqL6+vt3++vp6Sy/SS05O7vTzf4/blZh/iDbndrkY+zhg3OMjocf9K+/J43FLNgtXt9stj8eTkKEfynuyzTxHUlKS8vLyVFlZGdxnGIYqKys1fvz4OFYGAAhJ22N6I/2K9HkA6MA2Z/qSVFJSolmzZmns2LHKz8/XkiVL1NTUpMLCQknSzJkzNXjwYJWVlUk6ffHfZ599Fnz9pz/9SVu2bFGfPn2UnZ0dt/cBAI7GY3pty1ahP336dDU0NGju3Lmqq6tTbm6u1q5dG7y4r7a2tt1FGF9++aWuuuqq4PaiRYu0aNEiTZo0SVVVVbEuHwCcq+0xvfXbreuz7TG9Kf2t69PhbBX6klRcXKzi4uJOv/fXQZ6VlSWT6R8AiD8e09st2C70AQDdFI/ptT3bXMgHAACii9AHAMAhCH0AAByC0AcAwCG4kA8AYF9tC/1EIhBgoZ8zCH0AgH1ZsNCPR9LX0kdLl62PvJ5ujul9AIC9tC30Y6GUQ1ul5qOW9tkdcaYPALAXFvqJGkIfAGA/LPQTFUzvAwDgEIQ+AAAOQegDAOAQhD4AAA5B6AMA4BBcvQ8AcAYrVveTpJ4pp+8u6IYIfQCAI3iW5ljT0ZB8qeidbhn8TO8DABJXD6/MAdau7qf9G6XmI9b2GSOc6QMAEpfLJeOGp1VTe0DZA1LkcUdwdp4Aq/sR+gCAxOZyyUzqLfXqJ0US+gmwuh/T+wAAOAShDwCAQxD6AAA4BKEPAIBDEPoAADgEV+8DABAqK1b3i8PKfoQ+AACheurKyPuIw8p+TO8DANAVPbxShoWr+8VhZT/O9AEA6AqXS5q2VPL7JNMMv584ruxH6AMA0FUul+TtF1kfcVzZj+l9AAAcgjN9AADixYq7AEI4ntAHACBerLgLwN/16wuY3gcAIJasvgsglB8dl58KAIBTWXUXQJujR6QnbulSU0IfAIBYs+IugDYnA11uyvQ+AAAOQegDAOAQhD4AAA5B6AMA4BC2DP1ly5YpKytLXq9XBQUF2rhx4znb//rXv9aIESPk9Xp15ZVX6je/+U2MKgUAoPuwXeivXr1aJSUlmjdvnqqrq5WTk6OpU6fq4MGDnbb/6KOPdNttt6moqEibN2/WTTfdpJtuuknbtm2LceUAANib7UK/vLxcs2fPVmFhoUaOHKkVK1YoJSVFFRUVnbZ/6qmn9Hd/93d68MEHdfnll+vRRx/VmDFjtHTp0hhXDgCAvdnqPv2WlhZt2rRJpaWlwX1ut1tTpkzRhg0bOj1mw4YNKikpabdv6tSpeuONNzpt7/f75ff7g9vHjh2TJB09EttnGkMKGKZ8viYdSWqVx+2KdzmOwbjHB+MeP4k+9r6jRyVJZhcW+rFV6B86dEiBQEAZGRnt9mdkZGjHjh2dHlNXV9dp+7q6uk7bl5WVaf78+R32D500PcyqAQCIv+PHj6tfv3Mv+GOr0I+F0tLSdjMDR48e1SWXXKLa2trzDhas5fP5lJmZqX379ik1NTXe5TgG4x4fjHv8JPrYm6ap48ePa9CgQedta6vQT09Pl8fjUX19fbv99fX1GjhwYKfHDBw4MKT2ycnJSk5O7rC/X79+CfnL0B2kpqYy9nHAuMcH4x4/iTz2XT1ptdWFfElJScrLy1NlZWVwn2EYqqys1Pjx4zs9Zvz48e3aS9K777571vYAADiVrc70JamkpESzZs3S2LFjlZ+fryVLlqipqUmFhYWSpJkzZ2rw4MEqKyuTJN13332aNGmSFi9erOuvv16rVq3S//3f/+nnP/95PN8GAAC2Y7vQnz59uhoaGjR37lzV1dUpNzdXa9euDV6sV1tbK7f7LxMUEyZM0CuvvKJ///d/149+9CNddtlleuONN3TFFVd06eclJydr3rx5nU75I7oY+/hg3OODcY8fxv4vXGZXrvEHAADdnq0+0wcAANFD6AMA4BCEPgAADkHoAwDgEI4IfR7VGz+hjP3KlSvlcrnafXm93hhWmxh+97vf6YYbbtCgQYPkcrnO+hyKr6qqqtKYMWOUnJys7OxsrVy5Mup1JppQx72qqqrD77vL5TrrEuLoXFlZmb7xjW+ob9++GjBggG666Sbt3LnzvMc59f98woc+j+qNn1DHXjq9YtaBAweCX3v37o1hxYmhqalJOTk5WrZsWZfa//GPf9T111+vq6++Wlu2bNGcOXN055136u23345ypYkl1HFvs3Pnzna/8wMGDIhShYnpvffe0z333KOPP/5Y7777rlpbW3XttdeqqanprMc4+v+8meDy8/PNe+65J7gdCATMQYMGmWVlZZ22/4d/+Afz+uuvb7evoKDA/Od//ueo1pmIQh37F1980ezXr1+MqnMGSebrr79+zjYPPfSQOWrUqHb7pk+fbk6dOjWKlSW2roz7+vXrTUnmkSNHYlKTUxw8eNCUZL733ntnbePk//MJfabf9qjeKVOmBPd15VG9X20vnX5U79nao3PhjL0kNTY26pJLLlFmZqZuvPFGbd++PRblOhq/8/GVm5uriy++WNdcc40+/PDDeJfT7bU9Lr1///5nbePk3/mEDv1zPar3bJ+bhfqoXnQunLEfPny4Kioq9F//9V96+eWXZRiGJkyYoP3798eiZMc62++8z+dTc3NznKpKfBdffLFWrFihNWvWaM2aNcrMzNTkyZNVXV0d79K6LcMwNGfOHH3zm98856qsTv4/b7tleOFc48ePb/egpAkTJujyyy/Xs88+q0cffTSOlQHWGz58uIYPHx7cnjBhgnbv3q0nn3xSL730Uhwr677uuecebdu2TR988EG8S7GthD7Tj8WjetG5cMb+r/Xs2VNXXXWVampqolEizjjb73xqaqp69eoVp6qcKT8/n9/3MBUXF+t//ud/tH79eg0ZMuScbZ38fz6hQ59H9cZPOGP/1wKBgD799FNdfPHF0SoT4nfeTrZs2cLve4hM01RxcbFef/11rVu3Tpdeeul5j3H073y8rySMtlWrVpnJycnmypUrzc8++8z8wQ9+YF5wwQVmXV2daZqmeccdd5gPP/xwsP2HH35o9ujRw1y0aJH5+eefm/PmzTN79uxpfvrpp/F6C91WqGM/f/588+233zZ3795tbtq0yZwxY4bp9XrN7du3x+stdEvHjx83N2/ebG7evNmUZJaXl5ubN2829+7da5qmaT788MPmHXfcEWz/xRdfmCkpKeaDDz5ofv755+ayZctMj8djrl27Nl5voVsKddyffPJJ84033jB37dplfvrpp+Z9991nut1u87e//W283kK3dNddd5n9+vUzq6qqzAMHDgS/Tpw4EWzD//m/SPjQN03TfPrpp82vfe1rZlJSkpmfn29+/PHHwe9NmjTJnDVrVrv2v/rVr8yvf/3rZlJSkjlq1CjzrbfeinHFiSOUsZ8zZ06wbUZGhvn3f//3ZnV1dRyq7t7abgX766+2sZ41a5Y5adKkDsfk5uaaSUlJ5tChQ80XX3wx5nV3d6GO+4IFC8xhw4aZXq/X7N+/vzl58mRz3bp18Sm+G+tszCW1+x3m//xf8GhdAAAcIqE/0wcAAH9B6AMA4BCEPgAADkHoAwDgEIQ+AAAOQegDAOAQhD4AAA5B6AMA4BCEPgAADkHoAwDgEIQ+gLBUVVUpKyvL9n0C+AtCHwAAhyD0AQBwCEIfgCWGDBmiZ555pt2+jz76SCkpKdq7d2+cqgLwVYQ+AEsUFBTok08+CW6bpqk5c+bo/vvv1yWXXBLHygC0IfQBWGLcuHHtQv+ll17Svn37VFpaKkk6fvy47r77bm3durXTbQDRR+gDsMS4ceP0+eefq7GxUU1NTfrRj36kxx57TH369JEkLV++XH6/Xxs3bux0G0D0EfoALJGXlye3263q6motWLBAF110kQoLC4PfX7dunbKyspSbm9vpNoDoI/QBWCIlJUVXXnml1qxZo0WLFunJJ5+U2336X8zJkyfldrv12WefKS8vr8M2gNgg9AFYZty4cXr66ac1depUTZ48Obh/165dampqUl5enlwuV4dtALFB6AOwTE5Ojnr27KmFCxe229/Q0KAvvvhCd911V6fbAGKjR7wLAJA4Vq1apeLiYmVnZ7fbf+DAAd16661qbm6WYRgdtvv27RunigFn4UwfQEQMw1B9fb0ef/xx7dq1S/PmzWv3/UAgoOrqau3bt0/33ntv8GK/tu2ePXvGqXLAeTjTBxCR3/3ud/qbv/kbjRgxQmvWrFFqamq773s8Hi1evLjdvr/eBhAbhD6AsGRlZWnOnDmaPHmyDMOwtE8A0eEyTdOMdxEAACD6+EwfAACHIPQBAHAIQh8AAIcg9AEAcAhCHwAAhyD0AQBwCEIfAACHIPQBAHAIQh8AAIf4/5bmCcvpjkZOAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from matplotlib import pyplot as plt\n", + "\n", + "fig, ax = plt.subplots(figsize=(5.6, 3.9))\n", + "\n", + "# matplotlib's 'step' function requires the last value to be repeated\n", + "nbins = np.append(bins, bins[-1])\n", + "edges = np.arange(0.0, 2.4, 0.1)\n", + "\n", + "ax.step(edges, nbins, where='post', color=\"C1\")\n", + "plt.fill_between(np.arange(0.0, 2.4, 0.1), nbins, step=\"post\", color=\"C1\", alpha=0.2)\n", + "ax.set_xlabel(\"$|y_{\\ell\\ell}|$\")\n", + "ax.set_ylabel(\"$\\mathrm{d} \\sigma / \\mathrm{d} |y_{\\ell\\ell}|$ [pb]\")\n", + "ax.grid(True, alpha=0.5)\n", + "ax.set_ylim(bottom=0.0)\n", + "ax.set_xlim([edges[0], edges[-1]])\n", + "\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "PineAPPL", + "language": "python", + "name": "pineappl" + }, + "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.12.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pineappl_py/docs/source/conf.py b/pineappl_py/docs/source/conf.py index 8c3b56cf..182ce7f3 100644 --- a/pineappl_py/docs/source/conf.py +++ b/pineappl_py/docs/source/conf.py @@ -26,6 +26,7 @@ 'sphinx.ext.napoleon', 'sphinx.ext.todo', 'sphinx_rtd_theme', + 'nbsphinx', ] diff --git a/pineappl_py/docs/source/index.rst b/pineappl_py/docs/source/index.rst index 600c36b3..2b26d767 100644 --- a/pineappl_py/docs/source/index.rst +++ b/pineappl_py/docs/source/index.rst @@ -1,22 +1,26 @@ Welcome to PineAPPL =================== -This is the Python wrapper for the `Rust PineAPPL library `_ using `PyO3 `_. +This is the Python wrapper for the `Rust PineAPPL library `_ +using `PyO3 `_. -PineAPPL is a computer library that makes it possible to produce fast-interpolation grids for fitting parton distribution functions (PDFs) including corrections of strong and electroweak origin. +PineAPPL is a computer library that makes it possible to produce fast-interpolation +grids for fitting parton distribution functions (PDFs) including corrections of +strong and electroweak origin. -The :doc:`installation` instructions are given :doc:`here `. +The :doc:`installation` instructions are given :doc:`here ` and tutorials +are available :doc:`here `. In addition, some practical examples can be found +in the ``example/`` subfolder of the `repository `_. -A practical example can be found in the ``example/`` subfolder of the `repository `_. -The Python wrapper is also used in `yadism `_ and `pineko `_. -We also list some common :doc:`recipes` here. +The Python wrapper is also used in `yadism `_ and +`pineko `_. .. toctree:: :maxdepth: 1 :hidden: :caption: Contents: - installation - recipes + Installation + Tutorials API indices diff --git a/pineappl_py/docs/source/introduction.ipynb b/pineappl_py/docs/source/introduction.ipynb new file mode 100644 index 00000000..63808c6c --- /dev/null +++ b/pineappl_py/docs/source/introduction.ipynb @@ -0,0 +1,732 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1bc5cec5-dd52-4dbf-8e0d-86cb3996292c", + "metadata": {}, + "source": [ + "# A Basic Introduction\n", + "\n", + "In the following introduction, we showcase some common use cases of `PineAPPL`\n", + "that are available through the Python interface. For more details about all\n", + "the available functionalities, refer to the [API documentation](https://pineappl.readthedocs.io/en/latest/modules/pineappl.html#)." + ] + }, + { + "cell_type": "markdown", + "id": "fb1fbb97-5222-40e5-b2e4-e737212f4c5b", + "metadata": {}, + "source": [ + "## Setting up\n", + "\n", + "In order to start using `PineAPPL`, one first needs to install it. \n", + "The easiest way to install the latest stable version is via `pip` \n", + "using the following command:" + ] + }, + { + "cell_type": "raw", + "id": "ecb04362-1f02-4d0b-81f6-df1b95134c67", + "metadata": {}, + "source": [ + "pip install pineappl" + ] + }, + { + "cell_type": "markdown", + "id": "318119e2-2be1-48c1-8448-8b4b12a08191", + "metadata": {}, + "source": [ + "In order to check that `PineAPPL` was installed properly, one can try to import it:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "5b870d1b-1325-4903-b885-c528a42c869c", + "metadata": {}, + "outputs": [], + "source": [ + "import pineappl" + ] + }, + { + "cell_type": "markdown", + "id": "5e3411d9-fcca-438f-8717-64c509461711", + "metadata": {}, + "source": [ + "If the above command was successful (did not result in an error), we can\n", + "download the `PineAPPL` grid that will be used throughout this tutorial." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d67f35c1-2ad6-4a5a-a651-2f271258875d", + "metadata": {}, + "outputs": [], + "source": [ + "!wget --no-verbose --no-clobber 'https://data.nnpdf.science/pineappl/test-data/LHCB_DY_8TEV.pineappl.lz4'" + ] + }, + { + "cell_type": "markdown", + "id": "767c24c7-e1f6-4ec4-ac42-25badf7dc363", + "metadata": {}, + "source": [ + "We can now load the downloaded grid by running the following:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "7ed7bdf0-d28a-4830-ac56-1fa9c7f81d56", + "metadata": {}, + "outputs": [], + "source": [ + "# Load the PineAPPL grid\n", + "grid = pineappl.grid.Grid.read(\"./LHCB_DY_8TEV.pineappl.lz4\")" + ] + }, + { + "cell_type": "markdown", + "id": "d0449d9b-26fc-4847-9689-5a4361cfef8b", + "metadata": {}, + "source": [ + "**NOTE:** For the `pineappl.grid.Grid.read()` function, both `.pineappl` and `.pineappl.lz4` \n", + "extensions are acceptable, as long as they are consistent. Without the `.lz4` extension, the\n", + "grid is assumed to not be compressed." + ] + }, + { + "cell_type": "markdown", + "id": "097137f5-4f46-4d20-af7d-ff5ad3c31ced", + "metadata": {}, + "source": [ + "## How can I convolve a given PineAPPL grid with a PDF?" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "fc2f1cd5-28d3-43db-8991-557287780281", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "LHAPDF 6.5.4 loading /home/felix/local/share/LHAPDF/NNPDF40_nnlo_as_01180/NNPDF40_nnlo_as_01180_0000.dat\n", + "NNPDF40_nnlo_as_01180 PDF set, member #0, version 1; LHAPDF ID = 331100\n" + ] + } + ], + "source": [ + "# We first need to load the PDF set with LHAPDF\n", + "import lhapdf\n", + "import numpy as np\n", + "# `Polars` is a better alternative to Pandas (written in Rust!)\n", + "import polars as pl\n", + "\n", + "# Choose the PDF set to be used\n", + "pdfset = \"NNPDF40_nnlo_as_01180\"\n", + "# Use the central PDF, ie. member ID=0\n", + "pdf = lhapdf.mkPDF(pdfset, 0)" + ] + }, + { + "cell_type": "markdown", + "id": "640a1efd-94bb-4c22-a6d7-785e940a0013", + "metadata": {}, + "source": [ + "Our grid can now be convolved with our PDF set using the `convolve_with_one()` function:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "d41e391c-affd-49ad-8d70-483bf952d4f3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "shape: (18, 2)
binspredictions
i64f64
08.159231
123.890855
238.306994
351.278627
462.349204
……
1322.640245
1412.927123
156.138797
161.435175
170.052598
" + ], + "text/plain": [ + "shape: (18, 2)\n", + "┌──────┬─────────────┐\n", + "│ bins ┆ predictions │\n", + "│ --- ┆ --- │\n", + "│ i64 ┆ f64 │\n", + "╞══════╪═════════════╡\n", + "│ 0 ┆ 8.159231 │\n", + "│ 1 ┆ 23.890855 │\n", + "│ 2 ┆ 38.306994 │\n", + "│ 3 ┆ 51.278627 │\n", + "│ 4 ┆ 62.349204 │\n", + "│ … ┆ … │\n", + "│ 13 ┆ 22.640245 │\n", + "│ 14 ┆ 12.927123 │\n", + "│ 15 ┆ 6.138797 │\n", + "│ 16 ┆ 1.435175 │\n", + "│ 17 ┆ 0.052598 │\n", + "└──────┴─────────────┘" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "predictions = grid.convolve_with_one(2212, pdf.xfxQ2, pdf.alphasQ2)\n", + "df_preds = pl.DataFrame({\n", + " \"bins\": range(predictions.size),\n", + " \"predictions\": predictions,\n", + "})\n", + "df_preds" + ] + }, + { + "cell_type": "markdown", + "id": "231dd694-e068-4e62-8e19-b93c96f4d937", + "metadata": {}, + "source": [ + "The length of `predictions` is the same as the number of bins that define the\n", + "observable in question. In our example, it is the inclusive cross-section for\n", + "Z boson production in bins of the muon pseudorapidity at 8 TeV. See \n", + "[arXiv:1511.08039](https://arxiv.org/abs/1511.08039) for more details, with the\n", + "corresponding [HepData](https://www.hepdata.net/record/ins1406555) tables.\n", + "\n", + "To see that our predictions are consistent with the experimental measurements,\n", + "we can plot a comparison of the two using **only** the central values." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "fee97bfd-e392-4a78-8850-a4093bbcb330", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAFrCAYAAABbtho0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABC4ElEQVR4nO3dfVhUZf4/8PeZAQblUUNRFHzKUkLBAMlWNzUSsVyt3Fzzl9iDpju0JVrpdxOzttxdlSx3Vsryqb5l61XpblqauD58zRAxKkVYcUkMFVNBZNQBzjm/P2hGkAPMDHPmyffruriGOXPmns89Z+B85j73gyDLsgwiIiKiG2hcHQARERG5JyYJREREpIhJAhERESlikkBERESKmCQQERGRIiYJREREpIhJAhERESnycXUA7kiSJJw+fRpBQUEQBMHV4RAREdlElmVcvnwZERER0Gjsbw9gkqDg9OnTiIyMdHUYRERE7XLq1Cn07NnT7uczSVAQFBQEoOHNDQ4Obnd5oijixIkT6NevH7RabbvLc1esp3dhPb0L6+ld2qpndXU1IiMjLeczezFJUGC+xBAcHOywJCEwMBDBwcFe/6FlPb0H6+ldWE/vYm0923vJnB0XiYiISBGTBCIiIlLEJIGIiIgUsU8CERE5hCiKqKurc3kMkiTh2rVrXt0noT3DGm3BJIHIC5VXXUWlsdbq/TsF+KFHaAcVIyJvJssyzp49i6qqKleHAlmWUV9fj5MnT3r1PDfmesqyrOrreHWScOXKFQwcOBC//e1vsWzZMleHQ9SEWify8qqrGL1sN0z1EgBAh1qM0+RijPYQQlGDKgRih5iAbVISTPBr2MdHg13zRlpdPhMQasycIHTt2hUdO3Z06clZlmWYTCbodDqvTRJkWYbRaMSZM2dQUVGBHj16qPZaXp0kvPbaa7jrrrtcHQZRM2qeyCuNtZZykzX5WOabjVDBCFEWoBVkiLKAVG0eFskbMLduFnKkeJjqJVQaa9ssW+0EhDyPKIqWBOGWW25xdTiWb9b+/v5emyQADfWrq6tDZWUlunXrptqlFa9NEo4fP46ioiKMHz8eR44ccXU4RE2oeSI3S9bk4x3fLAAN/zS1QtPbYBix2jcLM+sysFOKd5u42VLhWcx9EDp27OjiSG4+/v7+ABqOwU2VJOzduxdLly5Ffn4+zpw5g88++wwTJ05sso/BYMDSpUtx9uxZxMbGYuXKlRg6dKjl8Xnz5mHp0qX4+uuvnRw9eROlE5Yoiii7YIKp46Vmf5i2nrDUOJEDDd/wl/lmA5ChaeHLlEYAJFnGMt9sJJkMVpetZtw3tlRYgy0V7sGbv7W7K2e8526ZJBiNRsTGxuKJJ57AQw891Ozxjz/+GBkZGcjOzkZSUhJWrFiBlJQUFBcXo2vXrtiyZQtuu+023HbbbVYlCSaTCSaTyXK/uroaQMPJQBTFdtfH3NvWEWW5M2+r5+mqq0h+Y59i03oYanC6hab1nXNGIKKNE5YoinadyK35TIqiiHGaXIQKxjbrqBGAUBiRqjkIUbynSdlKx1PNuAHgfPXVJgmCNZczTPUSzldfRbcgvzbLV+Jtn9uWqFVPURQhy7Llx1pqtRiZY1C7Q5+rNa6f0t+Xo46zWyYJqampSE1NbfHxrKwszJgxA48//jgAIDs7G1u3bsWaNWswf/58fPPNN9i4cSM2bdqEmpoa1NXVITg4GJmZmYrlLVmyBIsXL262/cSJEwgMDGx3fSRJwsWLF1FSUuK0YSuu4G31PH7BZFfTekHRCRhv0bVadtkFk10n8rJTfaG7cq7NssdoD1nibIsoC0jR5qHs1KkmZSsdTzXjNpdvZu17DqBZ7Lbwts9tS9SqpyRJqK+vb/JFqy2nL11D6t++Qa0NLUZ+Php8kX4XIkL829y3vr7e6nIb27t3L8aOHYvTp08jNDTUrjKcSRRFy0iOG49pTU2NQ17DLZOE1tTW1iI/Px8LFiywbNNoNEhOTsaBAwcANJz0lyxZAgBYt24djhw50mKCAAALFixARkaG5b55YYx+/fo5bO2GkpIS3HrrrV49btfb6mnqeAnATzY3rUdFRqJ/j5A2y7bnRB4R+bxVZV9BjVXlmusRihp0vCFupeOpZtzm8tV6z1uM0cs+ty1Rq57Xrl3DyZMnodPpLNfI23LlgsmmBAEAauslXKkX2nwN8zdsa0Y3jBo1CrGxsVixYgUAwM+voTXK39/f6rq4inkUh4+PD3r16tUsXnOLeHt5XJJw/vx5iKKI8PDwJtvDw8NRVFRkV5k6nQ46XfNvflqt1mF/TBqNxqHluStvqqdWq7Wrad2a+mu1WoTacSK3tuwqBNp0Iq9CIIIUyr7xeKoZt7l8td7z1njT57Y1atRTq9VCEATLjzXsvZZu62tYs2/j/RrfOrOPRW1trSVBsYfSMXXYucshpbix6dOnc44EL1ZedRVHyi/hSPklHD1ZgZ92r8Gldb9DzdspuLTud/hp9xocPVlh2ae86qpN5Zub1ls6WZlpBCBUaGhat5b5RG4N84ncWjvEBJtO5NvFRKvLVjNuQN33nMhs+vTp2LNnD958801LUvDjjz8CAPLz85GQkICOHTvi7rvvRnFxcZPnbtmyBXfeeSf8/f3Rt29fLF68uMkljrKyMkyYMMGySuMjjzyCiooKy+Mvv/wy4uLi8O6776JPnz7w9/fHhg0bcMsttzS7bDNx4kQ89thj6r0RbfC4loSwsDBotdombzgAVFRUoFu3bi6KilyhcU/4lq5fh/z4BQL//ZLl+rWtPeHtaVq31g4xAalW7m8+kT9lZdnbpCQskjcgGK2fbCUZqEYAvpCGWl22mnED6r7nao9WIc/x5ptv4j//+Q9iYmLwyiuvAACOHj0KAPjjH/+I5cuXo0uXLpg1axaeeOIJ7N+/HwCwb98+TJs2DW+99RZGjBiBEydOYObMmQCARYsWQZIkS4KwZ88e1NfXQ6/XY/Lkydi9e7fl9UtKSvDJJ5/g008/hVarRf/+/fGHP/wB//znP/Hb3/4WAHDu3Dls3boVO3bscOI705THJQl+fn6Ij49HTk6OZVikJEnIyclBenp6u8o2GAwwGAxe38vZW5jH7Nt0/bo+3qYx+/Y0rVtLzRO5CX6YWzcLq32zIMnKzfaSDAAC5tbNsowUcHXcgHrveWsTQbU2WoXDK71TSEgI/Pz80LFjR8sXTPMl69deew333HMPAGD+/Pm4//77ce3aNfj7+2Px4sWYP38+0tLSAAB9+/bFq6++ihdeeAGLFi1CTk4OfvjhB5SWliIyMhIAsGHDBtxxxx3Iy8tDYmJDq11tbS02bNiALl26WGJ69NFHsXbtWkuS8MEHHyAqKgojR450ynuixC0vN9TU1KCgoAAFBQUAgNLSUhQUFKCsrAwAkJGRgdWrV2P9+vU4duwYZs+eDaPRaBntYC+9Xo/CwkLk5Vn/zYRcy9rr10DD9WsdrB9yBajbtG4+kQPCLyfs5uw9kQNAjhSPmXUZqEaAJb7Gt9UIwIy6DMsIAXeJW633/MaJoHJ1erzhtwpjNIcwTHsMYzSH8IbfKuTq9LhXkw8Alomg6OYyePBgy+/du3cH0PCtHgC+++47vPLKKwgMDLT8zJgxA2fOnMGVK1dw7NgxREZGWhIEAIiOjkZoaCiOHTtm2darV68mCQIAzJgxAzt27EB5eTmAho7306dPd+kcFG7ZknDo0CGMGjXKct888iAtLQ3r1q3D5MmT8fPPPyMzMxNnz55FXFwcvvzyy2adGcn72TMcDxhtdflqN62bT+TLfLMRiqaXSrSCjGoENBnqZ41OAX7Q+WhgqpewU4pHksmAVM1BpGjzLHMNbBcT8YU0tMk35k4B1p/M1YjbTO33XK2JoMh7+Pr6Wn43n6AlqSHBrKmpweLFixXn8LFlRERAQECzbUOGDEFsbCw2bNiAMWPG4OjRo9i6daut4TuUWyYJI0eObHMijPT09HZfXiDPp+b1a0C9pnU1T+Q9Qjtg17yRN3wDvp4YBQF46pefxvFY06TujAREzcsZas9ESZ7Fz8/P5svLd955J4qLi3HrrbcqPj5w4ECcOnUKp06dsrQmFBYWoqqqCtHR0W2W/9RTT2HFihUoLy9HcnJykxYJV3DLJMFV2CfB86jZZwBQ79q+midyc/lqXEdXO25A3f4Uarc8kWfp3bs3cnNz8eOPPyIwMNDSWtCazMxMPPDAA4iKisKkSZOg0Wjw3Xff4ciRI/jTn/6E5ORkDBo0CFOnTsWKFStQX1+P3//+97jnnnuQkJDQZvmPPvoo5s2bh9WrV2PDhg2OqGa7MEloRK/XQ6/Xo7q6GiEh9k3MQm2ouwYUbgaKPgeuVAIdOwEDHgCiJwK+tk9eYtd8ADa+hlpN62qdyNWmZtzmloqceuvfc1taKtRueSLPMm/ePKSlpSE6OhpXr17F2rVr23xOSkoKPv/8c7zyyiv4y1/+Al9fXwwYMABPPdWQGguCgC1btuCZZ57Br3/9a2g0GowdOxYrV660KqaQkBA8/PDD2Lp1a7M1i1yBSQKpzjzsLOjkDvTYPRc+tZcgQwMBUsPtsX+hfusLKB+Zhcu97rPpm6ea16+d0bROTTVtqRiO8vqncbl0G4J/3A6tqRKirhOqe6egus84zPHxxxzY1lKhdssTeZbbbrvNMlOv2fTp05vcj4uLa3b5OyUlBSkpKS2WGxUVhS1btrT4+Msvv4yXX365xcfLy8sxdepUxUn+nI1JAqnKPOxshJR3vbOYAAhoaNYz32pMlxC5/SnMrMvAPk2i1cPO1Lx+3VLTuiiKKDt1ClGRkXhKq21X0zo117SlIgTo9TiA6yOX2tPG54yWJ2pb4wTcWjdDAl5ZWYndu3dj9+7d+Pvf/+7qcAAwSSCVVRprgfprWKazpbPYIKvnMlDz+jWg3LQuiiJ0V86hf48Qr5/G19uo2fKk1qqG3kg5AW/dzfB+DRkyBJWVlfjLX/6C22+/3dXhAGCSQE6gdmcxNYfjkXdRq+XpxomarHGzT9TkqX1y1GSeFtqdMElohKMb1OGMzmLW9hmgm5taLU+NJ2oCms7maP4s3jibo3miJp4oyZ0xSWiEoxvUoVZnsRuva5rgh83ScGyWhrf4nJvhuia1Tu2Wp5bWEUnV5mGRvIGtWuRRmCSQ6tTqLMbrmmQLZ4xW4WyO5G2YJJDq1OwsxuuaZC21R6twNkfyRkwSSHVqrxpIZC01R6twNkfyRm65CiR5F7VXDSRyB+YOutbgbI7kKZgkNGIwGBAdHW1Z75scR61li4ncBWdzbIe6a8B3G4GP/x+w9v6G2+82Nmx3spEjR+K5555z+uu6K15uaISjG37h4PUVzDhMkbwZZ3O0U9E2YPNs4FoVIGgAWWq4PfYv4IsXgQezgdtTXR2lot27d2PUqFGorKxEaGioq8NRBZMEAqDu+gqNWTNMkcgTqdlB12sVbQM2Pnr9viw1vb12CfhoCvC7D4EB45wfH/FyA12fLW6F4U1Ebp8BjekSgJbXV1hheBOjl+1GedXVNss2DzuzBecyIE+0TUpClRzQYr8bM0kGquSGDro3tbprDS0IAMxDRpv7Zfvm2apcejAajZg2bRoCAwPRvXt3LF++vMnj77//PhISEhAUFIRu3brh0Ucfxblz5wA0zI44atQoAECnTp0gCIJlcagvv/wSw4cPR2hoKG655RY88MADOHHihMPjdwa2JJCq6ytwLgO6Wai9jojXKdzccImhTXLDfoVbgNjJDg3h+eefx549e7BlyxZ07doV//M//4PDhw8jLi4OAFBXV4dXX30Vt99+O86dO4eMjAxMnz4d27ZtQ2RkJD755BM8/PDDKC4uRnBwMDp0aPi/ZTQakZGRgcGDB6OmpgaZmZl48MEHUVBQAI3Gs76bM0kgAOoO3+JcBuTtzC1mOfXWz+Z407eYFX1+vQ9CWwQNUPQvhyYJNTU1eO+99/DBBx/g3nvvBQCsX78ePXv2tOzzxBNPWH7v27cv3nrrLSQmJqKmpgaBgYHo3LkzAKBr165N+iQ8/PDDTV5rzZo16NKlCwoLCxETE+OwOjgDkwQC4Jz1FYi8VdMWs+Eor38al0u3IfjH7dCaKiHqOqG6dwqq+4zDHB9/zIFtLWZeucLklUrrEgSgYb+rlQ59+RMnTqC2thZJSUmWbZ07d26y+mJ+fj5efvllfPfdd6isrIQkNcRbVlaG6OjoFss+fvw4MjMzkZubi/Pnzzd5HpMED3YzL/DE4VtE7dO0xSwE6PU4gMctj9s7XurGFSatWTzKI1aY7NjJtpaEDp3Uj6kRo9GIlJQUpKSk4H//93/RpUsXlJWVISUlBbW1rSds48ePR69evbB69WpERERAkiTExMS0+Tx35FkXR1Sm1+tRWFiIvLyb71uyefiWNczDt4hIfY1XmEzW5CNXp8cbfqswRnMIw7THMEZzCG/4rUKuTo97NfkArq8w6dYGPGBbS8KA8Q59+X79+sHX1xe5ubmWbZWVlfjPf/4DACgqKsKFCxfw5z//GSNGjMCAAQMsnRbN/PwakrLGXywvXLiA4uJivPTSS7j33nsxcOBAVFY6thXEmZgkEICG4Vu2tCRsFznhFJEzmRePCkZD36GWFo9K/iVRcHvREwH/UABtfTkRGvaLnuDQlw8MDMSTTz6J559/Hrt27cKRI0cwffp0S8fCqKgo+Pn5YeXKlfjvf/+Lf/7zn3j11VeblNGrVy8IgoDPP/8cP//8M2pqatCpUyfccssteOedd1BSUoJdu3YhIyPDobE7E5MEAsDhW0TuzNrFo4CG0Uc6uHkrAtAwMduD2b/caSlR+GX7g9ntmsitJUuXLsWIESMwfvx4JCcnY/jw4YiPb5j1tUuXLli3bh02bdqE6Oho/PnPf8ayZcuaPL9Hjx5YvHgx5s+fj/DwcKSnp0Oj0WDjxo3Iz89HTEwM5syZg6VLlzo8dmdhnwQCwOFbRO7MaxePuj21YaIkpRkXZQnwD1F1xsXAwEC8//77eP/99y3bnn/+ecvvU6ZMwZQpU5o8R5abfpNauHAhFi5c2GRbcnIyCgsLW32ep2CSQBbm9RWsGb5FRM7j1aOPBowD5hY3zINQ9K+GUQwdOjX0QYieoEoLAlmPSQI1wfUViNyP148+8vVvmAPBwZMlUfsxSaBmuL4CkXvh4lHkKuy4SFxfgcjNcfQRuQpbEhq5WSdT4voKRO5tm5SERfIGBMPY4ugGoKFzcTUaRh/d9CtMkkMwSWhEr9dDr9ejuroaISH2zo/mmbi+ApH78oTRR+aph8l5nDFigkmCp6q71rCKWtHnDXOgd+zUMINZ9ET2BibyQu46+sjPzw8ajQanT59Gly5d4OfnB0GwbvZWNciyDJPJBAAujUNNsiyjtrYWFRUV0Gg0lpkf1cAkwRMVbVMeV3zsX8AXL6o6rpiInMvcZ8hUL1k9+siZfYY0Gg369OmDM2fO4PTp0055zdbIsoz6+nr4+Ph4bZIANNSzrq4Ot912m6rLTzNJ8DRF24CNj16/b5773Hx77RLw0ZSGCUoGjHN+fETkUMp9hq5PlBQE4Klffsyc3WfIz88PUVFRqK+vd3mfLlEUcfLkSfTq1QtardalsahJEAT897//VbUVAWCS4FHKz1ci/NNZ0AIQ0NK1KBkyBIifzkLFzO/QI8y5K6cRkeN5Qp8hQRDg6+sLX19fl8YhiiI0Gg38/f29OkkQRdEpLSUcAukhyquu4s03/wqf2kutJAgNBMjwqb2EFW8uRXnVVSdFSERE3oYtCR6i0liLUcizaUKV0TiISmOt238DISLXKa+6yuHP1CImCR7E66dmJSKnKq+6itHLdsNUb/3wRZ2PBrvmjWSicJPg5QYPYp6a1RrmqVmJiFpSaay1KUEAAFO9ZFPLA3k2tiR4kB1iAlKtXN3NPDUrZ10jImvpUItxmlyM0R6yDK/cISZgm5TExd1uUmxJ8CDbpCRUyQG/zKzWMkkGquSGqVmJiKyRrMlHrk6PN/xWYYzmEIZpj2GM5hDe8FuFXJ0e92ryXR0iuQCThEYMBgOio6ORmOiei6OYp2YFhBYTBVdPzUpEnidZk493fLMQDCMAWPo+mW+DYcRq3ywkM1G46TBJaESv16OwsBB5edY16buCeWrWagQAgKWPgvm2GgGYUZfh9KlZicgz6VCLZb7ZAJTXhADwy3YZy3yzoQP7I9xM2CfBA1k7NSsRUVvGaXIRKhjb3E8jAKEwIlVzEI1nfCTvxiTBQ5ngh83ScGyWhrs6FCLyYGO0h2yafyXFys7T5B14uYGI6CbG+VeoNUwSiIhuYpx/hVrDJMFDmJeLtYUzl4slIs+0Q0ywqSVhu+ieo79IHeyT4CGUl4ttHedYJ6K2bJOSsEjegGAYWxzdADQMr65Gw/wrnKTt5sEkwYN4wnKxROQ5OgX4AT7+mFs3C6t9syDJysMgG8+/Ah9/tlDeRJgkEBHdpK63UN6NUycHosfuudDUXoIMDQRIlltJF4LykVmY0+s+vMIWypsKkwQiopuYpYWyx2+BxPFA4RYIRf8CrlZC6NAJGDAePtET0MvX39WhkgswSSAioga+/kDs5IYfInB0AxEREbWASQIREREpYpJAREREipgkEBERkSImCY0YDAZER0cjMZEzihERETFJaESv16OwsBB5eVzljIiIiEMgiYhINeVVVzmdvAdjkkBERKoor7qK0ct2w1QvAQB0qMU4TS7GaA8hFDWoQiB2iAnYJiXBhIapnnU+GuyaN5KJgptgkkBERKqoNNZaEoRkTT6W+WYjVDBClAVoBRmiLCBVm4dF8gbMrZuFHCkepnoJlcZaJglugn0SiIhIVcmafLzjm4VgGAHAsjS1+TYYRqz2zUKyJt9lMZIyJglERKQaHWqxzDcbgPIKkwB+2S5jmW82dLC+/wKpj0kCERGpZpwmF6GCscUEwUwjAKGCEamag84JjKzCJIGIiFQzRnsIotxGhvALURaQouUQdHfCJIGIiFQTihpL34O2aAUZoahROSKyBZMEIiJSTRUCbWpJqEKgyhGRLZgkEBGRanaICTa1JGwXOS2+O2GSQEREqtkmJaFKDoDURp4gyUCVHIAvpKHOCYyswiSBiIhUY4If5tbNAiC0mCg0bBcwt26WZeZFcg9MEtRUdw34biM0m6YhMmc2NJumAd9tbNhORHSTyJHiMbMuA9UIAABLHwXzbTUCMKMuAzlSvMtiJGWcllktRduAzbOBa1WAoEGALEE+XwAUfQ588SLwYDZwe6qroyQiUk2nAD/ofDQw1UvYKcUjyWRAquYgUrR5lrUbtouJ+EIa2mTthk4BbE1wF0wS1FC0Ddj4qOWuIEtNbnHtEvDRFOB3HwIDxrkiQiIi1fUI7YBd80besArkaMtvQQCe+uXHjKtAuhcmCY5Wd62hBQEA0FJPHRmA0LDf3GLA199JwREROVeP0A486Xsw9klwtMLNDZcYWkwQzOSG/Qq3qB4SERGRPZgkOFrR54Bg5dsqaICif6kbDxERkZ2YJDiYqfo8YO570BZZwrXLF9QNiIiIyE5MEhyovOoqdp8SbZqCdE9ZPcqrrqocGRERke28MkmoqqpCQkIC4uLiEBMTg9WrVzvldSuNtfiyPt6mKUi/qE+4oecvERGRe/DKJCEoKAh79+5FQUEBcnNz8frrr+PCBec063MKUiIi8hZemSRotVp07NgRAGAymSDLMmTZum/37cUpSImIyFu4ZZKwd+9ejB8/HhERERAEAZs3b262j8FgQO/eveHv74+kpCQcPHiwyeNVVVWIjY1Fz5498fzzzyMsLMxJ0XMKUiIi8g5umSQYjUbExsbCYDAoPv7xxx8jIyMDixYtwuHDhxEbG4uUlBScO3fOsk9oaCi+++47lJaW4sMPP0RFRYWzwgcAyxSkz9X+HjukBBwQB2KHlIDnan+PJJOBCQIREbk9t5xxMTU1FampLa9rkJWVhRkzZuDxxx8HAGRnZ2Pr1q1Ys2YN5s+f32Tf8PBwxMbGYt++fZg0aZJieSaTCSaTyXK/uroaACCKIkRRtDruG/c1wQ+bpeHYLA1v83m2vI67EkURkiR5RV1aw3p6F9bTu7Ce1x93BLdMElpTW1uL/Px8LFiwwLJNo9EgOTkZBw4cAABUVFSgY8eOCAoKwqVLl7B3717Mnj27pSKxZMkSLF68uNn2EydOIDAw0OrYyi6Y2t5J6XmnTkF35VzbO7o5SZJw8eJFlJSUQKNxy0Yqh2A9vQvr6V1YzwY1NTUOeR2PSxLOnz8PURQRHh7eZHt4eDiKiooAACdPnsTMmTMtHRafeeYZDBo0qMUyFyxYgIyMDMv96upqREZGol+/fggODrY6NlPHSwB+sq1CAKIiI9G/R4jNz3M3oiiipKQEt956K7RaravDUQ3r6V1YT+/CejYwt4i3l8clCdYYOnQoCgoKrN5fp9NBp9M1267Vam36kNn7gbT1ddyZRqPxqvq0hPX0Lqynd2E97T8fNXsNh5TiRGFhYdBqtc06IlZUVKBbt24uioqIiMj7eFxLgp+fH+Lj45GTk4OJEycCaLg2k5OTg/T09HaVbTAYYDAYvL7DCxGRNyivutpsxlpRFFF2wQRTx0vNvk13CvDjstU2csskoaamBiUlJZb7paWlKCgoQOfOnREVFYWMjAykpaUhISEBQ4cOxYoVK2A0Gi2jHeyl1+uh1+tRXV2NkBDb+wh0CvCDzkcDU72VCzwB0Plo0CmAEyoREdmivOoqRi/b3cr/2+b9w3Q+GuyaN5KJgg3cMkk4dOgQRo0aZblv7lSYlpaGdevWYfLkyfj555+RmZmJs2fPIi4uDl9++WWzzozO1iO0A3bNG6mc2Z46hajISGa2REQOUGmsbZIg6FCLcZpcjNEeQihqUIVA7BATsE1Kssxsa6qXUGms5f9cG7hlkjBy5Mg2p1FOT09v9+UFNfQI7dDsAyiKInRXzqF/jxCv70hDRORsyZp8LPPNRqhghCgL0AoyRFlAqjYPi+QNmFs3ixPY2cnjOi4SERGZJWvy8Y5vFoJhBADLKrzm22AYsdo3C8mafJfF6MnsThLq6upw6tQpFBcX4+LFi46MyWUMBgOio6ORmJjo6lCIiKgNOtRimW82ABkaQXmfhu0ylvlmQ4da5Z2oRTYlCZcvX8aqVatwzz33IDg4GL1798bAgQPRpUsX9OrVCzNmzEBeXp5asapOr9ejsLDQo+tARHSzGKfJRahgbDFBMNMIQKhgRKrmYOs7UjNWJwlZWVno3bs31q5di+TkZGzevBkFBQX4z3/+gwMHDmDRokWor6/HmDFjMHbsWBw/flzNuImI6CY3RnvIsrpuW0RZQIqWXwBtZXXHxby8POzduxd33HGH4uNDhw7FE088gezsbKxduxb79u1D//79HRYoERFRY6GosfQ9aItWkBEKx6xncDOxOkn46KOPrNpPp9Nh1qxZdgdERERkjSoEWkYztEWUBVQhEEFOiMubtHt0g3kRJW/AjotERJ5jh5hgU0vCdpH/221ld5Lw3nvvISYmBv7+/vD390dMTAzeffddR8bmdOy4SETkObZJSaiSAyC1kSdIMlAlB+ALaahzAvMidiUJmZmZePbZZzF+/Hhs2rQJmzZtwvjx4zFnzhxkZmY6OkYiIqJmTPDD3LpZAIQWE4WG7QLm1s2yzLxI1rNrxsVVq1Zh9erVmDJlimXbb37zGwwePBjPPPMMXnnlFYcFSEREdCPzWjk59fGYWZfRMOMims64qBVkVCPAMuMi18qxnV1JQl1dHRISEpptj4+PR319fbuDIiIiak3TtXKGo7z+aVwu3Yag0i9RX30WPsHdcLnPWFT3GYc5Pv6YA66VYw+7koTHHnsMq1atQlZWVpPt77zzDqZOneqQwIiIiFrTdK2cEKDX4xBHTMPx48fRv39/hHKtnHazOkkwr8QIAIIg4N1338WOHTtw1113AQByc3NRVlaGadOmOT5KJzEYDDAYDBBF0dWhEBERuZzVScK3337b5H58fMOKWidOnAAAhIWFISwsDEePHnVgeM6l1+uh1+tRXV2NkJAQV4dDRETkUlYnCW+++SbuuOMOLnVMRER0k7B6COSQIUMsqz327dsXFy5cUC0oIiIicj2rk4TQ0FD897//BQD8+OOPkCRJtaCIiIjI9ay+3PDwww/jnnvuQffu3SEIAhISElq89GBOJoiIiMhzWZ0kvPPOO3jooYdQUlKCP/zhD5gxYwaCgrhUBhERkbeyaZ6EsWPHAgDy8/Px7LPPMkkgIiLyYlb3SSgrK7P8vnbt2jYThPLycvujchGuAklERHSd1UlCYmIinn766VZXSLx06RJWr16NmJgYfPLJJw4J0Jm4CiQREdF1Vl9uKCwsxGuvvYb77rsP/v7+iI+PR0REBPz9/VFZWYnCwkIcPXoUd955J/76179i3LhxasZNREREKrO6JeGWW25BVlYWzpw5g7/97W/o378/zp8/j+PHjwMApk6divz8fBw4cIAJAhERkReweYGnDh06YNKkSZg0aZIa8RAREbmXumtA4Wag6HPgSiXQsRMw4AEgeiLg6+/q6FRl1yqQREREN4WibcDm2cC1KkDQALLUcHvsX8AXLwIPZgO3p7o6StVYfbmBiIjoZlFedRUnv94EeeOjkK9datgoS01u5WuXIH80BSe/3oTyqqsuilRdTBKIiIgaKa+6irHLdiBk+7OQZRkCZMX9BMiQZRkh25/F2GU7vDJRYJJARETUSKWxFsnSAYQKRmiE1vfVCECoYMS90jeoNNY6J0AnYpLQCCdTIiIiABijPQRRbiND+IUoC0jReuf8Og5JEiZNmoRXX30VmzdvxokTJwAAH374oSOKdipOpkRERAAQihpoBeXLDDfSCjJCUaNyRK7hkCThj3/8I3r16oX9+/fj6aefRlBQEBYvXuyIoomIiJyuCoE2tSRUIVDliFzDIUMghwwZgiFDhlju5+XlYf369Y4omoiIyOl2iAlItfISglaQsV1MxFMqx+QKDmlJOHfuXJP7iYmJ2L9/vyOKJiIicrptUhKq5ABIbVxxkGSgSg7AF9JQ5wTmZA5pSZg4cSJ+/vln9OnTB4MGDcK1a9cQExODuro6+Pr6OuIliIiInMYEP8ytm4XVvlmQZFlxlENDAiFgbt0smODn7BCdwq6WhFOnTjW5//XXX+P48ePIzs7GiBEjEB4ejqtXryI2NhaDBw92SKBERETOlCPFY2ZdBqoRAACWPgrm22oEYEZdBnKkeJfFqDa7WhIGDBiAuXPnYv78+ejYsaNle9++fdG3b19MnDjRsu3KlSvtDpKIiMgVdkrxSDIZkKo5iBRtHkJRgyoEYruYiC+koV7bgmBmV0vCV199he3bt6N///5Yt25dq/s2TiKIiIg8jQl+2CwNx+y6OZhStxCz6+ZgszTc6xMEwM4k4e6770Zubi6WLFmChQsXIj4+Hvv27XN0bERERORC7RrdMG3aNBQXF+P+++9HamoqJk2ahNLSUkfFRkRE5HSdAvyg87Ht9Kjz0aBTgPe1LDhkdMOYMWNQXV2NlStXYuvWrXjmmWeQmZmJwEDPmlzCYDDAYDBAFEVXh0JERC7SI7QDds0badNaDJ0C/NAjtIOKUbmGXUlCdnY28vLykJeXh2PHjkGj0SAmJgazZs1CbGwsNm7ciOjoaHz66adISEhwdMyq0ev10Ov1qK6uRkhIiKvDISIiF+kR2sErT/q2sitJeO2115CUlIRp06bhrrvuQnx8PDp0uP5mzpw5E6+//jqmT5+OI0eOOCxYIiIich67koQb50lQ8uSTT2LhwoX2FE9ERERuwOokYdiwYRgyZAji4uIwZMgQDBo0CP7+/i3u37VrV+zatcshQRIREZHzWZ0k3H///fj++++xfPlynDhxAoIgoH///oiLi2vy07VrVwCAIAi45557VAuciIiI1GV1kvDSSy9Zfj948CAmTpyImJgYCIKAdevWoaioCIIgIDw8HKdPn1YlWCIiInIeu/okzJ49GwaDAQ8++KBl27Zt2zBz5kykpaU5LDgiIiJyHbsmUzp27Bji4uKabBs3bhz+/ve/4+uvv3ZEXERERORidiUJiYmJWL9+fbPtgwYNwsGDB9sdFBEREbmeXZcbsrKyMHr0aJw8eRJz5sxBTEwMamtrsXz5coSFhTk6RiIiInIBu5KE+Ph45ObmIj09HXFxcfD19YUkSfDx8cF7773n6BiJiIjIBexeu2HAgAHYuXMnysrKUFBQAI1Gg/j4eHTv3t2R8REREZGL2JQkZGZmYsKECYiPj7dsi4qKQlRUlMMDIyIiIteyqePiTz/9hNTUVPTs2ROzZ8/GF198gdpa61fJIiIiIs9hU5KwZs0anD17Fh999BGCgoLw3HPPISwsDA8//DA2bNiAixcvqhUnEREROZnNQyA1Gg1GjBiBv/71ryguLkZubi6SkpLw9ttvIyIiAr/+9a+xbNkylJeXqxEvEREROYld8yTU1NRYfh84cCBeeOEF7N+/H6dOnUJaWhr27duHjz76yGFBOovBYEB0dDQSExNdHQoREZHL2ZUkhISE4JNPPmm2vUuXLnjyySexZcsWzJs3r93BOZter0dhYSHy8vJcHQoREZHL2ZUkyLKMt99+G7/61a8wfPhwPPfcczyxEhEReRm7kgQA+Pbbb3HnnXdi+PDhOHr0KEaMGOGRrQdERESkzO7JlD788EPcd999lvvff/89JkyYgB49emDOnDkOCY6IiIhcx66WhM6dOyMyMrLJtsGDB+Nvf/sbVq1a5ZDAiIiIyLXsShLi4uKwdu3aZttvvfVWlJWVtTsoIiIicj27Ljf86U9/wqhRo3D69Gn8/ve/x+DBg2E0GvH666+jT58+jo6RiIiIXMCuJOGuu+7CN998gz/84Q8YMWIEZFkGAPj7+2PTpk0ODZCIiIhcw+okYdiwYRgyZAji4uIQFxeHwYMHY8+ePTh37hzy8/MhSRKSkpIQFhamZrxERETkJFYnCffffz++//57LF++HCdOnIAgCOjfv78laYiLi4MkSWrGSkRERE5kdZLw0ksvWX4/ePAgJk6ciJiYGAiCgHXr1qGoqAiCICA8PBynT59WJVgiIiJyHrv6JMyePRsGgwEPPvigZdu2bdswc+ZMpKWlOSw4IiIich27hkAeO3YMcXFxTbaNGzcOf//73/H11187Ii4iIiJyMbuShMTERKxfv77Z9kGDBuHgwYPtDoqIiIhcz67LDVlZWRg9ejROnjyJOXPmICYmBrW1tVi+fDlHNxAREXkJu5KE+Ph45ObmIj09HXFxcfD19YUkSfDx8cF7773n6BiJiIjIBexe4GnAgAHYuXMnysrKUFBQAI1Gg/j4eHTv3t2R8REREZGL2J0kmEVFRSEqKsoRsRAREZEbsavjIhEREXk/JglERESkiEkCERERKfLKJOHUqVMYOXIkoqOjMXjwYK5MSUREZId2d1x0Rz4+PlixYgXi4uJw9uxZxMfHY9y4cQgICHB1aERERB7DK5OE7t27W4ZiduvWDWFhYbh48SKTBCIiIhu45eWGvXv3Yvz48YiIiIAgCNi8eXOzfQwGA3r37g1/f38kJSW1OB10fn4+RFFEZGSkylETERF5F7dsSTAajYiNjcUTTzyBhx56qNnjH3/8MTIyMpCdnY2kpCSsWLECKSkpKC4uRteuXS37Xbx4EdOmTcPq1atbfT2TyQSTyWS5X11dDQAQRRGiKLa7PqIoQpIkh5TlzlhP78J6ehfW07u0VU9H1V+QZVl2SEkqEQQBn332GSZOnGjZlpSUhMTERPztb38DAEiShMjISDzzzDOYP38+gIYT/3333YcZM2bgsccea/U1Xn75ZSxevLjZ9ry8PAQGBra7DpIk4eLFi+jcuTM0GrdsvHEI1tO7sJ7ehfX0Lm3Vs6amBomJibh06RKCg4Ptfh23bEloTW1tLfLz87FgwQLLNo1Gg+TkZBw4cAAAIMsypk+fjtGjR7eZIADAggULkJGRYblfXV2NyMhI9OvXr11vrpkoiigpKcGtt94KrVbb7vLcFevpXVhP78J6epe26mluEW8vj0sSzp8/D1EUER4e3mR7eHg4ioqKAAD79+/Hxx9/jMGDB1v6M7z//vsYNGiQYpk6nQ46na7Zdq1W67APmUajcWh57or19C6sp3dhPb1La/V0VN09LkmwxvDhwyFJkqvDICIi8mged8EmLCwMWq0WFRUVTbZXVFSgW7duLoqKiIjI+3hckuDn54f4+Hjk5ORYtkmShJycHAwbNqxdZRsMBkRHRyMxMbG9YRIREXk8t7zcUFNTg5KSEsv90tJSFBQUoHPnzoiKikJGRgbS0tKQkJCAoUOHYsWKFTAajXj88cfb9bp6vR56vR7V1dUICQlpbzWIiIg8mlsmCYcOHcKoUaMs980jD9LS0rBu3TpMnjwZP//8MzIzM3H27FnExcXhyy+/bNaZkYiIiOznlknCyJEj0db0Denp6UhPT3dSRERERDcfj+uToCb2SSAiIrqOSUIjer0ehYWFyMvLc3UoRERELsckgYiIiBQxSSAiIiJFTBKIiIhIEZMEIiIiUsQkoRGObiAiIrqOSUIjHN1ARER0HZMEIiIiUsQkgYiIiBQxSSAiIiJFTBKIiIhIEZOERji6gYiI6DomCY1wdAMREdF1TBKIiIhIEZMEIiIiUsQkgYiIiBQxSSAiIiJFTBKIiIhIEZOERjgEkoiI6DomCY1wCCQREdF1TBKIiIhIEZMEIiIiUsQkgYiIiBQxSSAiIiJFTBKIiIhIEZMEIiIiUsQkgYiIiBQxSWiEkykRERFdxyShEU6mREREdB2TBCIiIlLEJIGIiIgUMUkgIiIiRUwSiIiISBGTBCIiIlLEJIGIiIgUMUkgIiIiRUwSiIiISBGTBCIiIlLEJIGIiIgUMUlohGs3EBERXcckoRGu3UBERHQdkwQiIiJSxCSBiIiIFDFJICIiIkVMEoiIiEgRkwQiIiJSxCSBiIiIFDFJICIiIkVMEoiIiEgRkwQiIiJSxCSBiIiIFDFJICIiIkVMEoiIiEgRkwQiIiJSxCSBiIiIFDFJaMRgMCA6OhqJiYmuDoWIiMjlmCQ0otfrUVhYiLy8PFeHQkRE5HJMEoiIiEgRkwQiIiJSxCSBiIiIFDFJICIiIkVMEoiIiEgRkwQiIiJSxCSBiIiIFDFJICIiIkU+rg6AiIjoZlJedRWVxlqr9+8U4IceoR1UjKhlTBKIiIicpLzqKkYv2w1TvWT1c3Q+GuyaN9IliQIvNxARETlJpbHWpgQBAEz1kk0tD47EJIGIiIgUMUkgIiIiRUwSiIiISBGTBCIiIlLEJIGIiIgUcQgkERGRi+hQi3GaXIzRHkIoalCFQOwQE7BNSoIJfq4Oj0kCERGRKyRr8rHMNxuhghGiLEAryBBlAanaPCySN2Bu3SzkSPEujZGXG4iIiJwsWZOPd3yzEAwjAEAryE1ug2HEat8sJGvyXRYjwCSBiIjIqXSoxTLfbAAyNILyPg3bZSzzzYYOrplICfDiJOHBBx9Ep06dMGnSJFeHQkREZDFOk4tQwdhigmCmEYBQwYhUzUHnBKYUg8teWWXPPvssNmzY4OowiIiILDoF+GGsTz5EuY0M4ReiLCDV5xA6BbimE6PXJgkjR45EUFCQq8MgIiKy6BHaASMjtZa+B23RCjLuifJx2SqQbpkk7N27F+PHj0dERAQEQcDmzZub7WMwGNC7d2/4+/sjKSkJBw+6rjmGiIjIWrrgMECw8vQraOAfdIu6AbXCLZMEo9GI2NhYGAwGxcc//vhjZGRkYNGiRTh8+DBiY2ORkpKCc+fOOTlSIiIiGw14AJCtXAlSloAB49WNpxVuOU9CamoqUlNTW3w8KysLM2bMwOOPPw4AyM7OxtatW7FmzRrMnz/f5tczmUwwmUyW+9XV1QAAURQhiqLN5d1IFEVIkuSQstwZ6+ldWE/vwnq6kQHjofF/AbhWDQEtX3aQIQD+IZAGPADcUJ+26umo+rtlktCa2tpa5OfnY8GCBZZtGo0GycnJOHDggF1lLlmyBIsXL262/cSJEwgMDLQ7VjNJknDx4kWUlJRAo3HLxhuHYD29C+vpXVhP9xKYuBA99j0PGYJioiCjoWNjeeJLqCk91ezxtupZU1PjkDg9Lkk4f/48RFFEeHh4k+3h4eEoKiqy3E9OTsZ3330Ho9GInj17YtOmTRg2bJhimQsWLEBGRoblfnV1NSIjI9GvXz8EBwe3O2ZRFFFSUoJbb70VWq223eW5K9bTu7Ce3oX1dDP9+0Pq3h2af+qBa1WQBQ0EWbLcwj8E0oS/o/ttYxWf3lY9zS3i7eVxSYK1du7cafW+Op0OOp2u2XatVuuwD5lGo3Foee6K9fQurKd3YT3dTPQDQP9koHALhKJ/AVcrIXToBAwYDyF6ArS+/q0+vbV6OqruHpckhIWFQavVoqKiosn2iooKdOvWzUVRERER2cHXH4id3PDjhtz3gk0L/Pz8EB8fj5ycHMs2SZKQk5PT4uUEaxkMBkRHRyMxMbG9YRIREXk8t2xJqKmpQUlJieV+aWkpCgoK0LlzZ0RFRSEjIwNpaWlISEjA0KFDsWLFChiNRstoB3vp9Xro9XpUV1cjJCSkvdUgIiLyaG6ZJBw6dAijRo2y3Dd3KkxLS8O6deswefJk/Pzzz8jMzMTZs2cRFxeHL7/8sllnRiIiIrKfWyYJI0eOhCy3PmVleno60tPTnRQRERHRzcfj+iQQERGRc7hlS4KrGAwGGAwG1NfXA3DcOFNRFFFTU4Pq6mr3H5LTDqynd2E9vQvr6V3aqqf5/NVWq3xbBLm9JXihn376CZGRka4Og4iIqF1OnTqFnj172v18JgkKJEnC6dOnERQUBEGwbs3v1phncDx16pRDZnB0V6ynd2E9vQvr6V3aqqcsy7h8+TIiIiLaNT01Lzco0Gg07cq8WhIcHOzVH1oz1tO7sJ7ehfX0Lq3V0xFD+dlxkYiIiBQxSSAiIiJFTBKcQKfTYdGiRYqLSHkT1tO7sJ7ehfX0Ls6qJzsuEhERkSK2JBAREZEiJglERESkiEkCERERKWKSQERERIqYJDiIwWBA79694e/vj6SkJBw8eLDV/Tdt2oQBAwbA398fgwYNwrZt25wUqX2WLFmCxMREBAUFoWvXrpg4cSKKi4tbfc66desgCEKTH39/fydFbJ+XX365WcwDBgxo9TmediwBoHfv3s3qKQgC9Hq94v6eciz37t2L8ePHIyIiAoIgYPPmzU0el2UZmZmZ6N69Ozp06IDk5GQcP368zXJt/ftWW2v1rKurw4svvohBgwYhICAAERERmDZtGk6fPt1qmfZ89tXW1vGcPn16s5jHjh3bZrmedDwBKP6tCoKApUuXtlimo44nkwQH+Pjjj5GRkYFFixbh8OHDiI2NRUpKCs6dO6e4/9dff40pU6bgySefxLfffouJEydi4sSJOHLkiJMjt96ePXug1+vxzTff4KuvvkJdXR3GjBkDo9HY6vOCg4Nx5swZy8/JkyedFLH97rjjjiYx/9///V+L+3risQSAvLy8JnX86quvAAC//e1vW3yOJxxLo9GI2NhYGAwGxcf/+te/4q233kJ2djZyc3MREBCAlJQUXLt2rcUybf37dobW6nnlyhUcPnwYCxcuxOHDh/Hpp5+iuLgYv/nNb9os15bPvjO0dTwBYOzYsU1i/uijj1ot09OOJ4Am9Ttz5gzWrFkDQRDw8MMPt1quQ46nTO02dOhQWa/XW+6LoihHRETIS5YsUdz/kUceke+///4m25KSkuSnn35a1Tgd6dy5czIAec+ePS3us3btWjkkJMR5QTnAokWL5NjYWKv394ZjKcuy/Oyzz8r9+vWTJUlSfNwTjyUA+bPPPrPclyRJ7tatm7x06VLLtqqqKlmn08kfffRRi+XY+vftbDfWU8nBgwdlAPLJkydb3MfWz76zKdUzLS1NnjBhgk3leMPxnDBhgjx69OhW93HU8WRLQjvV1tYiPz8fycnJlm0ajQbJyck4cOCA4nMOHDjQZH8ASElJaXF/d3Tp0iUAQOfOnVvdr6amBr169UJkZCQmTJiAo0ePOiO8djl+/DgiIiLQt29fTJ06FWVlZS3u6w3Hsra2Fh988AGeeOKJVhc088Rj2VhpaSnOnj3b5HiFhIQgKSmpxeNlz9+3O7p06RIEQUBoaGir+9ny2XcXu3fvRteuXXH77bdj9uzZuHDhQov7esPxrKiowNatW/Hkk0+2ua8jjieThHY6f/48RFFEeHh4k+3h4eE4e/as4nPOnj1r0/7uRpIkPPfcc/jVr36FmJiYFve7/fbbsWbNGmzZsgUffPABJEnC3XffjZ9++smJ0domKSkJ69atw5dffolVq1ahtLQUI0aMwOXLlxX39/RjCQCbN29GVVUVpk+f3uI+nngsb2Q+JrYcL3v+vt3NtWvX8OKLL2LKlCmtLnhk62ffHYwdOxYbNmxATk4O/vKXv2DPnj1ITU2FKIqK+3vD8Vy/fj2CgoLw0EMPtbqfo44nV4Ekm+n1ehw5cqTN61vDhg3DsGHDLPfvvvtuDBw4EG+//TZeffVVtcO0S2pqquX3wYMHIykpCb169cI//vEPqzJ3T/Tee+8hNTUVERERLe7jiceSGjoxPvLII5BlGatWrWp1X0/87P/ud7+z/D5o0CAMHjwY/fr1w+7du3Hvvfe6MDL1rFmzBlOnTm2z47CjjidbEtopLCwMWq0WFRUVTbZXVFSgW7duis/p1q2bTfu7k/T0dHz++ef497//bfNy2r6+vhgyZAhKSkpUis7xQkNDcdttt7UYsycfSwA4efIkdu7ciaeeesqm53nisTQfE1uOlz1/3+7CnCCcPHkSX331lc3LJrf12XdHffv2RVhYWIsxe/LxBIB9+/ahuLjY5r9XwP7jySShnfz8/BAfH4+cnBzLNkmSkJOT0+SbV2PDhg1rsj8AfPXVVy3u7w5kWUZ6ejo+++wz7Nq1C3369LG5DFEU8cMPP6B79+4qRKiOmpoanDhxosWYPfFYNrZ27Vp07doV999/v03P88Rj2adPH3Tr1q3J8aqurkZubm6Lx8uev293YE4Qjh8/jp07d+KWW26xuYy2Pvvu6KeffsKFCxdajNlTj6fZe++9h/j4eMTGxtr8XLuPZ7u7PpK8ceNGWafTyevWrZMLCwvlmTNnyqGhofLZs2dlWZblxx57TJ4/f75l//3798s+Pj7ysmXL5GPHjsmLFi2SfX195R9++MFVVWjT7Nmz5ZCQEHn37t3ymTNnLD9Xrlyx7HNjPRcvXixv375dPnHihJyfny//7ne/k/39/eWjR4+6ogpWmTt3rrx79265tLRU3r9/v5ycnCyHhYXJ586dk2XZO46lmSiKclRUlPziiy82e8xTj+Xly5flb7/9Vv72229lAHJWVpb87bffWnr1//nPf5ZDQ0PlLVu2yN9//708YcIEuU+fPvLVq1ctZYwePVpeuXKl5X5bf9+u0Fo9a2tr5d/85jdyz5495YKCgiZ/ryaTyVLGjfVs67PvCq3V8/Lly/K8efPkAwcOyKWlpfLOnTvlO++8U+7fv7987do1SxmefjzNLl26JHfs2FFetWqVYhlqHU8mCQ6ycuVKOSoqSvbz85OHDh0qf/PNN5bH7rnnHjktLa3J/v/4xz/k2267Tfbz85PvuOMOeevWrU6O2DYAFH/Wrl1r2efGej733HOW9yQ8PFweN26cfPjwYecHb4PJkyfL3bt3l/38/OQePXrIkydPlktKSiyPe8OxNNu+fbsMQC4uLm72mKcey3//+9+Kn1NzXSRJkhcuXCiHh4fLOp1Ovvfee5vVv1evXvKiRYuabGvt79sVWqtnaWlpi3+v//73vy1l3FjPtj77rtBaPa9cuSKPGTNG7tKli+zr6yv36tVLnjFjRrOTvacfT7O3335b7tChg1xVVaVYhlrHk0tFExERkSL2SSAiIiJFTBKIiIhIEZMEIiIiUsQkgYiIiBQxSSAiIiJFTBKIiIhIEZMEIiIiUsQkgYiIiBQxSSCidhs5ciSee+65Fh/v3bs3VqxY4bR4iMgxuFQ0EakuLy8PAQEBrg6DiGzEJIGIVNelSxdXh0BEduDlBiJyiPr6eqSnpyMkJARhYWFYuHAhzEvD3Hi5QRAEvPvuu3jwwQfRsWNH9O/fH//85z8tj1dWVmLq1Kno0qULOnTogP79+2Pt2rXOrhLRTY9JAhE5xPr16+Hj44ODBw/izTffRFZWFt59990W91+8eDEeeeQRfP/99xg3bhymTp2KixcvAgAWLlyIwsJCfPHFFzh27BhWrVqFsLAwZ1WFiH7Byw1E5BCRkZF44403IAgCbr/9dvzwww944403MGPGDMX9p0+fjilTpgAAXn/9dbz11ls4ePAgxo4di7KyMgwZMgQJCQkAGloiiMj52JJARA5x1113QRAEy/1hw4bh+PHjEEVRcf/Bgwdbfg8ICEBwcDDOnTsHAJg9ezY2btyIuLg4vPDCC/j666/VDZ6IFDFJICKX8PX1bXJfEARIkgQASE1NxcmTJzFnzhycPn0a9957L+bNm+eKMIluakwSiMghcnNzm9z/5ptv0L9/f2i1WrvK69KlC9LS0vDBBx9gxYoVeOeddxwRJhHZgH0SiMghysrKkJGRgaeffhqHDx/GypUrsXz5crvKyszMRHx8PO644w6YTCZ8/vnnGDhwoIMjJqK2MEkgIoeYNm0arl69iqFDh0Kr1eLZZ5/FzJkz7SrLz88PCxYswI8//ogOHTpgxIgR2Lhxo4MjJqK2CLJ5IDMRERFRI+yTQERERIqYJBAREZEiJglERESkiEkCERERKWKSQERERIqYJBAREZEiJglERESkiEkCERERKWKSQERERIqYJBAREZEiJglERESk6P8DARA6T/7YuagAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "# Experimental central values as provided by HepData\n", + "data_central = np.array([\n", + " 1223.0, 3263.0, 4983.0, 6719.0, 8051.0, 8967.0, 9561.0, 9822.0, 9721.0, 9030.0, 7748.0, 6059.0, 4385.0, 2724.0, 1584.0, 749.0, 383.0, 11.0\n", + "])\n", + "\n", + "# Normalization for each bin. See Section below for more details.\n", + "bin_norm = np.array([0.125 for _ in range(predictions.size - 2)] + [0.250, 0.250])\n", + "\n", + "fig, ax = plt.subplots(figsize=(5.6, 3.9))\n", + "# Factor of `1e3` takes into account the unit conversion into `fb`\n", + "ax.plot(df_preds[\"bins\"], 1e3 * bin_norm * df_preds[\"predictions\"], 's', markersize=8, label=\"theory\")\n", + "ax.plot(df_preds[\"bins\"], data_central, 'o', markersize=8, label=\"data\")\n", + "ax.grid(True, alpha=0.5)\n", + "ax.set_yscale(\"log\")\n", + "ax.set_xlabel(\"bins\")\n", + "ax.set_ylabel(r\"$d\\sigma / dy_\\mu$ (fb)\")\n", + "ax.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "d4fbfdc3-a35b-44bf-bcc9-d6864f20bfca", + "metadata": {}, + "source": [ + "As we can see, the theory predictions agree fairly well with the experimental\n", + "measurements. However, in order to have a meaningful comparison, one needs to\n", + "include all sources of experimental and theory uncertainties. **Notice** that\n", + "in the plot above, we are accounting for the change in normalization related\n", + "to each bin, something that we can directly modify inside the PineAPPL grid as \n", + "we will see later." + ] + }, + { + "cell_type": "markdown", + "id": "ba2b1e4c-bc60-48dd-9b04-d349b4118707", + "metadata": {}, + "source": [ + "**NOTE:** If the two initial state hadrons are different (as is the case in\n", + "$pp$ collisions in which one of the protons is polarized), then one can convolve\n", + "the grid with **two** different PDF sets using the `convolve_with_two()` function:" + ] + }, + { + "cell_type": "markdown", + "id": "3c84ba6b-8c52-4e8c-8861-b6b8d0d39cef", + "metadata": {}, + "source": [ + "```python\n", + "# Convolve the two initial state hadrons with different PDF sets\n", + "predictions = grid.convolve_with_two(2212, polarized_pdf.xfxQ2, 2212, unpolarized_pdf.xfxQ2, unpolarized_pdf.alphasQ2)\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "3d56b36c-b888-4f2d-b3fb-2f8f689ff003", + "metadata": {}, + "source": [ + "**NOTE:** The same functions `convolve_with_one()` and `convolve_with_two()` also work for convolving FK tables with PDF sets." + ] + }, + { + "cell_type": "markdown", + "id": "c890454d-a432-45cd-9ac9-e7da9a2eca31", + "metadata": {}, + "source": [ + "## How can I get information on the contents of a given PineAPPL grid?\n", + "\n", + "Once the `PineAPPL` grid is loaded, one can look into its contents and modify them.\n", + "Below we illustrate how to extract the most common information from a given grid." + ] + }, + { + "cell_type": "markdown", + "id": "07cf5156-96ae-49d1-9e24-934e3f94ca40", + "metadata": {}, + "source": [ + "#### Contributing channels" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "b0d1a772-3be5-47df-99e2-96daa5ebf380", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0: [(2, -2, 1.0), (4, -4, 1.0)]\n", + "1: [(0, -4, 1.0), (0, -2, 1.0)]\n", + "2: [(22, -4, 1.0), (22, -2, 1.0)]\n", + "3: [(2, 0, 1.0), (4, 0, 1.0)]\n", + "4: [(2, 22, 1.0), (4, 22, 1.0)]\n", + "5: [(1, -1, 1.0), (3, -3, 1.0)]\n", + "6: [(0, -3, 1.0), (0, -1, 1.0)]\n", + "7: [(22, -3, 1.0), (22, -1, 1.0)]\n", + "8: [(1, 0, 1.0), (3, 0, 1.0)]\n", + "9: [(1, 22, 1.0), (3, 22, 1.0)]\n", + "10: [(5, -5, 1.0)]\n", + "11: [(0, -5, 1.0)]\n", + "12: [(22, -5, 1.0)]\n", + "13: [(5, 0, 1.0)]\n", + "14: [(5, 22, 1.0)]\n", + "15: [(22, 22, 1.0)]\n", + "16: [(-5, 22, 1.0), (-3, 22, 1.0), (-1, 22, 1.0)]\n", + "17: [(1, 22, 1.0), (3, 22, 1.0), (5, 22, 1.0)]\n" + ] + } + ], + "source": [ + "# To get all the contributing channels\n", + "for idx, c in enumerate(grid.channels()):\n", + " print(f\"{idx}: {c}\")" + ] + }, + { + "cell_type": "markdown", + "id": "967ef36d-865c-4a3f-b41c-151615c05a6a", + "metadata": {}, + "source": [ + "The above prints out all the contributing channels where the first\n", + "two elements of the tuple represent the particle PID (following the\n", + "PDG conventions). For instance, in the first entry, the up-antiup\n", + "`(2,-2)` and charm-anticharm `(4,-4)` initial states are merged\n", + "together, each with a multiplicative factor of `1`. In general this \n", + "number can be different from `1`, if the Monte Carlo decides to factor \n", + "out CKM values or electric charges." + ] + }, + { + "cell_type": "markdown", + "id": "61acaa44-0cba-4267-8611-6ff1c545ccb5", + "metadata": {}, + "source": [ + "#### Perturbative orders" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "6658a6c3-bb88-42fa-a567-9f5b6cebb1ac", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "shape: (7, 5)
indexasalflr
u32i64i64i64i64
00200
11200
21210
31201
40300
50310
60301
" + ], + "text/plain": [ + "shape: (7, 5)\n", + "┌───────┬─────┬─────┬─────┬─────┐\n", + "│ index ┆ as ┆ a ┆ lf ┆ lr │\n", + "│ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", + "│ u32 ┆ i64 ┆ i64 ┆ i64 ┆ i64 │\n", + "╞═══════╪═════╪═════╪═════╪═════╡\n", + "│ 0 ┆ 0 ┆ 2 ┆ 0 ┆ 0 │\n", + "│ 1 ┆ 1 ┆ 2 ┆ 0 ┆ 0 │\n", + "│ 2 ┆ 1 ┆ 2 ┆ 1 ┆ 0 │\n", + "│ 3 ┆ 1 ┆ 2 ┆ 0 ┆ 1 │\n", + "│ 4 ┆ 0 ┆ 3 ┆ 0 ┆ 0 │\n", + "│ 5 ┆ 0 ┆ 3 ┆ 1 ┆ 0 │\n", + "│ 6 ┆ 0 ┆ 3 ┆ 0 ┆ 1 │\n", + "└───────┴─────┴─────┴─────┴─────┘" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# To get the information on the orders:\n", + "orders = []\n", + "for idx, o in enumerate(grid.orders()):\n", + " orders.append(o.as_tuple())\n", + "\n", + "df_orders = pl.DataFrame(\n", + " np.array(orders),\n", + " schema=[\"as\", \"a\", \"lf\", \"lr\"]\n", + ")\n", + "df_orders.with_row_index()" + ] + }, + { + "cell_type": "markdown", + "id": "ef2ea5d4-e4b7-4760-aca1-fe7c2d8205f0", + "metadata": {}, + "source": [ + "The table above lists the perturbative orders contained in the\n", + "grid where the powers of the strong coupling $a_s$, the electroweak\n", + "coupling $a$, the factorization $\\ell_F = \\log(\\mu_F^2/Q^2)$ and renormalization $\\ell_R=\\log(\\mu_R^2/Q^2)$ \n", + "logs are shown. For instance, the first index shows that the grid \n", + "contains a leading-order (LO) which has the coupling $a_s^2$." + ] + }, + { + "cell_type": "markdown", + "id": "ac93839a-320c-4fd6-a758-6d0a677121ec", + "metadata": {}, + "source": [ + "#### Observable binning" + ] + }, + { + "cell_type": "markdown", + "id": "c7696983-778e-4d20-9c74-bb7e210cb530", + "metadata": {}, + "source": [ + "Grids often correspond to an dataset, which correspond to an observable measured at several kinematic configurations, or bins.\n", + "The bins might be just 1 dimensional, as is the case here with rapidity, or for more complicated measurements higher dimensional.\n", + "Each bin is characterized by its left and right border." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "4514380c-65e3-4116-a835-238212d68d10", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "shape: (18, 2)
dim 0 leftdim 0 right
f64f64
2.02.125
2.1252.25
2.252.375
2.3752.5
2.52.625
……
3.6253.75
3.753.875
3.8754.0
4.04.25
4.254.5
" + ], + "text/plain": [ + "shape: (18, 2)\n", + "┌────────────┬─────────────┐\n", + "│ dim 0 left ┆ dim 0 right │\n", + "│ --- ┆ --- │\n", + "│ f64 ┆ f64 │\n", + "╞════════════╪═════════════╡\n", + "│ 2.0 ┆ 2.125 │\n", + "│ 2.125 ┆ 2.25 │\n", + "│ 2.25 ┆ 2.375 │\n", + "│ 2.375 ┆ 2.5 │\n", + "│ 2.5 ┆ 2.625 │\n", + "│ … ┆ … │\n", + "│ 3.625 ┆ 3.75 │\n", + "│ 3.75 ┆ 3.875 │\n", + "│ 3.875 ┆ 4.0 │\n", + "│ 4.0 ┆ 4.25 │\n", + "│ 4.25 ┆ 4.5 │\n", + "└────────────┴─────────────┘" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# To get the bin configurations\n", + "bin_dims = grid.bin_dimensions()\n", + "\n", + "# Each element of bins is an object with a left and right limit and\n", + "# an associated bin normalization.\n", + "df = pl.DataFrame({})\n", + "for bin_dim in range(bin_dims):\n", + " df = pl.concat([df,pl.DataFrame({\n", + " f\"dim {bin_dim} left\": grid.bin_left(bin_dim),\n", + " f\"dim {bin_dim} right\": grid.bin_right(bin_dim),\n", + " })],how=\"vertical\",)\n", + "df" + ] + }, + { + "cell_type": "markdown", + "id": "1040c3e6-4e66-41e8-a45a-890c0f6749c7", + "metadata": {}, + "source": [ + "## How can I edit a grid?\n", + "\n", + "Sometimes it is required to modify the contents of a grid for various\n", + "reasons (like plotting for example, as we saw earlier). The contents \n", + "of PineAPPL grids can be edited in various ways. In our example, we \n", + "are going to change the normalization related to each bin in order for\n", + "our observable to match the experimental measurements (see Section above).\n", + "\n", + "Let us first check the normalizations before we modify them:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "16f768ea-4007-4ee9-be97-c862fbb8dfab", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "shape: (18, 2)
indexbin normalization
u32f64
00.125
10.125
20.125
30.125
40.125
……
130.125
140.125
150.125
160.25
170.25
" + ], + "text/plain": [ + "shape: (18, 2)\n", + "┌───────┬───────────────────┐\n", + "│ index ┆ bin normalization │\n", + "│ --- ┆ --- │\n", + "│ u32 ┆ f64 │\n", + "╞═══════╪═══════════════════╡\n", + "│ 0 ┆ 0.125 │\n", + "│ 1 ┆ 0.125 │\n", + "│ 2 ┆ 0.125 │\n", + "│ 3 ┆ 0.125 │\n", + "│ 4 ┆ 0.125 │\n", + "│ … ┆ … │\n", + "│ 13 ┆ 0.125 │\n", + "│ 14 ┆ 0.125 │\n", + "│ 15 ┆ 0.125 │\n", + "│ 16 ┆ 0.25 │\n", + "│ 17 ┆ 0.25 │\n", + "└───────┴───────────────────┘" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df_bins = pl.DataFrame({\"bin normalization\": grid.bin_normalizations()})\n", + "df_bins.with_row_index()" + ] + }, + { + "cell_type": "markdown", + "id": "4df9626d-918a-406a-bba5-754b47b2d398", + "metadata": {}, + "source": [ + "The column `bin normalization` shows the factor that all convolutions are\n", + "divided with. In our case, these values should be `1` in order to match\n", + "the experimental measurements." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "60745b45-0cb5-413d-8f84-cea84e557998", + "metadata": {}, + "outputs": [], + "source": [ + "# Extract the left & right bin limits\n", + "bin_limits = [\n", + " (left, right)\n", + " for left, right in zip(\n", + " grid.bin_left(bin_dims - 1), grid.bin_right(bin_dims - 1)\n", + " )\n", + "]\n", + "\n", + "# Multiply the normalization by a factor of `2`\n", + "normalizations = [1.0 for _ in grid.bin_normalizations()]\n", + "\n", + "# NOTE: `normalizations` have to be of type np.ndarray\n", + "remapper = pineappl.bin.BinRemapper(np.array(normalizations), bin_limits)\n", + "\n", + "# Modify the bin normalization\n", + "grid.set_remapper(remapper)\n", + "\n", + "# Save the modified grid\n", + "grid.write_lz4(\"./LHCB_DY_8TEV_custom_normalizations.pineappl.lz4\")" + ] + }, + { + "cell_type": "markdown", + "id": "79b623ae-3a2f-420b-b20c-7b6f8ad15e5c", + "metadata": {}, + "source": [ + "We can now check that the normalizations have indeed been changed:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "82efebec-966b-4a8f-89b1-c78077857752", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "shape: (18, 2)
indexbin normalization
u32f64
01.0
11.0
21.0
31.0
41.0
……
131.0
141.0
151.0
161.0
171.0
" + ], + "text/plain": [ + "shape: (18, 2)\n", + "┌───────┬───────────────────┐\n", + "│ index ┆ bin normalization │\n", + "│ --- ┆ --- │\n", + "│ u32 ┆ f64 │\n", + "╞═══════╪═══════════════════╡\n", + "│ 0 ┆ 1.0 │\n", + "│ 1 ┆ 1.0 │\n", + "│ 2 ┆ 1.0 │\n", + "│ 3 ┆ 1.0 │\n", + "│ 4 ┆ 1.0 │\n", + "│ … ┆ … │\n", + "│ 13 ┆ 1.0 │\n", + "│ 14 ┆ 1.0 │\n", + "│ 15 ┆ 1.0 │\n", + "│ 16 ┆ 1.0 │\n", + "│ 17 ┆ 1.0 │\n", + "└───────┴───────────────────┘" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Load our modified grids\n", + "grid_nrm = pineappl.grid.Grid.read(\"./LHCB_DY_8TEV_custom_normalizations.pineappl.lz4\")\n", + "df_nbins = pl.DataFrame({\"bin normalization\": grid_nrm.bin_normalizations()})\n", + "df_nbins.with_row_index()" + ] + }, + { + "cell_type": "markdown", + "id": "243c36ed-2adf-44e4-bd96-0b391815bc3c", + "metadata": {}, + "source": [ + "## Fast-Kernel (FK) tables" + ] + }, + { + "cell_type": "markdown", + "id": "61bbee09-dd78-4a9d-9830-8eb417545f72", + "metadata": {}, + "source": [ + "Fast-kernel (FK) tables are a special kind of grid, which contain the\n", + "Evolution Kernel Operators ([EKO](https://eko.readthedocs.io/))\n", + "which encode the DGLAP evolution equations. Furthermore, FK tables do not\n", + "have perturbative orders while PineAPPL grids do.\n", + "\n", + "You can load them as well, using the FK interface:" + ] + }, + { + "cell_type": "markdown", + "id": "866767df-261e-46c9-9038-acaa0ac7c214", + "metadata": {}, + "source": [ + "```python\n", + "# Loading a FK table\n", + "fk = pineappl.fk_table.FkTable.read(\"path/to/fktable.pineappl.lz4\")\n", + "```" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "PineAPPL", + "language": "python", + "name": "pineappl" + }, + "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.12.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pineappl_py/docs/source/recipes.rst b/pineappl_py/docs/source/recipes.rst deleted file mode 100644 index a599a888..00000000 --- a/pineappl_py/docs/source/recipes.rst +++ /dev/null @@ -1,74 +0,0 @@ -Recipes -======= - -Below we list some common use cases with their solutions. - - -How can I convolve a given PineAPPL grid with my PDF? ------------------------------------------------------- - -.. code:: python - - import pineappl - import lhapdf - g = pineappl.grid.Grid.read("path/to/grid.pineappl.lz4") - pdf = lhapdf.mkPDF("YourPDF", 0) - bins = g.convolve_with_one(2212, pdf.xfxQ2, pdf.alphasQ2) - -If the grid is actually an FkTable just replace - -.. code:: python - - g = pineappl.fk_table.FKTable.read("path/to/grid.pineappl.lz4") - -.. note:: - - For the :meth:`pineappl.grid.Grid.read` function, both ``.pineappl`` - and ``.pineappl.lz4`` extensions are acceptable, as long as they are - consistent (without ``.lz4`` the grid is assumed not to be compressed, with - it is assumed compressed). - -How can I edit a grid? ----------------------- - -.. code:: python - - import pineappl - g = pineappl.grid.Grid.read("path/to/grid.pineappl.lz4") - # edit the way you prefer - g = pineappl.grid.Grid.write_lz4("path/to/grid.pineappl.lz4") - -You can edit your grid in several ways. A few are listed in this section. - -Change normalizations -~~~~~~~~~~~~~~~~~~~~~ - -One possible option is to change the normalization related to each bin, or -change even the bins themselves. - -.. code:: python - - # each element of `bins` is an object with a left and right limit, and an - # associated normalization - limits = [(bin.left, bin.right) for bin in bins] - normalizations = [bin.norm for bin in bins] - remapper = pineappl.bin.BinRemapper(normalizations, limits) - g.set_remapper(remapper) - -For more details about :class:`pineappl.bin.BinRemapper` check also -the Rust documentation of :rustdoc:`pineappl::bin::BinRemapper `, e.g. -on how to treat multidimensional distributions. - -How can I get the bin configurations from a given PineAPPL grid? ----------------------------------------------------------------- - -.. code:: python - - import pineappl - g = pineappl.grid.Grid.read("path/to/grid.pineappl.lz4") - # first get the number of dimensions - bin_dims = g.bin_dimensions() - # now you can get each of them - for bin_dim in range(bin_dims): - bin_left = g.bin_left(bin_dim) - bin_right = g.bin_right(bin_dim) diff --git a/pineappl_py/docs/source/tutorials.rst b/pineappl_py/docs/source/tutorials.rst new file mode 100644 index 00000000..ee2c1941 --- /dev/null +++ b/pineappl_py/docs/source/tutorials.rst @@ -0,0 +1,17 @@ +Tutorials +========= + +The following tutorials provide a brief introduction to some of the +main features of the `PineAPPL` Python API. These tutorials showcase +common use cases of `PineAPPL`, from the convolution of a grid with +a PDF set to compute prediction observables to the editing of its +contents. A (slightly) more advanced tutorial is also covered in +which we implement a very simple and short Monte Carlo program and +produce leading-order (LO) predictions. + + +.. toctree:: + :hidden: + + A basic introduction + A more advanced tutorial diff --git a/pineappl_py/pyproject.toml b/pineappl_py/pyproject.toml index 2b0864f2..f8e19cb3 100644 --- a/pineappl_py/pyproject.toml +++ b/pineappl_py/pyproject.toml @@ -26,8 +26,11 @@ dependencies = ["numpy>=1.16.0,<2.0.0"] [project.optional-dependencies] cli = ["pineappl-cli"] docs = [ - "sphinx>=6.2.1", + "sphinx>=7.0.0", "sphinx_rtd_theme>=1.2.2", + "nbsphinx>=0.8.8", + "ipykernel>=6.13.0", + "polars>=1.8.0" ] test = ["pytest", "pytest-cov"]