{ "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", "\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(\"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", "\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\n", " # 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\n", " # 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", " n_tuple = [\n", " q2,\n", " x1,\n", " x2,\n", " ] # Pass kinematics as list; order has to follow `[q2, x1, x2, ..., xn]`\n", " grid.fill(\n", " order=0,\n", " observable=np.abs(yll),\n", " channel=0,\n", " ntuple=n_tuple,\n", " weight=weight,\n", " )" ] }, { "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": [ "from pineappl.boc import (\n", " BinsWithFillLimits,\n", " Channel,\n", " Kinematics,\n", " ScaleFuncForm,\n", " Scales,\n", " Order,\n", ")\n", "from pineappl.convolutions import Conv, ConvType\n", "from pineappl.grid import Grid\n", "from pineappl.interpolation import (\n", " Interp,\n", " InterpolationMethod,\n", " MappingMethod,\n", " ReweightingMethod,\n", ")\n", "from pineappl.pids import PidBasis\n", "\n", "\n", "def grid_specs(\n", " orders: list[Order],\n", " channels: list[Channel],\n", " bins: np.ndarray,\n", ") -> Grid:\n", " \"\"\"Construct the PineAPPL grid based on various specifications. These include\n", " the types of kinematics involved, the types of convolutions required by the\n", " involved hadrons, and the interpolations required by each kinematic variables.\n", " \"\"\"\n", " ### Define the specs that define the Grid ###\n", " kinematics = [\n", " Kinematics.Scale(0), # Scale\n", " Kinematics.X(0), # momentum fraction x1\n", " Kinematics.X(1), # momentum fraction x2\n", " ]\n", " # Define the interpolation specs for each item of the Kinematics\n", " interpolations = [\n", " Interp(\n", " min=1e2,\n", " max=1e8,\n", " nodes=40,\n", " order=3,\n", " reweight_meth=ReweightingMethod.NoReweight,\n", " map=MappingMethod.ApplGridH0,\n", " interpolation_meth=InterpolationMethod.Lagrange,\n", " ), # Interpolation on the Scale\n", " Interp(\n", " min=2e-7,\n", " max=1,\n", " nodes=50,\n", " order=3,\n", " reweight_meth=ReweightingMethod.ApplGridX,\n", " map=MappingMethod.ApplGridF2,\n", " interpolation_meth=InterpolationMethod.Lagrange,\n", " ), # Interpolation on momentum fraction x1\n", " Interp(\n", " min=2e-7,\n", " max=1,\n", " nodes=50,\n", " order=3,\n", " reweight_meth=ReweightingMethod.ApplGridX,\n", " map=MappingMethod.ApplGridF2,\n", " interpolation_meth=InterpolationMethod.Lagrange,\n", " ), # Interpolation on momentum fraction x2\n", " ]\n", "\n", " # Construct the `Scales` object\n", " scale_funcs = Scales(\n", " ren=ScaleFuncForm.Scale(0),\n", " fac=ScaleFuncForm.Scale(0),\n", " frg=ScaleFuncForm.NoScale(0),\n", " )\n", "\n", " # Construct the type of convolutions and the convolution object\n", " # In our case we have symmetrical unpolarized protons in the initial state\n", " conv_type = ConvType(polarized=False, time_like=False)\n", " conv_object = Conv(convolution_types=conv_type, pid=2212)\n", " convolutions = [conv_object, conv_object]\n", "\n", " # Construct the Bin object\n", " bin_obj = BinsWithFillLimits.from_fill_limits(fill_limits=bins)\n", "\n", " return Grid(\n", " pid_basis=PidBasis.Evol,\n", " channels=channels,\n", " orders=orders,\n", " bins=bin_obj,\n", " convolutions=convolutions,\n", " interpolations=interpolations,\n", " kinematics=kinematics,\n", " scale_funcs=scale_funcs,\n", " )\n", "\n", "\n", "def generate_grid(calls: int) -> Grid:\n", " \"\"\"Generate the grid.\"\"\"\n", " # create a new luminosity function for the $\\gamma\\gamma$ initial state\n", " channels = [Channel([([22, 22], 1.0)])]\n", " # only LO $\\alpha_\\mathrm{s}^0 \\alpha^2 \\log^0(\\xi_\\mathrm{R})\n", " # \\log^0(\\xi_\\mathrm{F}) \\log^0(\\xi_\\mathrm{A})$$\n", " orders = [Order(0, 2, 0, 0, 0)]\n", " bins = np.arange(0, 2.4, 0.1)\n", "\n", " # Instantiate the PineAPPL Grid\n", " grid = grid_specs(orders, channels, bins)\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" ] } ], "source": [ "import lhapdf\n", "\n", "lhapdf.setVerbosity(0)\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(\n", " pdg_convs=[grid.convolutions[0]],\n", " xfxs=[pdf.xfxQ2],\n", " alphas=pdf.alphasQ2,\n", ")" ] }, { "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": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAFvCAYAAABAYhLAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsVUlEQVR4nO3df3RU9Z3/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": "pineko", "language": "python", "name": "pineko" }, "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.11.9" } }, "nbformat": 4, "nbformat_minor": 5 }