{
"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": "e904212d-ee57-4664-8e78-fd0e8f2a665e",
"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": [],
"source": [
"# We first need to load the PDF set with LHAPDF\n",
"import lhapdf\n",
"import numpy as np\n",
"\n",
"# `Polars` is a better alternative to Pandas (written in Rust!)\n",
"# If you don't have it installed. Simply run `pip install polars`\n",
"import polars as pl\n",
"\n",
"lhapdf.setVerbosity(0)\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": [
"In order to convolve a grid, we need to specify the types of convolutions that are required.\n",
"This includes the polarization and the PDG IDs of the involved hadrons, as well as wether or\n",
"not the hadron is in the initial- or final-state.\n",
"\n",
"In our example, the grid involves two initial-state unpolarized protons. We can therefore\n",
"construct the convolution types as follows:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "a10346c8-cda8-4d71-9b10-518cab0ec38b",
"metadata": {},
"outputs": [],
"source": [
"from pineappl.convolutions import Conv, ConvType\n",
"\n",
"conv_type = ConvType(polarized=False, time_like=False)\n",
"conv_object = Conv(convolution_types=conv_type, pid=2212)"
]
},
{
"cell_type": "markdown",
"id": "813d3ba0-df32-4f2d-8364-73c625855bce",
"metadata": {},
"source": [
"Our grid can now be convolved with our PDF set using the `convolve()` function."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "d41e391c-affd-49ad-8d70-483bf952d4f3",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"shape: (18, 2)
bins
predictions
i64
f64
0
8.159231
1
23.890855
2
38.306994
3
51.278627
4
62.349204
…
…
13
22.640245
14
12.927123
15
6.138797
16
1.435175
17
0.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": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"predictions = grid.convolve(\n",
" pdg_convs=[conv_object, conv_object], # Similar convolutions for symmetric protons\n",
" xfxs=[pdf.xfxQ2, pdf.xfxQ2], # Similar PDF sets for symmetric protons\n",
" alphas=pdf.alphasQ2,\n",
")\n",
"df_preds = pl.DataFrame(\n",
" {\n",
" \"bins\": range(predictions.size),\n",
" \"predictions\": predictions,\n",
" }\n",
")\n",
"df_preds"
]
},
{
"cell_type": "markdown",
"id": "6ad4724b-e303-4200-a212-c24c81cb1d3f",
"metadata": {},
"source": [
"We can see that `convolve()` can perform convolutions with an arbitrary number of distributions. This is\n",
"why `pdf_convs` and `xfxs` are lists that respectively take all the types of convolutions and distributions\n",
"corresponding to the involved hadrons.\n",
"\n",
"**NOTE:** If the hadrons have the same type of convolutions and require the convolution to the same distribution,\n",
"then only one single element can be passed to the list:"
]
},
{
"cell_type": "markdown",
"id": "1cb9a279-cb32-4817-8d5a-0fde5213085b",
"metadata": {},
"source": [
"```python\n",
"# Pass the shared convolution type and distribution to all hadrons\n",
"predictions = grid.convolve(\n",
" pdg_convs=[conv_object],\n",
" xfxs=[pdf.xfxQ2],\n",
" alphas=pdf.alphasQ2,\n",
")\n",
"```"
]
},
{
"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": 7,
"id": "fee97bfd-e392-4a78-8850-a4093bbcb330",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAFrCAYAAABbtho0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABC4ElEQVR4nO3dfVhUZf4/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",
" [\n",
" 1223.0,\n",
" 3263.0,\n",
" 4983.0,\n",
" 6719.0,\n",
" 8051.0,\n",
" 8967.0,\n",
" 9561.0,\n",
" 9822.0,\n",
" 9721.0,\n",
" 9030.0,\n",
" 7748.0,\n",
" 6059.0,\n",
" 4385.0,\n",
" 2724.0,\n",
" 1584.0,\n",
" 749.0,\n",
" 383.0,\n",
" 11.0,\n",
" ]\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(\n",
" df_preds[\"bins\"],\n",
" 1e3 * bin_norm * df_preds[\"predictions\"],\n",
" \"s\",\n",
" markersize=8,\n",
" label=\"theory\",\n",
")\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:** As mentioned before, if the two initial state hadrons are different \n",
"(as is the case in $pp$ collisions in which one of the protons is polarized),\n",
"then one can convolve the grid with **two** different PDF sets using the `convolve()` \n",
"function:"
]
},
{
"cell_type": "markdown",
"id": "3c84ba6b-8c52-4e8c-8861-b6b8d0d39cef",
"metadata": {},
"source": [
"```python\n",
"# Define the convolution types for each of the (un)polarized hadrons\n",
"pol_type = ConvType(polarized=True, time_like=False) # `polarized = True`\n",
"pol_object = Conv(convolution_types=pol_type, pid=2212)\n",
"\n",
"unpol_type = ConvType(polarized=False, time_like=False)\n",
"unpol_object = Conv(convolution_types=conv_type, pid=2212)\n",
"\n",
"# Convolve the two initial state hadrons with different PDF sets\n",
"predictions = predictions = grid.convolve(\n",
" pdg_convs=[pol_object, conv_object],\n",
" xfxs=[polarized_pdf.xfxQ2, 2212, unpolarized_pdf.xfxQ2],\n",
" alphas=pdf.alphasQ2,\n",
")\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "3d56b36c-b888-4f2d-b3fb-2f8f689ff003",
"metadata": {},
"source": [
"**NOTE:** The same function `convolve` also works 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": 8,
"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": 9,
"id": "6658a6c3-bb88-42fa-a567-9f5b6cebb1ac",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"shape: (7, 6)
Index
as
a
lf
lr
la
u32
i64
i64
i64
i64
i64
0
0
2
0
0
0
1
1
2
0
0
0
2
1
2
1
0
0
3
1
2
0
1
0
4
0
3
0
0
0
5
0
3
1
0
0
6
0
3
0
1
0
"
],
"text/plain": [
"shape: (7, 6)\n",
"┌───────┬─────┬─────┬─────┬─────┬─────┐\n",
"│ Index ┆ as ┆ a ┆ lf ┆ lr ┆ la │\n",
"│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n",
"│ u32 ┆ i64 ┆ i64 ┆ i64 ┆ i64 ┆ i64 │\n",
"╞═══════╪═════╪═════╪═════╪═════╪═════╡\n",
"│ 0 ┆ 0 ┆ 2 ┆ 0 ┆ 0 ┆ 0 │\n",
"│ 1 ┆ 1 ┆ 2 ┆ 0 ┆ 0 ┆ 0 │\n",
"│ 2 ┆ 1 ┆ 2 ┆ 1 ┆ 0 ┆ 0 │\n",
"│ 3 ┆ 1 ┆ 2 ┆ 0 ┆ 1 ┆ 0 │\n",
"│ 4 ┆ 0 ┆ 3 ┆ 0 ┆ 0 ┆ 0 │\n",
"│ 5 ┆ 0 ┆ 3 ┆ 1 ┆ 0 ┆ 0 │\n",
"│ 6 ┆ 0 ┆ 3 ┆ 0 ┆ 1 ┆ 0 │\n",
"└───────┴─────┴─────┴─────┴─────┴─────┘"
]
},
"execution_count": 9,
"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(np.array(orders), schema=[\"as\", \"a\", \"lf\", \"lr\", \"la\"])\n",
"df_orders.with_row_index(\"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)$, renormalization $\\ell_R=\\log(\\mu_R^2/Q^2)$,\n",
"and fragmentation $\\ell_A=\\log(\\mu_A^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": 10,
"id": "4514380c-65e3-4116-a835-238212d68d10",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"shape: (18, 2)
dim 0 left
dim 0 right
f64
f64
2.0
2.125
2.125
2.25
2.25
2.375
2.375
2.5
2.5
2.625
…
…
3.625
3.75
3.75
3.875
3.875
4.0
4.0
4.25
4.25
4.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": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# To get the bin configurations\n",
"bin_dims = grid.bin_dimensions()\n",
"\n",
"# Get the bin specifications.\n",
"bin_specs = np.array(grid.bin_limits())\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(\n",
" [\n",
" df,\n",
" pl.DataFrame(\n",
" {\n",
" f\"dim {bin_dim} left\": bin_specs[:, bin_dim, 0],\n",
" f\"dim {bin_dim} right\": bin_specs[:, bin_dim, 1],\n",
" }\n",
" ),\n",
" ],\n",
" how=\"vertical\",\n",
" )\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": 11,
"id": "16f768ea-4007-4ee9-be97-c862fbb8dfab",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"shape: (18, 2)
index
bin normalization
u32
f64
0
0.125
1
0.125
2
0.125
3
0.125
4
0.125
…
…
13
0.125
14
0.125
15
0.125
16
0.25
17
0.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": 11,
"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": 12,
"id": "60745b45-0cb5-413d-8f84-cea84e557998",
"metadata": {},
"outputs": [],
"source": [
"# Import the Module to set bin specifications\n",
"from pineappl.boc import BinsWithFillLimits\n",
"\n",
"# Extract the left & right bin limits\n",
"bin_limits = [\n",
" [(bin_specs[b, d, 0], bin_specs[b, d, 1]) for d in range(bin_dims)]\n",
" for b in range(grid.len())\n",
"]\n",
"\n",
"# Set the normalization factor to `1`\n",
"normalizations = [1.0 for _ in grid.bin_normalizations()]\n",
"\n",
"# Instantiate the bin spec object\n",
"bin_configs = BinsWithFillLimits.from_limits_and_normalizations(\n",
" limits=bin_limits,\n",
" normalizations=normalizations,\n",
")\n",
"\n",
"# Set the updated bin configurations\n",
"grid.set_bwfl(bin_configs)\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": 13,
"id": "82efebec-966b-4a8f-89b1-c78077857752",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"