{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f40b7065-5763-4abc-91ff-40437cf86517",
   "metadata": {},
   "source": [
    "# Simple $e^+e^- \\to \\mu^+ \\mu^- $ Monte Carlo Event Generator example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "215df39b-69fc-4f74-b882-3c96f5526783",
   "metadata": {},
   "source": [
    "Differential cross section:\n",
    "$$\n",
    "\\frac{d \\sigma}{ d \\Omega} = \\frac{\\alpha^2}{4 s } ( 1+ \\cos^2(\\theta))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f7124c5f-9f86-4a0a-ac08-6e30cc55e01d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import random\n",
    "\n",
    "import hist\n",
    "import numpy as np\n",
    "\n",
    "import pylhe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f2553a00-1d72-4d4b-ba6d-c552ca0a7b59",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 1 / 127.4  # alpha QED\n",
    "aqcd = 0.1075  # alpha QCD\n",
    "EB = 209\n",
    "s = (2 * EB) ** 2\n",
    "theta_min = 0\n",
    "theta_max = math.pi\n",
    "phi_min = 0\n",
    "phi_max = 2 * math.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7242869e-7fbb-4eee-9d49-8a9d63724e35",
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://equation-database.readthedocs.io/en/latest/_autosummary/equation_database.isbn_9780471887416.html#equation_database.isbn_9780471887416.equation_6_32\n",
    "\n",
    "\n",
    "def dsigma(s, theta, _phi):\n",
    "    return (math.cos(theta) ** 2 + 1) / 4 * alpha**2 / s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "147dfe17-a5d9-4b68-a07f-6f8d8aa390a1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated integral: 1.47792134606482e-09\n",
      "Real integral:      1.477056777521761e-09\n"
     ]
    }
   ],
   "source": [
    "# Monte Carlo integration\n",
    "\n",
    "\n",
    "def monte_carlo_integration(\n",
    "    func, s, theta_min, theta_max, phi_min, phi_max, num_samples\n",
    "):\n",
    "    theta_samples = [random.uniform(theta_min, theta_max) for _ in range(num_samples)]\n",
    "    phi_samples = [random.uniform(phi_min, phi_max) for _ in range(num_samples)]\n",
    "    func_values = [\n",
    "        (phi_max - phi_min)\n",
    "        * (theta_max - theta_min)\n",
    "        * func(s, theta, phi)\n",
    "        * math.sin(theta)\n",
    "        for theta, phi in zip(theta_samples, phi_samples)\n",
    "    ]\n",
    "    maximum = np.max(func_values)\n",
    "    integral = np.mean(func_values)\n",
    "    return integral, maximum\n",
    "\n",
    "\n",
    "# Parameters\n",
    "\n",
    "num_samples = 1_000_000\n",
    "\n",
    "# Perform integration\n",
    "result, maximum = monte_carlo_integration(\n",
    "    dsigma, s, theta_min, theta_max, phi_min, phi_max, num_samples\n",
    ")\n",
    "print(f\"Estimated integral: {result}\")\n",
    "\n",
    "# https://equation-database.readthedocs.io/en/latest/_autosummary/equation_database.isbn_9780471887416.html#equation_database.isbn_9780471887416.equation_6_33\n",
    "sigma = 4 * math.pi / 3 * alpha**2 / s\n",
    "print(f\"Real integral:      {sigma}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "692b7fad-a519-49b8-9dbc-d7d8c27eaba9",
   "metadata": {},
   "source": [
    "## Weighted Events\n",
    "\n",
    "The information about the distribution is carried by the weights of the events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5c26e5e3-8e1e-470f-96a6-376a6e294bf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = 1_000_000\n",
    "integ = 0.0\n",
    "lheevents = []\n",
    "for _i in range(num_samples):\n",
    "    theta = random.uniform(theta_min, theta_max)\n",
    "    phi = random.uniform(phi_min, phi_max)\n",
    "    sig = (\n",
    "        (phi_max - phi_min)\n",
    "        * (theta_max - theta_min)\n",
    "        * dsigma(s, theta, phi)\n",
    "        * math.sin(theta)\n",
    "        / num_samples\n",
    "    )\n",
    "    integ += sig\n",
    "    # Fill the LHE event\n",
    "    e = pylhe.LHEEvent(\n",
    "        eventinfo=pylhe.LHEEventInfo(\n",
    "            nparticles=5,\n",
    "            pid=0,\n",
    "            weight=sig,  # The individual weight per event\n",
    "            scale=s,\n",
    "            aqed=alpha,\n",
    "            aqcd=aqcd,\n",
    "        ),\n",
    "        particles=[\n",
    "            pylhe.LHEParticle(\n",
    "                id=11,\n",
    "                status=-1,\n",
    "                mother1=0,\n",
    "                mother2=0,\n",
    "                color1=0,\n",
    "                color2=0,\n",
    "                px=0.0,\n",
    "                py=0.0,\n",
    "                pz=EB,\n",
    "                e=EB,\n",
    "                m=0.0,\n",
    "                lifetime=0,\n",
    "                spin=9.0,\n",
    "            ),\n",
    "            pylhe.LHEParticle(\n",
    "                id=-11,\n",
    "                status=-1,\n",
    "                mother1=0,\n",
    "                mother2=0,\n",
    "                color1=0,\n",
    "                color2=0,\n",
    "                px=0.0,\n",
    "                py=0.0,\n",
    "                pz=-EB,\n",
    "                e=EB,\n",
    "                m=0.0,\n",
    "                lifetime=0,\n",
    "                spin=9.0,\n",
    "            ),\n",
    "            pylhe.LHEParticle(\n",
    "                id=22,\n",
    "                status=2,\n",
    "                mother1=1,\n",
    "                mother2=2,\n",
    "                color1=0,\n",
    "                color2=0,\n",
    "                px=0.0,\n",
    "                py=0.0,\n",
    "                pz=EB - EB,\n",
    "                e=EB + EB,\n",
    "                m=0.0,\n",
    "                lifetime=0,\n",
    "                spin=9.0,\n",
    "            ),\n",
    "            pylhe.LHEParticle(\n",
    "                id=13,\n",
    "                status=1,\n",
    "                mother1=3,\n",
    "                mother2=3,\n",
    "                color1=0,\n",
    "                color2=0,\n",
    "                px=EB * math.sin(theta) * math.cos(phi),\n",
    "                py=EB * math.sin(theta) * math.sin(phi),\n",
    "                pz=EB * math.cos(theta),\n",
    "                e=EB,\n",
    "                m=0,\n",
    "                lifetime=0,\n",
    "                spin=9.0,\n",
    "            ),\n",
    "            pylhe.LHEParticle(\n",
    "                id=-13,\n",
    "                status=1,\n",
    "                mother1=3,\n",
    "                mother2=3,\n",
    "                color1=0,\n",
    "                color2=0,\n",
    "                px=-EB * math.sin(theta) * math.cos(phi),\n",
    "                py=-EB * math.sin(theta) * math.sin(phi),\n",
    "                pz=-EB * math.cos(theta),\n",
    "                e=EB,\n",
    "                m=0,\n",
    "                lifetime=0,\n",
    "                spin=9.0,\n",
    "            ),\n",
    "        ],\n",
    "    )\n",
    "    lheevents.append(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2fca49b0-2d93-476a-9b6b-835b72a24bb0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<event>\n",
      "  5      0  1.8831139316e-15  1.7472400000e+05  7.8492935636e-03  1.0750000000e-01\n",
      "   11  -1   0   0   0   0  0.00000000e+00  0.00000000e+00  2.09000000e+02  2.09000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "  -11  -1   0   0   0   0  0.00000000e+00  0.00000000e+00 -2.09000000e+02  2.09000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "   22   2   1   2   0   0  0.00000000e+00  0.00000000e+00  0.00000000e+00  4.18000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "   13   1   3   3   0   0  1.42473560e+02 -7.23202632e+01  1.34729597e+02  2.09000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "  -13   1   3   3   0   0 -1.42473560e+02  7.23202632e+01 -1.34729597e+02  2.09000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "</event>\n"
     ]
    }
   ],
   "source": [
    "# Fill the LHE init\n",
    "lheinit = pylhe.LHEInit(\n",
    "    initInfo=pylhe.LHEInitInfo(\n",
    "        beamA=11,\n",
    "        beamB=-11,\n",
    "        energyA=EB,\n",
    "        energyB=EB,\n",
    "        PDFgroupA=0,\n",
    "        PDFgroupB=0,\n",
    "        PDFsetA=0,\n",
    "        PDFsetB=0,\n",
    "        weightingStrategy=3,\n",
    "        numProcesses=1,\n",
    "    ),\n",
    "    procInfo=[\n",
    "        pylhe.LHEProcInfo(\n",
    "            xSection=integ,\n",
    "            error=0,\n",
    "            unitWeight=1,\n",
    "            procId=1,\n",
    "        )\n",
    "    ],\n",
    "    generators=[\n",
    "        pylhe.LHEGenerator(\n",
    "            name=\"pylhe\", version=pylhe.__version__, description=\"Monte Carlo Example\"\n",
    "        )\n",
    "    ],\n",
    ")\n",
    "lhef = pylhe.LHEFile(lheinit, lheevents)\n",
    "print(lhef.events[0].tolhe(rwgt=False, weights=False))\n",
    "lhef.tofile(\"weighted.lhe.gz\", rwgt=False, weights=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c94b9451-e7cc-4f40-9c1e-63dd8e180595",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html>\n",
       "<div style=\"display:flex; align-items:center;\">\n",
       "<div style=\"width:290px;\">\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"-10 -105 270 120\">\n",
       "<line x1=\"-5\" y1=\"0\" x2=\"255\" y2=\"0\" style=\"fill:none;stroke-width:2;stroke:currentColor\"/>\n",
       "<text text-anchor=\"middle\" x=\"0\" y=\"15\" style=\"fill:currentColor;\">\n",
       "0\n",
       "</text>\n",
       "<text text-anchor=\"middle\" x=\"250\" y=\"15\" style=\"fill:currentColor;\">\n",
       "3.14\n",
       "</text>\n",
       "<text text-anchor=\"middle\" x=\"125.0\" y=\"15\" style=\"fill:currentColor;\">\n",
       "Axis 0\n",
       "</text>\n",
       "<polyline points=\"  0,0   0,-14.4 12.5,-14.4 12.5,-41.5  25,-41.5  25,-64.9 37.5,-64.9 37.5,-82.6  50,-82.6  50,-94.2 62.5,-94.2 62.5,-99.3  75,-99.3  75,-99.1 87.5,-99.1 87.5,-98 100,-98 100,-94.5 112.5,-94.5 112.5,-93.3 125,-93.3 125,-92.8 137.5,-92.8 137.5,-93.5 150,-93.5 150,-96.8 162.5,-96.8 162.5,-100 175,-100 175,-98.8 187.5,-98.8 187.5,-93.4 200,-93.4 200,-83.1 212.5,-83.1 212.5,-64.6 225,-64.6 225,-41.4 237.5,-41.4 237.5,-14.3 250,-14.3 250,0\" style=\"fill:none; stroke:currentColor;\"/>\n",
       "</svg>\n",
       "</div>\n",
       "<div style=\"flex=grow:1;\">\n",
       "Regular(20, 0, 3.14159, label='Axis 0')<br/>\n",
       "<hr style=\"margin-top:.2em; margin-bottom:.2em;\"/>\n",
       "Weight() Σ=WeightedSum(value=1.47739e-09, variance=2.46116e-24)\n",
       "\n",
       "</div>\n",
       "</div>\n",
       "</html>"
      ],
      "text/plain": [
       "Hist(Regular(20, 0, 3.14159, label='Axis 0'), storage=Weight()) # Sum: WeightedSum(value=1.47739e-09, variance=2.46116e-24)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr = pylhe.to_awkward(\n",
    "    pylhe.LHEFile.fromfile(\"weighted.lhe.gz\", with_attributes=True).events\n",
    ")\n",
    "hist1 = hist.Hist.new.Reg(20, 0, math.pi).Weight()\n",
    "hist1.fill(arr.particles.vector[:, -1].theta, weight=arr.eventinfo.weight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b43498cd-1fe7-49c1-a24c-5b66cf4e257f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7fec93658c20>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAHACAYAAAD+yCF8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHnpJREFUeJzt3X9wVeWd+PHPJSQBKkS0RaCmoK1CKdVSY1cUbXehIP5Y67Rb21FqbTtTFfw5uyu0rB2pbWTGWvpjFwdL0a5VaItYZ1sraPlRR+0WBXXLLm6tLOwqw1htouAkJpz9Y8d8zZckcC/PSXJzX6+Z+0duzj3nuYeH3Pd97k1uIcuyLAAAEhjU1wMAAAYOYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACTTZ2GxadOmOP/882Ps2LFRKBTi/vvv7/Pj3XfffTFr1qx45zvfGYVCIbZu3ZrrmABgoOmzsNi7d2+cfPLJ8f3vf7/fHG/v3r1xxhlnxC233NIrYwKAgWZwXx149uzZMXv27G6/39raGgsXLowf//jH8ec//zkmT54cixcvjo997GO5HC8iYs6cORERsWPHjpKOAQCVrs/C4mAuu+yy2LFjR6xcuTLGjh0ba9asibPPPjueffbZOOGEE/p6eABAF/rlmzeff/75uPfee+OnP/1pnHnmmfHe9743/vZv/zamTZsWK1as6OvhAQDd6Jdh8dRTT0WWZXHiiSfGEUcc0XHZuHFjPP/88xHxfy9XFAqFHi/z5s3r43sCAJWlX74Usn///qiqqoonn3wyqqqqOn3viCOOiIiId7/73fHv//7vPe5n5MiRuY0RADhQvwyLKVOmRHt7e+zZsyfOPPPMLreprq6OiRMn9vLIAICe9FlYvP766/GHP/yh4+sXXnghtm7dGkcddVSceOKJcfHFF8fnPve5+Na3vhVTpkyJl19+OX7961/HBz/4wTjnnHOSHu8973lPRES88sorsXPnznjxxRcjImL79u0RETF69OgYPXr04dxdAKgMWR9Zv359FhEHXC699NIsy7KstbU1u/HGG7Px48dn1dXV2ejRo7MLL7wwe+aZZ3I5XpZl2YoVK7rc5mtf+9rh32EAqACFLMuyPugZAGAA6pe/FQIAlCdhAQAk0+tv3ty/f3+8+OKLMXz48CgUCr19eACgBFmWxWuvvRZjx46NQYO6X5fo9bB48cUXo76+vrcPCwAksGvXrjj22GO7/X6vh8Xw4cMj4v8GNmLEiN4+PABQgubm5qivr+94HO9Or4fFWy9/jBgxQlgAQJk52NsYvHkTAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJBMr39sOsBAl2VZvPFme+7HGVpdddCPsIbeJiwAEnvjzfaYdONDuR9n26JZMazGj3H6Fy+FAADJSF36hWKWjve1tkXDzY9ERMTmhdP7xTM2S9J0Z/PCGTGspqrHbYqZ0/ta26Ph5oeTjhFS6vufyBClLx2/9cO4r1mSpjvDaqoOOjeG1QyOHbec20sjgnz5SQhQpva19rzKV+rqnhU4DoewoN85lKXjvFiSppwUM/+KWd2zAsfhMHPodw5l6Ti/Y5e2JH2wZ46Hw7NHoJwIC0ggz5ULzx55u6HVVbFt0azk+7UCRyp+WgGUkUKhIDTp18xOKFExzxyLfROdZ49AuRIWUKJinjn6dUKgUggLoGLl9YfZ8nwzL/R3wgKoWOX+h9mgP/JZIQBAMlYsACK/P8w2tLpv/tgb9BVhARB9+4fZYCDxUggAkIywAACSERYAQDJeUIR+zkdjA+VEWHDIivljQhH+oFAqPhobKCd+qnDISv1jQhH+oBBApRAW0A/5aOz/p9iVsmJYKetanufFS3ADn7CgJHn9MaEIf1Aoonc+Gjuv925EpH3wOJyVMkqTZ3x6CW7g869LSfwxofKX13s3Ijx4QCXzPx9ILuVqyNv3ZaUsP8W+/Fbsv2G5vQRH6YQFVJC83rsR0fnBI6/VECtl+Sn25bdhNYNjxy3n5jgiypX/oVBBeuO9G0Bl8xMGSCLP1ZC3HwPo34QFkITVECDCZ4UAAAkJCwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkExRYdHW1hYLFy6M4447LoYOHRrHH398LFq0KPbv35/X+ACAMjK4mI0XL14ct99+e9x1113xgQ98IDZv3hyXXXZZ1NXVxTXXXJPXGAGAMlFUWDz++ONxwQUXxLnnnhsREePHj4977703Nm/enMvgAIDyUtRLIdOmTYtHHnkknnvuuYiIePrpp+PRRx+Nc845J5fBAQDlpagVixtuuCGamppi4sSJUVVVFe3t7fGNb3wjPvvZz3Z7m5aWlmhpaen4urm5ufTRAgD9WlErFqtWrYq777477rnnnnjqqafirrvuiltvvTXuuuuubm/T2NgYdXV1HZf6+vrDHjQA0D8VFRZ/93d/F/Pnz4/PfOYz8cEPfjDmzJkT1113XTQ2NnZ7mwULFkRTU1PHZdeuXYc9aACgfyrqpZB9+/bFoEGdW6SqqqrHXzetra2N2tra0kYHwICyr7U9l/0Ora6KQqGQy74pTlFhcf7558c3vvGNeM973hMf+MAHYsuWLXHbbbfFF77whbzGB8AA0nDzw7nsd9uiWTGspqiHNHJS1L/C9773vfiHf/iHuPLKK2PPnj0xduzY+PKXvxw33nhjXuOjBFmWxRtvpn9WkNczDQAGjkKWZVlvHrC5uTnq6uqiqakpRowY0ZuHrhj7Wtti0o0P5XoMzw6AQ1XMk519rW3RcPMjERGxeeH0Hn/O7Gtt71gB8TMpf4f6+O1fAYBcFQqFQ37QH1YzOHbccm7OIyJPwmKA27xwRgyrqer2+8U8O3i7odXd7xOAyiUsBrhhNVU9xoJnBwCk5GPTAYBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZAb39QAqVZZl8cab7bnse19rPvsFgIMRFn3kjTfbY9KND/X1MAAGhDyfUA2tropCoZDb/gcaYQFA2Wu4+eHc9r1t0awYVuPh8lAVfab+53/+J2644YZ48MEH44033ogTTzwxli9fHqecckoe46sImxfOiGE1Vbnse2h1PvsFgK4UFRavvvpqnHHGGfGXf/mX8eCDD8aoUaPi+eefjyOPPDKn4VWGYTVVahigSEOrq2Lbolm57Htfa3uuqyADWVGPZosXL476+vpYsWJFx3Xjx49PPSYAOKhCoeBJWT9U1K+bPvDAA9HQ0BB/8zd/E6NGjYopU6bEHXfckdfYAIAyU1RY/PGPf4ylS5fGCSecEA899FBcfvnlcfXVV8ePfvSjbm/T0tISzc3NnS4AwMBU1BrS/v37o6GhIb75zW9GRMSUKVPi97//fSxdujQ+97nPdXmbxsbGuOmmmw5/pABAv1fUisWYMWNi0qRJna57//vfHzt37uz2NgsWLIimpqaOy65du0obKQDQ7xW1YnHGGWfE9u3bO1333HPPxbhx47q9TW1tbdTW1pY2OgCgrBS1YnHdddfFE088Ed/85jfjD3/4Q9xzzz2xbNmymDt3bl7jAwDKSFFhceqpp8aaNWvi3nvvjcmTJ8fXv/71WLJkSVx88cV5jQ8AKCNF/wLweeedF+edd14eYwEAypyPTQcAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkDissGhsbo1AoxLXXXptoOABAOSs5LH73u9/FsmXL4qSTTko5HgCgjJUUFq+//npcfPHFcccdd8TIkSNTjwkAKFMlhcXcuXPj3HPPjRkzZqQeDwBQxgYXe4OVK1fGU089Fb/73e8OafuWlpZoaWnp+Lq5ubnYQwIAZaKosNi1a1dcc801sXbt2hgyZMgh3aaxsTFuuummkgbXH2RZFm+82Z58v/ta0+8TAPpaIcuy7FA3vv/+++PCCy+Mqqqqjuva29ujUCjEoEGDoqWlpdP3Irpesaivr4+mpqYYMWJEgruQr32tbTHpxodyPca2RbNiWE3Ri0cA5OTtP/s3L5wRw2qqDnKL0gytropCoZDLvlNrbm6Ourq6gz5+F/VoNn369Hj22Wc7XXfZZZfFxIkT44YbbjggKiIiamtro7a2tpjDAEC/0XDzw7nteyA+sSzq3gwfPjwmT57c6bp3vOMdcfTRRx9w/UCUV7UOrc6nhAGgtw2sTMrZsJqqAVeWABxoaHVVbFs0K5d972ttz3UVpK8d9qPkhg0bEgwDAPqPQqHgiWSJfFYIAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJFNUWDQ2Nsapp54aw4cPj1GjRsUnPvGJ2L59e15jAwDKTFFhsXHjxpg7d2488cQTsW7dumhra4uZM2fG3r178xofAFBGBhez8a9+9atOX69YsSJGjRoVTz75ZJx11llJBwYAlJ/Deo9FU1NTREQcddRRSQYDAJS3olYs3i7Lsrj++utj2rRpMXny5G63a2lpiZaWlo6vm5ubSz0kANDPlbxiMW/evHjmmWfi3nvv7XG7xsbGqKur67jU19eXekgAoJ8rKSyuuuqqeOCBB2L9+vVx7LHH9rjtggULoqmpqeOya9eukgYKAPR/Rb0UkmVZXHXVVbFmzZrYsGFDHHfccQe9TW1tbdTW1pY8QACgfBQVFnPnzo177rknfv7zn8fw4cNj9+7dERFRV1cXQ4cOzWWAAED5KOqlkKVLl0ZTU1N87GMfizFjxnRcVq1aldf4AIAyUvRLIQAA3fFZIQBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJDM4L4eQApZlsUbb7bnsu99rfnsFwAGogERFm+82R6Tbnyor4cBAEXJ68nr0OqqKBQKuez7YAZEWABAOWq4+eFc9rtt0awYVtM3D/EDLiw2L5wRw2qqctn30Op89gsAA8WAC4thNVV9VmkAcDBDq6ti26JZyfe7r7U9txWQYngEBoBeVCgUBvQTYL9uCgAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZEoKi3/6p3+K4447LoYMGRKnnHJK/OY3v0k9LgCgDBUdFqtWrYprr702vvrVr8aWLVvizDPPjNmzZ8fOnTvzGB8AUEaKDovbbrstvvjFL8aXvvSleP/73x9LliyJ+vr6WLp0aR7jAwDKSFFh0draGk8++WTMnDmz0/UzZ86Mxx57LOnAAIDyM7iYjV9++eVob2+PY445ptP1xxxzTOzevbvL27S0tERLS0vH101NTRER0dzcXOxYu7WvtS32t+zr2G9bTVF3CwDKXt6PhW89bmdZ1uN2JR21UCh0+jrLsgOue0tjY2PcdNNNB1xfX19fyqEPasySXHYLAGUjz8fC1157Lerq6rr9flFh8c53vjOqqqoOWJ3Ys2fPAasYb1mwYEFcf/31HV/v378/XnnllTj66KO7jZFSNDc3R319fezatStGjBiRbL8DhfPTM+ene85Nz5yfnjk/PSun85NlWbz22msxduzYHrcrKixqamrilFNOiXXr1sWFF17Ycf26deviggsu6PI2tbW1UVtb2+m6I488spjDFmXEiBH9/h+nLzk/PXN+uufc9Mz56Znz07NyOT89rVS8peiXQq6//vqYM2dONDQ0xNSpU2PZsmWxc+fOuPzyy0saJAAwcBQdFhdddFH86U9/ikWLFsVLL70UkydPjl/+8pcxbty4PMYHAJSRkt68eeWVV8aVV16ZeiyHpba2Nr72ta8d8LIL/8f56Znz0z3npmfOT8+cn54NxPNTyA72eyMAAIfIh5ABAMkICwAgGWEBACRTVmFR7Me1b9y4MU455ZQYMmRIHH/88XH77bf30kj7RjHnZ8OGDVEoFA64/Md//Ecvjrh3bNq0Kc4///wYO3ZsFAqFuP/++w96m0qaO8Wen0qaO42NjXHqqafG8OHDY9SoUfGJT3witm/fftDbVcr8KeX8VNL8Wbp0aZx00kkdf6Ni6tSp8eCDD/Z4m4Ewd8omLIr9uPYXXnghzjnnnDjzzDNjy5Yt8ZWvfCWuvvrqWL16dS+PvHeU+nH227dvj5deeqnjcsIJJ/TSiHvP3r174+STT47vf//7h7R9pc2dYs/PWyph7mzcuDHmzp0bTzzxRKxbty7a2tpi5syZsXfv3m5vU0nzp5Tz85ZKmD/HHnts3HLLLbF58+bYvHlz/NVf/VVccMEF8fvf/77L7QfM3MnKxEc+8pHs8ssv73TdxIkTs/nz53e5/d///d9nEydO7HTdl7/85ey0007LbYx9qdjzs379+iwisldffbUXRtd/RES2Zs2aHreptLnzdodyfip17mRZlu3ZsyeLiGzjxo3dblPJ8+dQzk8lz58sy7KRI0dmP/jBD7r83kCZO2WxYlHKx7U//vjjB2w/a9as2Lx5c7z55pu5jbUvHM7H2U+ZMiXGjBkT06dPj/Xr1+c5zLJRSXPncFTi3Hnr05mPOuqobrep5PlzKOfnLZU2f9rb22PlypWxd+/emDp1apfbDJS5UxZhUcrHte/evbvL7dva2uLll1/Obax9oZTzM2bMmFi2bFmsXr067rvvvpgwYUJMnz49Nm3a1BtD7tcqae6UolLnTpZlcf3118e0adNi8uTJ3W5XqfPnUM9Ppc2fZ599No444oiora2Nyy+/PNasWROTJk3qctuBMnfSflh7zor5uPbutu/q+oGimPMzYcKEmDBhQsfXU6dOjV27dsWtt94aZ511Vq7jLAeVNneKUalzZ968efHMM8/Eo48+etBtK3H+HOr5qbT5M2HChNi6dWv8+c9/jtWrV8ell14aGzdu7DYuBsLcKYsVi1I+rn306NFdbj948OA4+uijcxtrXyjl/HTltNNOi//8z/9MPbyyU0lzJ5WBPneuuuqqeOCBB2L9+vVx7LHH9rhtJc6fYs5PVwby/KmpqYn3ve990dDQEI2NjXHyySfHd77znS63HShzpyzC4u0f1/5269ati9NPP73L20ydOvWA7deuXRsNDQ1RXV2d21j7QinnpytbtmyJMWPGpB5e2amkuZPKQJ07WZbFvHnz4r777otf//rXcdxxxx30NpU0f0o5P10ZqPOnK1mWRUtLS5ffGzBzp4/eNFq0lStXZtXV1dny5cuzbdu2Zddee232jne8I9uxY0eWZVk2f/78bM6cOR3b//GPf8yGDRuWXXfdddm2bduy5cuXZ9XV1dnPfvazvroLuSr2/Hz729/O1qxZkz333HPZv/3bv2Xz58/PIiJbvXp1X92F3Lz22mvZli1bsi1btmQRkd12223Zli1bsv/6r//KsszcKfb8VNLcueKKK7K6urpsw4YN2UsvvdRx2bdvX8c2lTx/Sjk/lTR/FixYkG3atCl74YUXsmeeeSb7yle+kg0aNChbu3ZtlmUDd+6UTVhkWZb94z/+YzZu3LispqYm+/CHP9zpV5ouvfTS7KMf/Win7Tds2JBNmTIlq6mpycaPH58tXbq0l0fcu4o5P4sXL87e+973ZkOGDMlGjhyZTZs2LfvFL37RB6PO31u/3vb/Xy699NIsy8ydYs9PJc2drs5LRGQrVqzo2KaS508p56eS5s8XvvCFjp/J73rXu7Lp06d3REWWDdy549NNAYBkyuI9FgBAeRAWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAkigUCnH//ff39TCAPiYsgIiIeOyxx6KqqirOPvvskm7/0ksvxezZs0s+/urVq2PSpElRW1sbkyZNijVr1pS8L6DvCAsgIiJ++MMfxlVXXRWPPvpo7Ny5s+jbjx49Ompra0s69uOPPx4XXXRRzJkzJ55++umYM2dOfPrTn47f/va3Je0P6DvCAoi9e/fGT37yk7jiiivivPPOizvvvLPT9xctWhRjx46NP/3pTx3X/fVf/3WcddZZsX///ojo/FJIa2trzJs3L8aMGRNDhgyJ8ePHR2NjY7fHX7JkSXz84x+PBQsWxMSJE2PBggUxffr0WLJkSeq7CuRMWACxatWqmDBhQkyYMCEuueSSWLFiRbz98wm/+tWvxvjx4+NLX/pSRETcfvvtsWnTpvjnf/7nGDTowB8j3/3ud+OBBx6In/zkJ7F9+/a4++67Y/z48d0e//HHH4+ZM2d2um7WrFnx2GOPpbmDQK8Z3NcDAPre8uXL45JLLomIiLPPPjtef/31eOSRR2LGjBkREVFVVRV33313fOhDH4r58+fH9773vVi2bFmMGzeuy/3t3LkzTjjhhJg2bVoUCoVut3vL7t2745hjjul03THHHBO7d+9OcO+A3mTFAirc9u3b41//9V/jM5/5TEREDB48OC666KL44Q9/2Gm7448/Pm699dZYvHhxnH/++XHxxRd3u8/Pf/7zsXXr1pgwYUJcffXVsXbt2oOOo1AodPo6y7IDrgP6PysWUOGWL18ebW1t8e53v7vjuizLorq6Ol599dUYOXJkx/WbNm2Kqqqq2LFjR7S1tcXgwV3/CPnwhz8cL7zwQjz44IPx8MMPx6c//emYMWNG/OxnP+ty+9GjRx+wOrFnz54DVjGA/s+KBVSwtra2+NGPfhTf+ta3YuvWrR2Xp59+OsaNGxc//vGPO7ZdtWpV3HfffbFhw4bYtWtXfP3rX+9x3yNGjIiLLroo7rjjjli1alWsXr06XnnllS63nTp1aqxbt67TdWvXro3TTz/98O8k0KusWEAF+5d/+Zd49dVX44tf/GLU1dV1+t6nPvWpWL58ecybNy/++7//O6644opYvHhxTJs2Le68884499xzY/bs2XHaaacdsN9vf/vbMWbMmPjQhz4UgwYNip/+9KcxevToOPLII7scxzXXXBNnnXVWLF68OC644IL4+c9/Hg8//HA8+uijedxtIEdWLKCCLV++PGbMmHFAVEREfPKTn4ytW7fGk08+GZ///OfjIx/5SMybNy8iIj7+8Y/HvHnz4pJLLonXX3/9gNseccQRsXjx4mhoaIhTTz01duzYEb/85S+7/A2SiIjTTz89Vq5cGStWrIiTTjop7rzzzli1alX8xV/8Rdo7DOSukL39d8oAAA6DFQsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkMz/AtupipSSbcJuAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist1.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77d19d7a-07c7-4da3-81c6-31e716f3b15d",
   "metadata": {},
   "source": [
    "## Unweighted Events\n",
    "\n",
    "The information about the distribution is carried by the distribution/number of the events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a43da3e7-5263-4557-8024-cb08ae00f01b",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = 1_000_000\n",
    "integ = 0.0\n",
    "lheevents = []\n",
    "while len(lheevents) < num_samples:\n",
    "    theta = random.uniform(theta_min, theta_max)\n",
    "    phi = random.uniform(phi_min, phi_max)\n",
    "    sig = (\n",
    "        (phi_max - phi_min)\n",
    "        * (theta_max - theta_min)\n",
    "        * dsigma(s, theta, phi)\n",
    "        * math.sin(theta)\n",
    "    )\n",
    "    if sig / maximum > random.uniform(\n",
    "        0, 1\n",
    "    ):  # Pick events randomly according to their contribution\n",
    "        e = pylhe.LHEEvent(\n",
    "            eventinfo=pylhe.LHEEventInfo(\n",
    "                nparticles=5,\n",
    "                pid=0,\n",
    "                weight=result / num_samples,  # The same weight per event\n",
    "                scale=s,\n",
    "                aqed=alpha,\n",
    "                aqcd=aqcd,\n",
    "            ),\n",
    "            particles=[\n",
    "                pylhe.LHEParticle(\n",
    "                    id=11,\n",
    "                    status=-1,\n",
    "                    mother1=0,\n",
    "                    mother2=0,\n",
    "                    color1=0,\n",
    "                    color2=0,\n",
    "                    px=0.0,\n",
    "                    py=0.0,\n",
    "                    pz=EB,\n",
    "                    e=EB,\n",
    "                    m=0.0,\n",
    "                    lifetime=0,\n",
    "                    spin=9.0,\n",
    "                ),\n",
    "                pylhe.LHEParticle(\n",
    "                    id=-11,\n",
    "                    status=-1,\n",
    "                    mother1=0,\n",
    "                    mother2=0,\n",
    "                    color1=0,\n",
    "                    color2=0,\n",
    "                    px=0.0,\n",
    "                    py=0.0,\n",
    "                    pz=-EB,\n",
    "                    e=EB,\n",
    "                    m=0.0,\n",
    "                    lifetime=0,\n",
    "                    spin=9.0,\n",
    "                ),\n",
    "                pylhe.LHEParticle(\n",
    "                    id=22,\n",
    "                    status=2,\n",
    "                    mother1=1,\n",
    "                    mother2=2,\n",
    "                    color1=0,\n",
    "                    color2=0,\n",
    "                    px=0.0,\n",
    "                    py=0.0,\n",
    "                    pz=EB - EB,\n",
    "                    e=EB + EB,\n",
    "                    m=0.0,\n",
    "                    lifetime=0,\n",
    "                    spin=9.0,\n",
    "                ),\n",
    "                pylhe.LHEParticle(\n",
    "                    id=13,\n",
    "                    status=1,\n",
    "                    mother1=3,\n",
    "                    mother2=3,\n",
    "                    color1=0,\n",
    "                    color2=0,\n",
    "                    px=EB * math.sin(theta) * math.cos(phi),\n",
    "                    py=EB * math.sin(theta) * math.sin(phi),\n",
    "                    pz=EB * math.cos(theta),\n",
    "                    e=EB,\n",
    "                    m=0,\n",
    "                    lifetime=0,\n",
    "                    spin=9.0,\n",
    "                ),\n",
    "                pylhe.LHEParticle(\n",
    "                    id=-13,\n",
    "                    status=1,\n",
    "                    mother1=3,\n",
    "                    mother2=3,\n",
    "                    color1=0,\n",
    "                    color2=0,\n",
    "                    px=-EB * math.sin(theta) * math.cos(phi),\n",
    "                    py=-EB * math.sin(theta) * math.sin(phi),\n",
    "                    pz=-EB * math.cos(theta),\n",
    "                    e=EB,\n",
    "                    m=0,\n",
    "                    lifetime=0,\n",
    "                    spin=9.0,\n",
    "                ),\n",
    "            ],\n",
    "        )\n",
    "        lheevents.append(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a4a8faae-4025-44de-9513-821a489b5a7e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<event>\n",
      "  5      0  1.4779213461e-15  1.7472400000e+05  7.8492935636e-03  1.0750000000e-01\n",
      "   11  -1   0   0   0   0  0.00000000e+00  0.00000000e+00  2.09000000e+02  2.09000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "  -11  -1   0   0   0   0  0.00000000e+00  0.00000000e+00 -2.09000000e+02  2.09000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "   22   2   1   2   0   0  0.00000000e+00  0.00000000e+00  0.00000000e+00  4.18000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "   13   1   3   3   0   0  1.23644106e+02  4.11563954e+01 -1.63399162e+02  2.09000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "  -13   1   3   3   0   0 -1.23644106e+02 -4.11563954e+01  1.63399162e+02  2.09000000e+02  0.00000000e+00  0.0000e+00  9.0000e+00\n",
      "</event>\n"
     ]
    }
   ],
   "source": [
    "lheinit = pylhe.LHEInit(\n",
    "    initInfo=pylhe.LHEInitInfo(\n",
    "        beamA=11,\n",
    "        beamB=-11,\n",
    "        energyA=EB,\n",
    "        energyB=EB,\n",
    "        PDFgroupA=0,\n",
    "        PDFgroupB=0,\n",
    "        PDFsetA=0,\n",
    "        PDFsetB=0,\n",
    "        weightingStrategy=3,\n",
    "        numProcesses=1,\n",
    "    ),\n",
    "    procInfo=[\n",
    "        pylhe.LHEProcInfo(\n",
    "            xSection=result,\n",
    "            error=0,\n",
    "            unitWeight=1,\n",
    "            procId=1,\n",
    "        )\n",
    "    ],\n",
    "    generators=[\n",
    "        pylhe.LHEGenerator(\n",
    "            name=\"pylhe\", version=pylhe.__version__, description=\"Monte Carlo Example\"\n",
    "        )\n",
    "    ],\n",
    ")\n",
    "lhef = pylhe.LHEFile(lheinit, lheevents)\n",
    "print(lhef.events[0].tolhe(rwgt=False, weights=False))\n",
    "lhef.tofile(\"unweighted.lhe\", rwgt=False, weights=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e840fba3-e797-497c-9ad8-e48451928614",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html>\n",
       "<div style=\"display:flex; align-items:center;\">\n",
       "<div style=\"width:290px;\">\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"-10 -105 270 120\">\n",
       "<line x1=\"-5\" y1=\"0\" x2=\"255\" y2=\"0\" style=\"fill:none;stroke-width:2;stroke:currentColor\"/>\n",
       "<text text-anchor=\"middle\" x=\"0\" y=\"15\" style=\"fill:currentColor;\">\n",
       "0\n",
       "</text>\n",
       "<text text-anchor=\"middle\" x=\"250\" y=\"15\" style=\"fill:currentColor;\">\n",
       "3.14\n",
       "</text>\n",
       "<text text-anchor=\"middle\" x=\"125.0\" y=\"15\" style=\"fill:currentColor;\">\n",
       "Axis 0\n",
       "</text>\n",
       "<polyline points=\"  0,0   0,-14.4 12.5,-14.4 12.5,-41.5  25,-41.5  25,-65.2 37.5,-65.2 37.5,-82.9  50,-82.9  50,-94.1 62.5,-94.1 62.5,-99  75,-99  75,-100 87.5,-100 87.5,-97.3 100,-97.3 100,-94.6 112.5,-94.6 112.5,-92.3 125,-92.3 125,-92.5 137.5,-92.5 137.5,-94.5 150,-94.5 150,-97.9 162.5,-97.9 162.5,-99.7 175,-99.7 175,-99.8 187.5,-99.8 187.5,-93.9 200,-93.9 200,-83.3 212.5,-83.3 212.5,-64.8 225,-64.8 225,-41.4 237.5,-41.4 237.5,-14.4 250,-14.4 250,0\" style=\"fill:none; stroke:currentColor;\"/>\n",
       "</svg>\n",
       "</div>\n",
       "<div style=\"flex=grow:1;\">\n",
       "Regular(20, 0, 3.14159, label='Axis 0')<br/>\n",
       "<hr style=\"margin-top:.2em; margin-bottom:.2em;\"/>\n",
       "Weight() Σ=WeightedSum(value=1.47792e-09, variance=2.18425e-24)\n",
       "\n",
       "</div>\n",
       "</div>\n",
       "</html>"
      ],
      "text/plain": [
       "Hist(Regular(20, 0, 3.14159, label='Axis 0'), storage=Weight()) # Sum: WeightedSum(value=1.47792e-09, variance=2.18425e-24)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr = pylhe.to_awkward(\n",
    "    pylhe.LHEFile.fromfile(\"unweighted.lhe\", with_attributes=True).events\n",
    ")\n",
    "hist2 = hist.Hist.new.Reg(20, 0, math.pi).Weight()\n",
    "hist2.fill(arr.particles.vector[:, -1].theta, weight=arr.eventinfo.weight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c9ab56a1-50bc-405a-b7e5-c7d6b61881de",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7fec1e55a350>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAHACAYAAAD+yCF8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHjZJREFUeJzt3XuQ1XX9+PHXYdldIAHRQiA30FIISUOwRNEsEMRL1nTRRonsMl7AS0wllNlI6sqMF0oLByPUvECFmFOZoHHJUSsQ1KKwTIJShjENFJiFXT6/P/qx437ZXTiH99nds+fxmDl/7Gc/53Pe58ObPc/zOZdPLsuyLAAAEujS3gMAADoPYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACTTbmGxYsWKOOecc2LAgAGRy+Xi4Ycfbvfbe+ihh2L8+PHxzne+M3K5XKxZs6aoYwKAzqbdwmLbtm1x3HHHxR133NFhbm/btm1x8sknx0033dQmYwKAzqZre93whAkTYsKECS3+fufOnXHNNdfE/fffH//9739j2LBhMXPmzDjttNOKcnsRERMnToyIiPXr1xd0GwBQ7totLPbloosuivXr18f8+fNjwIABsWjRojjjjDPihRdeiKOOOqq9hwcANKNDvnnzpZdeigcffDB+9rOfxSmnnBLvfe9742tf+1qMHj065s2b197DAwBa0CHD4tlnn40sy+Loo4+Ogw46qPGyfPnyeOmllyLify9X5HK5Vi9Tpkxp53sCAOWlQ74Usnv37qioqIhVq1ZFRUVFk98ddNBBERHx7ne/O/7yl7+0up0+ffoUbYwAwN46ZFgMHz48GhoaYvPmzXHKKac0u05lZWUMGTKkjUcGALSm3cLirbfeir///e+NP7/88suxZs2aOOSQQ+Loo4+OCy64ID7/+c/HLbfcEsOHD4/XXnstfvvb38YHPvCBOPPMM5Pe3nve856IiHj99ddjw4YN8corr0RExLp16yIiol+/ftGvX78DubsAUB6ydrJ06dIsIva6TJo0KcuyLNu5c2d27bXXZoMGDcoqKyuzfv36ZZ/85Cez559/vii3l2VZNm/evGbX+c53vnPgdxgAykAuy7KsHXoGAOiEOuSnQgCA0iQsAIBk2vzNm7t3745XXnklevbsGblcrq1vHgAoQJZl8eabb8aAAQOiS5eWj0u0eVi88sorUVNT09Y3CwAksHHjxjj88MNb/H2bh0XPnj0j4n8D69WrV1vfPABQgK1bt0ZNTU3j43hL2jws9rz80atXL2EBACVmX29j8OZNACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAybX7adGhLWZbFjl0NRb+d7pUV+zyVMEA5EBZ0ajt2NcTQax8r+u2snTE+elT578T/tFXQFotQ5kD4SwiQWFsFbbEIZQ6EmUPZWHnN2OhRVdHqOtt31sfI65/4/+uPafWP6/adDTHy+seTjhGg1AkLOoR8Dh3n++C/R4+qin0+C+tR1TXW33TWfo0D9kfqoC0WoUwqwoIOodBDx3v+GENHJWgpNz5uCgAk44gFHc7+HDouRPfK9NuktBXr0xtvfwkOyo2woMPZn0PHkEKpf3oDOiIvhQAAyXhaCAns69D3gbzr35cVtY19vQRX6L+hl+AoN8ICEsjnY3r5fpLFlxW1jX29BOeTG7B/vBQCACTjaRAUqHtlRaydMb4o2/ZlRUCpEhZQoFwu5yUKgP/DSyEAQDLCAgBIRlgAAMkICwAgGe88gw6uWOedKJUv3irW+TwinNMDikFYQAdXrI+dlsoXbzmfB5SWjv9XhQ4j32eO+XwFsmeO0HEU8/9jqRwpo3DCgv12IM8c8/0a63JXrC/fKvUv3irW+TwinNPj7Yo5R0rlSBmF868LHZAv32qe83lAx+cvFwXZ1zPHA+GZI7S9fI+S5ftSZykfKSM/woKC7OuZI1Ba8j1K5ugQLfHIACRRrI+FemMvlBZhAWUq9QO2w91AhLCAsiUCgGIQFkByxfpYqDf2QscnLKCMFOv7MZq7nda+BMkb/6DzEhZQRnw/BlBszm4KACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSySss6uvr45prrokjjjgiunfvHkceeWTMmDEjdu/eXazxAQAlpGs+K8+cOTPuvPPOuOeee+KYY46JlStXxkUXXRS9e/eOK6+8slhjBABKRF5h8fTTT8e5554bZ511VkREDBo0KB588MFYuXJlUQYHAJSWvF4KGT16dDzxxBPx4osvRkTEc889F08++WSceeaZLV6nrq4utm7d2uQCAHROeR2xuPrqq2PLli0xZMiQqKioiIaGhrjhhhvic5/7XIvXqa2tjeuuu+6ABwoAdHx5HbFYsGBB3HffffHAAw/Es88+G/fcc0/cfPPNcc8997R4nenTp8eWLVsaLxs3bjzgQQMAHVNeRyy+/vWvx7Rp0+L888+PiIgPfOAD8c9//jNqa2tj0qRJzV6nuro6qqurD3ykAECHl1dYbN++Pbp0aXqQo6KiwsdNAdgv23c27OP39THy+iciImLlNWOiR9X+PUx1r6yIXC53wOPjwOUVFuecc07ccMMN8Z73vCeOOeaYWL16ddx6663xxS9+sVjjowBZlsWOXa3/5y3Evv4gAOzLyOsfz2PdJ/Z73bUzxu93hFBcef0r3H777fHtb387Lrvssti8eXMMGDAgLr744rj22muLNT4KsGNXQwy99rH2HgYAZSivsOjZs2fMmjUrZs2aVaThANDZdK+siLUzxiff7vadDXkdAaFtOG7Uya28Zmz0qKpIvt3ulem3CXROuVzOyxRlxL90J9ejqsJ/aADajLObAgDJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyXRt7wGUqyzLYseuhqJse/vO4mwXAPZFWLSTHbsaYui1j7X3MAA6hWI+oepeWRG5XK5o2+9shAUAJW/k9Y8XbdtrZ4yPHlUeLveXPdUBrLxmbPSoqijKtrtXFme7ANCcvMPi3//+d1x99dXx6KOPxo4dO+Loo4+OuXPnxogRI4oxvrLQo6pCDQPkqXtlRaydMb4o296+s6GoR0E6s7wezd544404+eST46Mf/Wg8+uij0bdv33jppZfi4IMPLtLwAKB5uVzOk7IOKK9/kZkzZ0ZNTU3MmzevcdmgQYNSjwkAKFF5fY/FI488EiNHjozPfOYz0bdv3xg+fHjcddddrV6nrq4utm7d2uQCAHROeYXFP/7xj5g9e3YcddRR8dhjj8Ull1wSV1xxRdx7770tXqe2tjZ69+7deKmpqTngQQMAHVNeYbF79+44/vjj48Ybb4zhw4fHxRdfHF/5yldi9uzZLV5n+vTpsWXLlsbLxo0bD3jQAEDHlFdY9O/fP4YOHdpk2fvf//7YsGFDi9eprq6OXr16NbkAAJ1TXmFx8sknx7p165ose/HFF2PgwIFJBwUAlKa8wuKrX/1qPPPMM3HjjTfG3//+93jggQdizpw5MXny5GKNDwAoIXmFxQknnBCLFi2KBx98MIYNGxbf/e53Y9asWXHBBRcUa3wAQAnJ+5tFzj777Dj77LOLMRYAoMTldcQCAKA1wgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkMwBhUVtbW3kcrm46qqrEg0HAChlBYfFH//4x5gzZ04ce+yxKccDAJSwgsLirbfeigsuuCDuuuuu6NOnT+oxAQAlqqCwmDx5cpx11lkxduzYfa5bV1cXW7dubXIBADqnrvleYf78+fHss8/GH//4x/1av7a2Nq677rq8BwYAlJ68wmLjxo1x5ZVXxuLFi6Nbt277dZ3p06fH1KlTG3/eunVr1NTU5DfKdpRlWezY1ZB8u9t3pt8mALS3vMJi1apVsXnz5hgxYkTjsoaGhlixYkXccccdUVdXFxUVFU2uU11dHdXV1WlG2w527GqIodc+1t7DAKCdFPOJYPfKisjlckXbfnvIKyzGjBkTL7zwQpNlF110UQwZMiSuvvrqvaICAErdyOsfL9q2184YHz2q8n5XQoeW173p2bNnDBs2rMmyd7zjHXHooYfutbwzWnnN2OhRlT6eulcKMgA6h86VSUXWo6qi05UlAHvrXlkRa2eML8q2t+9sKOpRkPZ2wI+Sy5YtSzAMAOg4crmcJ5IFcq4QACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASCavsKitrY0TTjghevbsGX379o1PfOITsW7dumKNDQAoMXmFxfLly2Py5MnxzDPPxJIlS6K+vj7GjRsX27ZtK9b4AIAS0jWflX/zm980+XnevHnRt2/fWLVqVZx66qlJBwYAlJ4Deo/Fli1bIiLikEMOSTIYAKC05XXE4u2yLIupU6fG6NGjY9iwYS2uV1dXF3V1dY0/b926tdCbBAA6uIKPWEyZMiWef/75ePDBB1tdr7a2Nnr37t14qampKfQmAYAOrqCwuPzyy+ORRx6JpUuXxuGHH97qutOnT48tW7Y0XjZu3FjQQAGAji+vl0KyLIvLL788Fi1aFMuWLYsjjjhin9eprq6O6urqggcIAJSOvMJi8uTJ8cADD8QvfvGL6NmzZ2zatCkiInr37h3du3cvygABgNKR10shs2fPji1btsRpp50W/fv3b7wsWLCgWOMDAEpI3i+FAAC0xLlCAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIJmu7T2AFLIsix27Goqy7e07i7NdACjWY0z3yorI5XJF2fa+dIqw2LGrIYZe+1h7DwMA8jLy+seLst21M8ZHj6r2eYj3UggAkEynOGLxdiuvGRs9qiqKsu3ulcXZLgDlo3tlRaydMT75drfvbCjaEZB8dLqw6FFV0W6HfwBgX3K5XKd+nPJSCACQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAywgIASEZYAADJCAsAIBlhAQAkIywAgGSEBQCQjLAAAJIRFgBAMsICAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsAAAkhEWAEAyBYXFD3/4wzjiiCOiW7duMWLEiPjd736XelwAQAnKOywWLFgQV111VXzrW9+K1atXxymnnBITJkyIDRs2FGN8AEAJyTssbr311vjSl74UX/7yl+P9739/zJo1K2pqamL27NnFGB8AUELyCoudO3fGqlWrYty4cU2Wjxs3Lp566qmkAwMASk/XfFZ+7bXXoqGhIQ477LAmyw877LDYtGlTs9epq6uLurq6xp+3bNkSERFbt27Nd6wt2r6zPnbXbW/cbn1VXncLAEpesR8L9zxuZ1nW6noF3Woul2vyc5Zley3bo7a2Nq677rq9ltfU1BRy0/vUf1ZRNgsAJaOYj4Vvvvlm9O7du8Xf5xUW73znO6OiomKvoxObN2/e6yjGHtOnT4+pU6c2/rx79+54/fXX49BDD20xRgqxdevWqKmpiY0bN0avXr2SbbezsH9aZ/+0zL5pnf3TOvundaW0f7IsizfffDMGDBjQ6np5hUVVVVWMGDEilixZEp/85Ccbly9ZsiTOPffcZq9TXV0d1dXVTZYdfPDB+dxsXnr16tXh/3Hak/3TOvunZfZN6+yf1tk/rSuV/dPakYo98n4pZOrUqTFx4sQYOXJkjBo1KubMmRMbNmyISy65pKBBAgCdR95hcd5558V//vOfmDFjRrz66qsxbNiw+PWvfx0DBw4sxvgAgBJS0Js3L7vssrjssstSj+WAVFdXx3e+8529Xnbhf+yf1tk/LbNvWmf/tM7+aV1n3D+5bF+fGwEA2E9OQgYAJCMsAIBkhAUAkExJhUW+p2tfvnx5jBgxIrp16xZHHnlk3HnnnW000vaRz/5ZtmxZ5HK5vS5//etf23DEbWPFihVxzjnnxIABAyKXy8XDDz+8z+uU09zJd/+U09ypra2NE044IXr27Bl9+/aNT3ziE7Fu3bp9Xq9c5k8h+6ec5s/s2bPj2GOPbfyOilGjRsWjjz7a6nU6w9wpmbDI93TtL7/8cpx55plxyimnxOrVq+Ob3/xmXHHFFbFw4cI2HnnbKPR09uvWrYtXX3218XLUUUe10YjbzrZt2+K4446LO+64Y7/WL7e5k+/+2aMc5s7y5ctj8uTJ8cwzz8SSJUuivr4+xo0bF9u2bWvxOuU0fwrZP3uUw/w5/PDD46abboqVK1fGypUr42Mf+1ice+658ec//7nZ9TvN3MlKxIc+9KHskksuabJsyJAh2bRp05pd/xvf+EY2ZMiQJssuvvji7MQTTyzaGNtTvvtn6dKlWURkb7zxRhuMruOIiGzRokWtrlNuc+ft9mf/lOvcybIs27x5cxYR2fLly1tcp5znz/7sn3KeP1mWZX369Ml+9KMfNfu7zjJ3SuKIRSGna3/66af3Wn/8+PGxcuXK2LVrV9HG2h4O5HT2w4cPj/79+8eYMWNi6dKlxRxmySinuXMgynHu7Dk78yGHHNLiOuU8f/Zn/+xRbvOnoaEh5s+fH9u2bYtRo0Y1u05nmTslERaFnK5906ZNza5fX18fr732WtHG2h4K2T/9+/ePOXPmxMKFC+Ohhx6KwYMHx5gxY2LFihVtMeQOrZzmTiHKde5kWRZTp06N0aNHx7Bhw1pcr1znz/7un3KbPy+88EIcdNBBUV1dHZdcckksWrQohg4d2uy6nWXupD1Ze5Hlc7r2ltZvbnlnkc/+GTx4cAwePLjx51GjRsXGjRvj5ptvjlNPPbWo4ywF5TZ38lGuc2fKlCnx/PPPx5NPPrnPdctx/uzv/im3+TN48OBYs2ZN/Pe//42FCxfGpEmTYvny5S3GRWeYOyVxxKKQ07X369ev2fW7du0ahx56aNHG2h4K2T/NOfHEE+Nvf/tb6uGVnHKaO6l09rlz+eWXxyOPPBJLly6Nww8/vNV1y3H+5LN/mtOZ509VVVW8733vi5EjR0ZtbW0cd9xx8b3vfa/ZdTvL3CmJsHj76drfbsmSJXHSSSc1e51Ro0bttf7ixYtj5MiRUVlZWbSxtodC9k9zVq9eHf379089vJJTTnMnlc46d7IsiylTpsRDDz0Uv/3tb+OII47Y53XKaf4Usn+a01nnT3OyLIu6urpmf9dp5k47vWk0b/Pnz88qKyuzuXPnZmvXrs2uuuqq7B3veEe2fv36LMuybNq0adnEiRMb1//HP/6R9ejRI/vqV7+arV27Nps7d25WWVmZ/fznP2+vu1BU+e6f2267LVu0aFH24osvZn/605+yadOmZRGRLVy4sL3uQtG8+eab2erVq7PVq1dnEZHdeuut2erVq7N//vOfWZaZO/nun3KaO5deemnWu3fvbNmyZdmrr77aeNm+fXvjOuU8fwrZP+U0f6ZPn56tWLEie/nll7Pnn38+++Y3v5l16dIlW7x4cZZlnXfulExYZFmW/eAHP8gGDhyYVVVVZccff3yTjzRNmjQp+8hHPtJk/WXLlmXDhw/PqqqqskGDBmWzZ89u4xG3rXz2z8yZM7P3vve9Wbdu3bI+ffpko0ePzn71q1+1w6iLb8/H2/7vZdKkSVmWmTv57p9ymjvN7ZeIyObNm9e4TjnPn0L2TznNny9+8YuNf5Pf9a53ZWPGjGmMiizrvHPH2U0BgGRK4j0WAEBpEBYAQDLCAgBIRlgAAMkICwAgGWEBACQjLACAZIQFAJCMsACSyOVy8fDDD7f3MIB2JiyAiIh46qmnoqKiIs4444yCrv/qq6/GhAkTCr79hQsXxtChQ6O6ujqGDh0aixYtKnhbQPsRFkBERPz4xz+Oyy+/PJ588snYsGFD3tfv169fVFdXF3TbTz/9dJx33nkxceLEeO6552LixInx2c9+Nn7/+98XtD2g/QgLILZt2xY//elP49JLL42zzz477r777ia/nzFjRgwYMCD+85//NC77+Mc/Hqeeemrs3r07Ipq+FLJz586YMmVK9O/fP7p16xaDBg2K2traFm9/1qxZcfrpp8f06dNjyJAhMX369BgzZkzMmjUr9V0FikxYALFgwYIYPHhwDB48OC688MKYN29evP38hN/61rdi0KBB8eUvfzkiIu68885YsWJF/OQnP4kuXfb+M/L9738/HnnkkfjpT38a69ati/vuuy8GDRrU4u0//fTTMW7cuCbLxo8fH0899VSaOwi0ma7tPQCg/c2dOzcuvPDCiIg444wz4q233oonnngixo4dGxERFRUVcd9998UHP/jBmDZtWtx+++0xZ86cGDhwYLPb27BhQxx11FExevToyOVyLa63x6ZNm+Kwww5rsuywww6LTZs2Jbh3QFtyxALK3Lp16+IPf/hDnH/++RER0bVr1zjvvPPixz/+cZP1jjzyyLj55ptj5syZcc4558QFF1zQ4ja/8IUvxJo1a2Lw4MFxxRVXxOLFi/c5jlwu1+TnLMv2WgZ0fI5YQJmbO3du1NfXx7vf/e7GZVmWRWVlZbzxxhvRp0+fxuUrVqyIioqKWL9+fdTX10fXrs3/CTn++OPj5ZdfjkcffTQef/zx+OxnPxtjx46Nn//8582u369fv72OTmzevHmvoxhAx+eIBZSx+vr6uPfee+OWW26JNWvWNF6ee+65GDhwYNx///2N6y5YsCAeeuihWLZsWWzcuDG++93vtrrtXr16xXnnnRd33XVXLFiwIBYuXBivv/56s+uOGjUqlixZ0mTZ4sWL46STTjrwOwm0KUcsoIz98pe/jDfeeCO+9KUvRe/evZv87tOf/nTMnTs3pkyZEv/617/i0ksvjZkzZ8bo0aPj7rvvjrPOOismTJgQJ5544l7bve2226J///7xwQ9+MLp06RI/+9nPol+/fnHwwQc3O44rr7wyTj311Jg5c2ace+658Ytf/CIef/zxePLJJ4txt4EicsQCytjcuXNj7Nixe0VFRMSnPvWpWLNmTaxatSq+8IUvxIc+9KGYMmVKREScfvrpMWXKlLjwwgvjrbfe2uu6Bx10UMycOTNGjhwZJ5xwQqxfvz5+/etfN/sJkoiIk046KebPnx/z5s2LY489Nu6+++5YsGBBfPjDH057h4Giy2Vv/0wZAMABcMQCAEhGWAAAyQgLACAZYQEAJCMsAIBkhAUAkIywAACSERYAQDLCAgBIRlgAAMkICwAgGWEBACTz/wCm3WBJir4ghgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist2.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c9be8352-6ab8-45c2-a4ac-fa458ee41652",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(([StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7fecb3a89f90>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)],\n",
       "  [StairsArtists(stairs=<matplotlib.patches.StepPatch object at 0x7fecb3a89e50>, errorbar=<ErrorbarContainer object of 3 artists>, legend_artist=<ErrorbarContainer object of 3 artists>)]),\n",
       " RatioErrorbarArtists(line=<matplotlib.lines.Line2D object at 0x7fecb3a8aad0>, errorbar=<ErrorbarContainer object of 3 artists>))"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAHACAYAAABaopmvAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAANfBJREFUeJzt3Xl4VNXh//HPJGQnBAQFAoEg+yZEwAqCiiCbouC3Lv0iIlUfRZHt8acsVgTRSOtCFcVKEdQWQQTEb13KHrRIK4EIioIomMhiZMtK9vP7g2YkZJuZ5GZm7rxfzzPPk7nLOeee3Mx8cu7mMMYYAQAA+LkgbzcAAACgNhBqAACALRBqAACALRBqAACALRBqAACALRBqAACALRBqAACALRBqAACALRBqAACALRBqAACALQRkqNm2bZtGjhyp2NhYORwOvf/++16vb82aNRo6dKiaNGkih8OhlJQUS9sEAIDdBGSoycnJUY8ePbRw4UKfqS8nJ0dXXXWVnn322TppEwAAdlPP2w3whuHDh2v48OGVzi8oKNDjjz+uv//97zpz5oy6deum+fPn69prr7WkPkkaO3asJOnw4cMe1QEAQKALyFBTnfHjx+vw4cNasWKFYmNjtXbtWg0bNkx79+5V+/btvd08AABQgYA8/FSV77//Xu+8845WrVqlAQMGqG3btnrkkUfUv39/LV261NvNAwAAlSDUXGDXrl0yxqhDhw6qX7++85WUlKTvv/9e0rlDRA6Ho8rXxIkTvbwlAAAEFg4/XaCkpETBwcFKTk5WcHBwmXn169eXJLVo0ULffPNNleU0atTIsjYCAIDyCDUXSEhIUHFxsdLT0zVgwIAKlwkJCVGnTp3quGUAAKAqARlqsrOzdfDgQef7Q4cOKSUlRRdddJE6dOigMWPG6K677tLzzz+vhIQEnThxQps3b1b37t01YsSIWq2vVatWkqRTp04pNTVVR48elSTt379fktSsWTM1a9asJpsLAEBAcBhjjLcbUde2bt2qgQMHlps+btw4LVu2TIWFhZo3b57eeustHTlyRI0bN1bfvn01Z84cde/evdbrk6Rly5Zp/Pjx5ZaZPXu2nnzySbfrBAAg0ARkqAEAAPbD1U8AAMAWCDUAAMAWAupE4ZKSEh09elTR0dFyOBzebg4AAHCBMUZZWVmKjY1VUFDl4zEBFWqOHj2quLg4bzcDAAB4IC0tTS1btqx0fkCFmujoaEnnOqVBgwZebg0AAHBFZmam4uLinN/jlQmoUFN6yKlBgwaEGgAA/Ex1p45wojAAALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALCFgHpMAlCrjJEKc+umrpBIiSfLo4aMMTpbWOzWOrkFReo9b5MkaefjgxQZ6vrXRkRIcLW3tQdqE6EGtuTJh7fbCnIU+Vwra+soNfOoFBpVN3XBts4WFKnX7A/cWsdIksIlSf3nfSh3Ikry44PdCkEeIfDjPIQa2NLZwmJ1eeKfltYRoTx9E25pFU7GGLe+TIAKFebqm/Dfu7VKrglVpKOg3M8uec6tqjxD4Md5CDWwJ2MUoTzXF5eU99//RsOV51KAiFS+8+deeYuUq7BarSNS+UoOnyDpXEiLrLp4wBLnhxi3Ak0dIfDjfIQa2JMH/5HWRK7CdFauD9vkubEsYIXcyd8qMqpB7ZdbUKRe8za6vLyn/1AQ+FERQg1QC5IfH1zrQ+C52ZnSS//9uaBYKiiq1fIvxEmdASYk0qLDNkVuBfzzEfZRU4Qa2J5V/5GeL9KKkxVDg50/DvjjFo+/KFy1b+5Q60/qROU8uZquIFd6rt25nx85KIVGVr18HVytFxESrH1zh1pax/mBHzgfn2CwP8v+I60755+/UxVPhvJ/Xdm43S7UosJc6ZlYz9cvDTdVqCby1AqHw2F9OD4v8APnI9QAPioi5NcP7tLzB6yUW5gqhcVYXg8qxgmvQM0RagAfxfktgeVsYbFzJMWVq+mkmo3MJYfUxbgNULcINYCvCok8dw8OC+XmZCryz50srQPuc/dqOsmDk2wJzbAhQg3gqxwO688FsviKKnjm00cHKrK+tSe3n394E7ALQg0A+JjI0GCuRAM8wFO6AQCALfCvAABJ0tmcLMvriIiMliPIT/+Xsvqp7HX1xHfAxgg1qHOmpERnc639Aj2bk1Un9+Swk8aLurq0XE0ecJj7SKoi6/vpZeM1vY9MNdhfgZoj1KDOnc3NUuRzrVxe3pMvUb4ggABRmCsVWPxVZsUdw2EJQg0QwCIio5X7SKrb6+VW8nNFzuZkOUeB/PkZVuffHC938rfnvuiqkFtQrN5//FyStPPRvoqs5i64uQXFGvDHLZK4h4w76uKWBGbGETnC6lteD2qOUAOvOjnha0VERVe7nDtfoheKiKy+/EDlCAqq08NBdfEMq52PD642QHgiN6dATf77c6/5293ajtJwU73/lsmoQK2qySFTiSeB+xNCDbwqIiraf8+xgE/qPW+jJeVGKE/f8BBp3xASqc55b7i8uCd3Xo5Ufp08ngS1i1ADwFJlnmH1WL9qD9t44vxDN67w9EuulNU3x+PGeFWLCK2n5Lk3W1oHTwL3T4QaAJY6//wWq85/iJTcGkWp6eGIxlGhcnBzPK/hSeCojJ/eMAIAPHd+iHE30Eg8bBTwVfyrAcBadfBgzjrH1UmATyLUALBWXTyYEwDE4ScAAGAThBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALfhNqioqK9Pjjj6tNmzaKiIjQpZdeqrlz56qkpMTbTQMAAD6gnrcb4Kr58+frtdde05tvvqmuXbtq586dGj9+vGJiYjR58mRvNw8AAHiZ34Sazz//XDfffLNuuOEGSVJ8fLzeeecd7dy508stAwAAvsBvDj/1799fmzZt0oEDByRJX375pT777DONGDGi0nXy8/OVmZlZ5gUAAOzJb0ZqHnvsMWVkZKhTp04KDg5WcXGxnn76af3ud7+rdJ3ExETNmTOnDlsJAAC8xW9GalauXKm//e1vWr58uXbt2qU333xTzz33nN58881K15kxY4YyMjKcr7S0tDpsMQAAqEt+M1Lz//7f/9P06dN1xx13SJK6d++uH3/8UYmJiRo3blyF64SFhSksLKwumwkAsJmzOVmWlh8RGS1HkN+MMfg0vwk1ubm5Crrglx4cHMwl3bXMlJTobK61f8Bnc7IUaWkNAFB7Gi/qWu0yuSZUkY6Ccj+7IveRVEXWj/G4ffiV34SakSNH6umnn1arVq3UtWtX7d69Wy+88IJ+//vfe7tptnI2N0uRz7Vyax13/5gJNAAAK/hNqHn55Zf1hz/8QQ8++KDS09MVGxur+++/X0888YS3mwYAsJmIyGjlPpLq1jq5lfxckbM5WS6NAME9fhNqoqOjtWDBAi1YsMDbTQkYJyd8rYioaJeWdeeP+XwRka6VDwB1yREUxCEhP+Q3oQZ1LyIqmj9qAIDf4HRrAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC4QaAABgC/W83QC4wRipMNfaOqwuHwAAixBq/IgpyJEjsYWldURaWjoAoJzCXKnA4q/jkEjJ4bC2Dh/gV6HmyJEjeuyxx/Txxx/r7Nmz6tChg5YsWaJevXp5u2l14mxhsduhI9eEKtJRUO5nAIBviPxzJ8vrMDOOyBFW3/J6vM1vQs3p06d11VVXaeDAgfr44491ySWX6Pvvv1fDhg293TSv6JW3SLkKq3Y5IylP4ZKkcOXJnZyeHMK4DQD4ipr8k3q2sFiR1X9l+D2/CTXz589XXFycli5d6pwWHx/vvQZ52T8fHabI+g0srSMiJNjS8gEgYIVEqnPeG5ZWEal8JYdPsLQOX+M3oeaDDz7Q0KFDdeuttyopKUktWrTQgw8+qPvuu8/bTfOKyNBgRYb6za8PAHCeiNB6Sp57s6V15GZnSi9ZWoXP8ZtvxR9++EGLFi3StGnTNHPmTP3nP//RpEmTFBYWprvuuqvCdfLz85Wfn+98n5mZWVfNBQCgUg6Hw/p/TEMDb7Tdb0JNSUmJevfurWeeeUaSlJCQoK+//lqLFi2qNNQkJiZqzpw5ddlMAADgJX5z873mzZurS5cuZaZ17txZqampla4zY8YMZWRkOF9paWlWNxMAAHiJ34zUXHXVVdq/f3+ZaQcOHFDr1q0rXScsLExhYQFwujcAAPCfkZqpU6dqx44deuaZZ3Tw4EEtX75cr7/+uh566CFvNw0AAPgAvwk1ffr00dq1a/XOO++oW7dueuqpp7RgwQKNGTPG200DAAA+wG8OP0nSjTfeqBtvvNHbzQAAAD7I8pGaXbt2ae/evc7369at06hRozRz5kwVFHDLfgAAUDssDzX333+/Dhw4IOncvWbuuOMORUZGatWqVXr00Uetrh4AAAQIy0PNgQMH1LNnT0nSqlWrdPXVV2v58uVatmyZVq9ebXX1AAAgQFgeaowxKikpkSRt3LhRI0aMkCTFxcXpxIkTVlcPAAAChOWhpnfv3po3b57efvttJSUl6YYbbpAkHTp0SE2bNrW6egAAECAsDzUvvviidu3apYkTJ2rWrFlq166dJOm9995Tv379rK4eAAAECMsv6e7Ro0eZq59K/elPf1K9en51RTkAAPBhlo/UXHrppTp58mS56Xl5eerQoYPV1QMAgABheag5fPiwiouLy03Pz8/XTz/9ZHX1AAAgQFh2/OeDDz5w/vzPf/5TMTExzvfFxcXatGmT2rRpY1X1AAAgwFgWakaNGiVJcjgcGjduXJl5ISEhio+P1/PPP29V9QAAIMBYFmpK703Tpk0bffHFF2rSpIlVVQEAAFh/9dOhQ4esrgIAAKBuntK9adMmbdq0Senp6c4RnFJvvPFGXTQBAADYnOWhZs6cOZo7d6569+6t5s2by+FwWF0lAAAIQJaHmtdee03Lli3T2LFjra4KAAAEMMvvU1NQUMDjEAAAgOUsDzX33nuvli9fbnU1AAAgwFl++CkvL0+vv/66Nm7cqMsuu0whISFl5r/wwgtWNwEAAAQAy0PNnj171LNnT0nSV199VWYeJw0DAIDaYnmo2bJli9VVAAAAWH9ODQAAQF2wfKRm4MCBVR5m2rx5s9VNAAAAAcDyUFN6Pk2pwsJCpaSk6Kuvvir3oEsAAABPWR5qXnzxxQqnP/nkk8rOzra6egAAECC8dk7NnXfeyXOfAABArfFaqPn8888VHh7ureoBAIDNWH746ZZbbinz3hijY8eOaefOnfrDH/5gdfUAACBAWB5qYmJiyrwPCgpSx44dNXfuXA0ZMsTq6gEAQICwPNQsXbrU6ioAAACsDzWlkpOT9c0338jhcKhLly5KSEioq6oBAEAAsDzUpKen64477tDWrVvVsGFDGWOUkZGhgQMHasWKFbr44outbgIAAAgAll/99PDDDyszM1Nff/21Tp06pdOnT+urr75SZmamJk2aZHX1AAAgQFg+UvPJJ59o48aN6ty5s3Naly5d9Morr3CiMAAAqDWWj9SUlJQoJCSk3PSQkBCVlJRYXT0AAAgQloea6667TpMnT9bRo0ed044cOaKpU6dq0KBBVlcPAAAChOWhZuHChcrKylJ8fLzatm2rdu3aqU2bNsrKytLLL79sdfUAACBAWH5OTVxcnHbt2qUNGzbo22+/lTFGXbp00eDBg62uGgAABBDLRmo2b96sLl26KDMzU5J0/fXX6+GHH9akSZPUp08fde3aVZ9++qlV1QMAgABjWahZsGCB7rvvPjVo0KDcvJiYGN1///164YUXrKoeAAAEGMtCzZdffqlhw4ZVOn/IkCFKTk62qnoAABBgLAs1P//8c4WXcpeqV6+efvnlF6uqBwAAAcayUNOiRQvt3bu30vl79uxR8+bNraoeAAAEGMtCzYgRI/TEE08oLy+v3LyzZ89q9uzZuvHGG62qHgAABBjLLul+/PHHtWbNGnXo0EETJ05Ux44d5XA49M033+iVV15RcXGxZs2aZVX1AAAgwFg2UtO0aVNt375d3bp104wZMzR69GiNGjVKM2fOVLdu3fSvf/1LTZs29bj8xMREORwOTZkypfYaDQAA/JalN99r3bq1PvroI50+fVoHDx6UMUbt27dXo0aNalTuF198oddff12XXXZZLbUUAAD4O8sfkyBJjRo1Up8+fXTFFVfUONBkZ2drzJgxWrx4cY3LAgAA9lEnoaY2PfTQQ7rhhhtcesxCfn6+MjMzy7wAAIA9Wf7sp9q0YsUK7dq1S1988YVLyycmJmrOnDkWt+pXpqREZ3OzLCv/bE6WIi0rHQAA/+Y3oSYtLU2TJ0/W+vXrFR4e7tI6M2bM0LRp05zvMzMzFRcXZ1UTdTY3S5HPtbKsfAINAMATZ3Os+4e7VERktBxB3j0A5DehJjk5Wenp6erVq5dzWnFxsbZt26aFCxcqPz9fwcHBZdYJCwtTWFhYXTfVZbkmVJGOgnI/AwBQmxov6urScjX5Xsp9JFWR9WM8al9t8ZtQM2jQoHJ3KB4/frw6deqkxx57rFyg8baTE75WRFR0tcvlVvJzdSIiqy8bAAB3nB9i/PEfbb8JNdHR0erWrVuZaVFRUWrcuHG56b4gIira64kVABC4IiKjlftIqqV1nM3JcnkUqC74TagBAACucwQFBdw/134darZu3ertJgAAAB/hd/epAQAAqAihBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2AKhBgAA2ILfhJrExET16dNH0dHRuuSSSzRq1Cjt37/f280CAAA+wm9CTVJSkh566CHt2LFDGzZsUFFRkYYMGaKcnBxvNw0AAPiAet5ugKs++eSTMu+XLl2qSy65RMnJybr66qu91CoAAOAr/Gak5kIZGRmSpIsuusjLLQEAAL7Ab0ZqzmeM0bRp09S/f39169at0uXy8/OVn5/vfJ+ZmVkXzQMAAF7glyM1EydO1J49e/TOO+9UuVxiYqJiYmKcr7i4uDpqIQAAqGt+F2oefvhhffDBB9qyZYtatmxZ5bIzZsxQRkaG85WWllZHrQQAAHXNbw4/GWP08MMPa+3atdq6davatGlT7TphYWEKCwurg9YBAABv85tQ89BDD2n58uVat26doqOjdfz4cUlSTEyMIiIivNw6AADgbX5z+GnRokXKyMjQtddeq+bNmztfK1eu9HbTAACAD/CbkRpjjLebAAAAfJjfjNQAAABUhVADAABsgVADAABsgVADAABsgVADAABsgVADAABsgVADAABsgVADAABsgVADAABswW/uKAwAAHxYYa5UYFGsKMhxabGADDW52RmqF1T7j104m5OlyFovFQAA3xf5507WFZ7v2nd2QIaayJe7KjLMUeUyuSZUkY6Ccj9XWW6ttA4AAHgiIEONK84PMa4EGgAAAk5IpDrnvWFpFZHK1xbHAy4tG5Ch5tQ9/1a9ZrGW1hERGW1p+QAAeFtEaD0lz73Z0jpyszOlP7m2bECGmvDI+oqsH+PtZgAA4NccDociQy2OEqHBKnJxUS7pBgAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAthCQoebzHf/WyJEjFRsbK4fDoffff9+t9fPy8nT33Xere/fuqlevnkaNGlVumWPHjul///d/1bFjRwUFBWnKlCm10nYAAFCxgAw1ubm56tGjhxYuXOjR+sXFxYqIiNCkSZM0ePDgCpfJz8/XxRdfrFmzZqlHjx41aS4AAHBBQIaaQdcN1Lx583TLLbdUOL+goECPPvqoWrRooaioKP3mN7/R1q1bnfOjoqK0aNEi3XfffWrWrFmFZcTHx+vPf/6z7rrrLsXExFixGQAA4Dz1vN0AXzR+/HgdPnxYK1asUGxsrNauXathw4Zp7969at++vbebBwAAKhCQIzVV+f777/XOO+9o1apVGjBggNq2batHHnlE/fv319KlS73dPAAAUAlGai6wa9cuGWPUoUOHMtPz8/PVuHFjL7UKAABUh1BzgZKSEgUHBys5OVnBwcFl5tWvX99LrQIAANUh1FwgISFBxcXFSk9P14ABA7zdHAAA4KKADDU5OTlKSUlxvj906JBSUlJ00UUXqUOHDhozZozuuusuPf/880pISNCJEye0efNmde/eXSNGjJAk7du3TwUFBTp16pSysrKc5fXs2dNZbum07Oxs/fLLL0pJSVFoaKi6dOlSR1sKAEDgcBhjjLcbUVcyMzMVExOj1atW6n9uvb3c/HHjxmnZsmUqLCzUvHnz9NZbb+nIkSNq3Lix+vbtqzlz5qh79+6Szl2y/eOPP5Yr4/zudDgc5ea3bt1ahw8frr2NAgDAxnKzM1T0dJxins1SRkaGGjRoUOmyfhdqXn31Vf3pT3/SsWPH1LVrVy1YsMDlw0SloebYkVQ1i42zuKUAAKCm3Ak1fnVJ98qVKzVlyhTNmjVLu3fv1oABAzR8+HClpqZ6u2kAAMDL/CrUvPDCC7rnnnt07733qnPnzlqwYIHi4uK0aNEibzcNAAB4md+cKFxQUKDk5GRNnz69zPQhQ4Zo+/btFa6Tn5+v/Px85/vTp09Lko4ePaqikvLnuwAAAN9yNidThRnFks49e7EqfhNqTpw4oeLiYjVt2rTM9KZNm+r48eMVrpOYmKg5c+aUm96rz5WWtBEAAFjn4MGD6tOnT6Xz/SbUlLrwiiJjTIVXGUnSjBkzNG3aNOf71NRUde/eXWlpaVWeaCSdu+w7NjZW0rmRnaioqBq2HACAwObpd2tmZqbi4uLUrl27Kpfzm1DTpEkTBQcHlxuVSU9PLzd6UyosLExhYWHO9w0bNpQkNWjQoNpQc/7dhBs0aECoAQCghmr63Xrhnf4v5DcnCoeGhqpXr17asGFDmekbNmxQv379vNQqAADgK/xmpEaSpk2bprFjx6p3797q27evXn/9daWmpuqBBx7wdtMAAICX+VWouf3223Xy5EnNnTtXx44dU7du3fTRRx+pdevWLq1//qEoAABgL353R+GaKL2jcHV3JJTOncxU+lTu7OxszqkBAKCGPP1udfX722/OqQEAAKiKX4Wabdu2aeTIkYqNjZXD4dD777/v7SYBAAAf4VehJicnRz169NDChQu93RQAAOBj/OpE4eHDh2v48OHebgYAAPBBfhVq3HXhs58yMzO92BoAAGAlvzr85K7ExETFxMQ4X3Fxcd5uEgAAsIitQ82MGTOUkZHhfKWlpXm7SQAAwCK2Pvx04bOfAACAfdl6pAYAAAQOvxqpyc7O1sGDB53vDx06pJSUFF100UVq1aqVF1sGAAC8za9Czc6dOzVw4EDn+2nTpkmSxo0bp2XLlnmpVQAAwBf4Vai59tprFUCPqgIAAG7gnBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAAGALhBoAACBJysnJkcPhkMPhUE5Ojreb4zZCjRf5+84DAIAvqXGo+emnn3TkyJHaaAsAAIDHPAo1JSUlmjt3rmJiYtS6dWu1atVKDRs21FNPPaWSkpLabiMAAEC16nmy0qxZs7RkyRI9++yzuuqqq2SM0b/+9S89+eSTysvL09NPP13b7QQAAKiSR6HmzTff1F//+lfddNNNzmk9evRQixYt9OCDDxJqAABAnfPo8NOpU6fUqVOnctM7deqkU6dO1bhRqD2cjAwA9sDnefU8CjU9evTQwoULy01fuHChevToUeNGAQAAuMujw09//OMfdcMNN2jjxo3q27evHA6Htm/frrS0NH300Ue13UYAAIBqeTRSc8011+jAgQMaPXq0zpw5o1OnTumWW27R/v37NWDAgNpuI3yc1UOiDLkGlrr4fbNPBRb2qcDh0UiNJMXGxnJCMGwjJydH9evXlyRlZ2crKiqKOrxYhx3Y5XdhdR3sT6hNLoeaPXv2qFu3bgoKCtKePXuqXPayyy6rccMAAADc4fLhp549e+rEiRPOnxMSEtSzZ89yr4SEBMsaK0mvvvqq2rRpo/DwcPXq1UuffvqpJfUUFxc7f962bVuZ99RRt3XYYRuow3fKpw7fqsMO20AdvlO+jIsOHz5sSkpKnD9X9bLKihUrTEhIiFm8eLHZt2+fmTx5somKijI//vijS+tnZGQYSSYjI6PK5VavXm1atGhhJDlfLVu2NKtXr66NzaAOHyqfOnyrDjtsA3X4TvnU4Vt11KR8V7+/XQ4150tKSjKFhYXlphcWFpqkpCRPinTJFVdcYR544IEy0zp16mSmT5/u0vqudMrq1auNw+Eo0+mSjMPhMA6Ho1Z+udThG+VTh2/VYYdtoA7fKZ86fKuOmpZvaagJCgoyP//8c7npJ06cMEFBQZ4UWa38/HwTHBxs1qxZU2b6pEmTzNVXX+1SGdV1SlFRkWnZsmW5Tj+/8+Pi4kxRUZHH20EdvlE+dfhWHXbYBurwnfKpw7fqqI3yXQ01Hl39ZIyRw+EoN/3kyZOWnbl+4sQJFRcXq2nTpmWmN23aVMePH69wnfz8fOXn5zvfZ2ZmSpJSUlKcZ9tLUqNGjdSmTRtt2rRJP/30U6VtMMYoLS1NS5YsUe/evSVJ8fHxuuiii/TLL78oLS2tzPLR0dFq3769iouL9eWXX0qSdu7c6VId69at0y233KLTp0/r0KFDZZaJiIhQ586dJUm7d++WMabM/F9++cXt7ZCkkJAQde/eXZK0d+9eFRYWllmvffv2io6O1pEjR/Thhx+6XEefPn2c51p98803Onv2bJll27Rpo0aNGunnn392PvHd1X7aunWrBg0apO+++05ZWVlllomLi9PFF1+sU6dO6fDhw2XmRUVF6dixYx71U/PmzdW8eXNlZmbq4MGDZdYJCwtT165dJZ07uX7Hjh0u1zFs2DC1atVKubm5+vbbb8ssFxQUpJ49e0qS9u3bp7y8POc8V/tq8+bNaty4cbn5PXv2VFBQkA4cOKDs7Owy81q1aqUmTZroH//4h9t91a1bN4WGhuqHH37QmTNnyiwfGxurZs2a6cyZM/rhhx9c3oYlS5bozjvvVGRkpFJTU53n+ZW65JJL1LJlS2VnZ+vAgQNl5u3evdulOt5//321adOmzLzSz4i8vDzt27ev3LqXX365JOnvf/+7W/1U0WfE+bp3766QkBB9//33ysjIkOT67/vTTz/VtddeW+FnROfOnRUREaEff/xRJ0+eLDOvadOm+u677zz62zj/M+Lnn38us07jxo3VunVrnT17Vm+//bZb5Vf0GVEqJiZGbdu2VWFhofbu3euc7k4/XXbZZRV+RnTs2FGStGvXrnLrd+nSxa2/79J+uvAzoqioqMw6HTp0UP369fXTTz8pPT3d5e1YunSp7r33XknlPyMk6dJLL1XDhg11/PhxHT16tMy8b7/91qU6li9f7mx7qdLPiBMnTig1NbXMvPr166tDhw5KSkpyu58u/Iy48LOpqsJcNnr0aDN69GgTFBRkRowY4Xw/evRoc9NNN5n4+HgzdOhQd4p02ZEjR4wks3379jLT582bZzp27FjhOrNnz640GZ7/GjNmjDHGmBdeeMGl5c9/vf3228YYYxYuXFhu3pAhQ4wxvyZMd16lbXr33XfLzUtISHBuY2hoaLn5f/zjH92uT5Jp0aKFs9wLj3tKMlu2bDHGGDN9+nS3yg0NDXWWm5CQUG7+u+++a4wx5vnnn3e7zUuWLDHGGDNkyJBy8xYuXGiMMebtt98uN+/KK680y5cv96ifZs+ebYwx5pNPPik3r23bts5tbdKkiVvlPvjgg8YYY5KTk8vNi46OdpbbpUsXj9r98ssvVzg9Ly/PGGPMNddcU27e4sWLjTHG3HvvvW7Xl5aWZowx5re//W25ec8884wxxph169a5XW5ycrIxxpgHH3yw3LypU6caY4zZvn17hX3oSvm33HJLuWmlf4/fffddheuUateunVvbUt1nRHp6ujHGmJEjR7rdT8uXL6/0M+Krr74yxhhzzz33lJs3ffp0j/82qvqMuOeee4wxxnz11Vdul1vVZ8TIkSONMcakp6d71Obly5dX+hlRqqL1vvvuO4/6qbrPiNLvuKlTp7pVbnh4eJWfEevWrTPGGPPMM8+Um3fFFVe4VEfnzp3LTSv9jFi8eHG5eddcc40xxpg333zT7X6q7DOiupEax39/YS4ZP368pHMPtLztttsUERHhnBcaGqr4+Hjdd999atKkiatFuqygoECRkZFatWqVRo8e7Zw+efJkpaSkKCkpqdw6FY3UxMXFKSkpqcKRmvXr12vo0KHVtuUvf/lLjUZq7r///mrrWL16dY1GaoYNG+bWdkjuj9S4sh1/+ctfPB6pcaX8jRs31mikZuDAgS5tQ01Galztp5qM1LhSx/r16z0eqVm3bp1GjRrl0nZ4OlLjaj/VZKSm9D/Zqrz33nsej9S89dZbGjdunEvbUZORGlf6asuWLTUaqfHkb8OdkRpXf981GalxtZ9qMlLjbj95MlLjynYsXry4RiM1Y8aMqbaOt956y6ORms2bN2vQoEHVln9+P1U0UnPNNdcoIyNDDRo0qLyQKiNPJZ588kmTnZ3tyao1csUVV5gJEyaUmda5c+daO1G49LhfRSczSbV77JI67L8N1OE75VOHb9Vhh22gjrot39IThb2l9JLuJUuWmH379pkpU6aYqKgoly8jd+fqpws734qzzKnD/ttAHb5TPnX4Vh122AbqqLvyLQ81q1atMrfeeqv5zW9+YxISEsq8rPTKK6+Y1q1bm9DQUHP55Ze7dQm5O/epufBM7bi4uFq/HwB1eL986vCtOuywDdThO+VTh2/VUZPyXf3+duucmlIvvfSSZs2apXHjxmnx4sUaP368vv/+e33xxRd66KGHfPaZUJmZmYqJian+mJzO3fXw008/1bFjx9S8eXMNGDBAwcHBtdoe6vCN8qnDt+qwwzZQh++UTx2+VYen5bv6/e1RqOnUqZNmz56t3/3ud4qOjtaXX36pSy+9VE888YROnTqlhQsXultknXAn1AAAAN/g6ve3y89+Ol9qaqr69esn6dyVOKVXnYwdO1bvvPOOJ0UCAADUiEehplmzZs7LAFu3bq0dO3ZIkg4dOlTu0kEAAIC64FGoue666/R///d/kqR77rlHU6dO1fXXX6/bb7+9zD1kAAAA6opH59SUlJSopKRE9eqde8rCu+++q88++0zt2rXT6NGjFRcXV+sNrQ2cUwMAgP+x9JyaoKAgZ6CRpNtuu00zZ87Ud999pw4dOnhSJAAAQI24FWrOnDmjMWPG6OKLL1ZsbKxeeukllZSU6IknnlDbtm21Y8cOvfHGG1a1FQAAoFJuPaV75syZ2rZtm8aNG6dPPvlEU6dO1SeffKK8vDx99NFHuuaaa6xqJwAAQJXcCjUffvihli5dqsGDB+vBBx9Uu3bt1KFDBy1YsMCi5gEAALjGrcNPR48eVZcuXSSde+JneHi4S0+9BQAAsJpboaakpEQhISHO98HBwYqKiqr1RgEAALjLrcNPxhjdfffdCgsLkyTl5eXpgQceKBds1qxZU3stBAAAcIFboWbcuHFl3t9555212hgAAABPuRVqli5dalU7AAAAasSjm+8BAAD4GkINAACwBUINAACwBUINAACwBUINAACwBUINAACwBUINAACwBUINAACwBUINAACwBb8JNU8//bT69eunyMhINWzY0NvNAQAAPsZvQk1BQYFuvfVWTZgwwdtNAQAAPsitZz9505w5cyRJy5Yt825DAACAT/KbkRoAAICq+M1IjSfy8/OVn5/vfJ+ZmenF1gAAACt5daTmySeflMPhqPK1c+dOj8tPTExUTEyM8xUXF1eLrQcAAL7EYYwx3qr8xIkTOnHiRJXLxMfHKzw83Pl+2bJlmjJlis6cOVNt+RWN1MTFxSkjI0MNGjTwuN0AAKDuZGZmKiYmptrvb68efmrSpImaNGliWflhYWEKCwuzrHwAAOA7/OacmtTUVJ06dUqpqakqLi5WSkqKJKldu3aqX7++dxsHAAC8zm9CzRNPPKE333zT+T4hIUGStGXLFl177bVeahUAAPAVXj2npq65ekwOAAD4Dle/v7lPDQAAsAVCDQAAsAVCDQAAsAVCDQAAsAVCDQAAsAW/uaS7NpRe6MUzoAAA8B+l39vVXbAdUKHm5MmTksQzoAAA8ENZWVmKiYmpdH5AhZqLLrpI0rm7E1fVKSiv9LlZaWlp3OPHTfSd5+g7z9F3nqHfPGdl3xljlJWVpdjY2CqXC6hQExR07hSimJgYdlYPNWjQgL7zEH3nOfrOc/SdZ+g3z1nVd64MRnCiMAAAsAVCDQAAsIWACjVhYWGaPXu2wsLCvN0Uv0PfeY6+8xx95zn6zjP0m+d8oe8C6oGWAADAvgJqpAYAANgXoQYAANgCoQYAANiC7ULNq6++qjZt2ig8PFy9evXSp59+WuXySUlJ6tWrl8LDw3XppZfqtddeq6OW+h53+m7r1q1yOBzlXt9++20dttj7tm3bppEjRyo2NlYOh0Pvv/9+teuwz53jbt+xz52TmJioPn36KDo6WpdccolGjRql/fv3V7se+51nfcd+d86iRYt02WWXOe9B07dvX3388cdVruONfc5WoWblypWaMmWKZs2apd27d2vAgAEaPny4UlNTK1z+0KFDGjFihAYMGKDdu3dr5syZmjRpklavXl3HLfc+d/uu1P79+3Xs2DHnq3379nXUYt+Qk5OjHj16aOHChS4tzz73K3f7rlSg73NJSUl66KGHtGPHDm3YsEFFRUUaMmSIcnJyKl2H/e4cT/quVKDvdy1bttSzzz6rnTt3aufOnbruuut088036+uvv65wea/tc8ZGrrjiCvPAAw+UmdapUyczffr0Cpd/9NFHTadOncpMu//++82VV15pWRt9lbt9t2XLFiPJnD59ug5a5x8kmbVr11a5DPtcxVzpO/a5iqWnpxtJJikpqdJl2O8q5krfsd9VrlGjRuavf/1rhfO8tc/ZZqSmoKBAycnJGjJkSJnpQ4YM0fbt2ytc5/PPPy+3/NChQ7Vz504VFhZa1lZf40nflUpISFDz5s01aNAgbdmyxcpm2gL7XM2xz5WVkZEh6ddn21WE/a5irvRdKfa7XxUXF2vFihXKyclR3759K1zGW/ucbULNiRMnVFxcrKZNm5aZ3rRpUx0/frzCdY4fP17h8kVFRTpx4oRlbfU1nvRd8+bN9frrr2v16tVas2aNOnbsqEGDBmnbtm110WS/xT7nOfa58owxmjZtmvr3769u3bpVuhz7XXmu9h373a/27t2r+vXrKywsTA888IDWrl2rLl26VList/Y52z3Q0uFwlHlvjCk3rbrlK5oeCNzpu44dO6pjx47O93379lVaWpqee+45XX311Za209+xz3mGfa68iRMnas+ePfrss8+qXZb9rixX+4797lcdO3ZUSkqKzpw5o9WrV2vcuHFKSkqqNNh4Y5+zzUhNkyZNFBwcXG5kIT09vVxaLNWsWbMKl69Xr54aN25sWVt9jSd9V5Err7xS3333XW03z1bY52pXIO9zDz/8sD744ANt2bJFLVu2rHJZ9ruy3Om7igTqfhcaGqp27dqpd+/eSkxMVI8ePfTnP/+5wmW9tc/ZJtSEhoaqV69e2rBhQ5npGzZsUL9+/Spcp2/fvuWWX79+vXr37q2QkBDL2uprPOm7iuzevVvNmzev7ebZCvtc7QrEfc4Yo4kTJ2rNmjXavHmz2rRpU+067HfneNJ3FQnE/a4ixhjl5+dXOM9r+5ylpyHXsRUrVpiQkBCzZMkSs2/fPjNlyhQTFRVlDh8+bIwxZvr06Wbs2LHO5X/44QcTGRlppk6davbt22eWLFliQkJCzHvvveetTfAad/vuxRdfNGvXrjUHDhwwX331lZk+fbqRZFavXu2tTfCKrKwss3v3brN7924jybzwwgtm9+7d5scffzTGsM9Vxd2+Y587Z8KECSYmJsZs3brVHDt2zPnKzc11LsN+VzFP+o797pwZM2aYbdu2mUOHDpk9e/aYmTNnmqCgILN+/XpjjO/sc7YKNcYY88orr5jWrVub0NBQc/nll5e5VG/cuHHmmmuuKbP81q1bTUJCggkNDTXx8fFm0aJFddxi3+FO382fP9+0bdvWhIeHm0aNGpn+/fubDz/80Aut9q7Syz0vfI0bN84Ywz5XFXf7jn3unIr6TJJZunSpcxn2u4p50nfsd+f8/ve/d34/XHzxxWbQoEHOQGOM7+xzPKUbAADYgm3OqQEAAIGNUAMAAGyBUAMAAGyBUAMAAGyBUAMAAGyBUAMAAGyBUAMAAGyBUAMAAGyBUAPA7zkcDr3//vvebgYALyPUAPC67du3Kzg4WMOGDfNo/WPHjmn48OEe17969Wp16dJFYWFh6tKli9auXetxWQC8h1ADwOveeOMNPfzww/rss8+Umprq9vrNmjVTWFiYR3V//vnnuv322zV27Fh9+eWXGjt2rG677Tb9+9//9qg8AN5DqAHgVTk5OXr33Xc1YcIE3XjjjVq2bFmZ+XPnzlVsbKxOnjzpnHbTTTfp6quvVklJiaSyh58KCgo0ceJENW/eXOHh4YqPj1diYmKl9S9YsEDXX3+9ZsyYoU6dOmnGjBkaNGiQFixYUNubCsBihBoAXrVy5Up17NhRHTt21J133qmlS5fq/Ofszpo1S/Hx8br33nslSa+99pq2bdumt99+W0FB5T/CXnrpJX3wwQd69913tX//fv3tb39TfHx8pfV//vnnGjJkSJlpQ4cO1fbt22tnAwHUmXrebgCAwLZkyRLdeeedkqRhw4YpOztbmzZt0uDBgyVJwcHB+tvf/qaePXtq+vTpevnll/X666+rdevWFZaXmpqq9u3bq3///nI4HJUuV+r48eNq2rRpmWlNmzbV8ePHa2HrANQlRmoAeM3+/fv1n//8R3fccYckqV69err99tv1xhtvlFnu0ksv1XPPPaf58+dr5MiRGjNmTKVl3n333UpJSVHHjh01adIkrV+/vtp2OByOMu+NMeWmAfB9jNQA8JolS5aoqKhILVq0cE4zxigkJESnT59Wo0aNnNO3bdum4OBgHT58WEVFRapXr+KPr8svv1yHDh3Sxx9/rI0bN+q2227T4MGD9d5771W4fLNmzcqNyqSnp5cbvQHg+xipAeAVRUVFeuutt/T8888rJSXF+fryyy/VunVr/f3vf3cuu3LlSq1Zs0Zbt25VWlqannrqqSrLbtCggW6//XYtXrxYK1eu1OrVq3Xq1KkKl+3bt682bNhQZtr69evVr1+/mm8kgDrFSA0Ar/jHP/6h06dP65577lFMTEyZeb/97W+1ZMkSTZw4UT/99JMmTJig+fPnq3///lq2bJluuOEGDR8+XFdeeWW5cl988UU1b95cPXv2VFBQkFatWqVmzZqpYcOGFbZj8uTJuvrqqzV//nzdfPPNWrdunTZu3KjPPvvMis0GYCFGagB4xZIlSzR48OBygUaS/ud//kcpKSlKTk7W3XffrSuuuEITJ06UJF1//fWaOHGi7rzzTmVnZ5dbt379+po/f7569+6tPn366PDhw/roo48qvFJKkvr166cVK1Zo6dKluuyyy7Rs2TKtXLlSv/nNb2p3gwFYzmHOv3YSAADATzFSAwAAbIFQAwAAbIFQAwAAbIFQAwAAbIFQAwAAbIFQAwAAbIFQAwAAbIFQAwAAbIFQAwAAbIFQAwAAbIFQAwAAbIFQAwAAbOH/A3GHqtWIoTToAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist1.plot_ratio(hist2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "932f8f9b-b7d9-4a6e-9f34-f2f90aeb86a3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}