binder/applying_corrections.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applying corrections to columnar data\n",
    "\n",
    "Here we will show how to apply corrections to columnar data using:\n",
    "\n",
    "- the `coffea.lookup_tools` package, which is designed to read in ROOT histograms and a variety of data file formats popular within CMS into a standardized lookup table format;\n",
    "- CMS-specific extensions to the above, for jet corrections (`coffea.jetmet_tools`) and b-tagging efficiencies/uncertainties (`coffea.btag_tools`);\n",
    "- the [correctionlib](https://cms-nanoaod.github.io/correctionlib/) package, which provides a experiment-agnostic serializable data format for common correction functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Test data**:\n",
    "We'll use NanoEvents to construct some test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/saransh/Code/HEP/coffea/src/coffea/nanoevents/methods/candidate.py:11: FutureWarning: In version 2024.7.0 (target date: 2024-06-30 11:59:59-05:00), this will be an error.\n",
      "To raise these warnings as errors (and get stack traces to find out where they're called), run\n",
      "    import warnings\n",
      "    warnings.filterwarnings(\"error\", module=\"coffea.*\")\n",
      "after the first `import coffea` or use `@pytest.mark.filterwarnings(\"error:::coffea.*\")` in pytest.\n",
      "Issue: coffea.nanoevents.methods.vector will be removed and replaced with scikit-hep vector. Nanoevents schemas internal to coffea will be migrated. Otherwise please consider using that package!.\n",
      "  from coffea.nanoevents.methods import vector\n"
     ]
    }
   ],
   "source": [
    "import awkward as ak\n",
    "from coffea.nanoevents import NanoEventsFactory, NanoAODSchema\n",
    "\n",
    "NanoAODSchema.warn_missing_crossrefs = False\n",
    "\n",
    "fname = \"https://raw.githubusercontent.com/CoffeaTeam/coffea/master/tests/samples/nano_dy.root\"\n",
    "events = NanoEventsFactory.from_root(\n",
    "    {fname: \"Events\"},\n",
    "    schemaclass=NanoAODSchema,\n",
    "    metadata={\"dataset\": \"DYJets\"},\n",
    ").events()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coffea lookup_tools\n",
    "\n",
    "The entrypoint for `coffea.lookup_tools` is the [extractor class](https://coffeateam.github.io/coffea/api/coffea.lookup_tools.extractor.html#coffea.lookup_tools.extractor)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from coffea.lookup_tools import extractor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "~/Code/HEP/coffea/binder/data ~/Code/HEP/coffea/binder\n",
      "~/Code/HEP/coffea/binder\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "# download some sample correction sources\n",
    "mkdir -p data\n",
    "pushd data\n",
    "PREFIX=https://raw.githubusercontent.com/CoffeaTeam/coffea/master/tests/samples\n",
    "curl -Os $PREFIX/testSF2d.histo.root\n",
    "curl -Os $PREFIX/Fall17_17Nov2017_V32_MC_L2Relative_AK4PFPuppi.jec.txt\n",
    "curl -Os $PREFIX/Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi.junc.txt\n",
    "curl -Os $PREFIX/DeepCSV_102XSF_V1.btag.csv.gz\n",
    "popd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Opening a root file and using it as a lookup table\n",
    "\n",
    "In [tests/samples](https://github.com/CoffeaTeam/coffea/tree/master/tests/samples), there is an example file with a `TH2F` histogram named `scalefactors_Tight_Electron`. The following code reads that histogram into an [evaluator](https://coffeateam.github.io/coffea/api/coffea.lookup_tools.evaluator.html#coffea.lookup_tools.evaluator) instance, under the key `testSF2d` and applies it to some electrons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "available evaluator keys:\n",
      "\t testSF2d\n",
      "testSF2d: <coffea.lookup_tools.dense_lookup.dense_lookup object at 0x313217210> 2 dimensional histogram with axes:\n",
      "\t1: [-2.5   -2.    -1.566 -1.444 -0.8    0.     0.8    1.444  1.566  2.\n",
      "  2.5  ]\n",
      "\t2: [ 10.  20.  35.  50.  90. 150. 500.]\n",
      "\n",
      "type of testSF2d: <class 'coffea.lookup_tools.dense_lookup.dense_lookup'>\n"
     ]
    }
   ],
   "source": [
    "ext = extractor()\n",
    "# several histograms can be imported at once using wildcards (*)\n",
    "ext.add_weight_sets([\"testSF2d scalefactors_Tight_Electron data/testSF2d.histo.root\"])\n",
    "ext.finalize()\n",
    "\n",
    "evaluator = ext.make_evaluator()\n",
    "\n",
    "print(\"available evaluator keys:\")\n",
    "for key in evaluator.keys():\n",
    "    print(\"\\t\", key)\n",
    "print(\"testSF2d:\", evaluator['testSF2d'])\n",
    "print(\"type of testSF2d:\", type(evaluator['testSF2d']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Electron eta: [[], [1.83], [-0.293, -0.904], [-2.19, 1.65], [], ..., [], [0.381], [], []]\n",
      "Electron pt: [[], [29.6], [60.1, 51.7], [10.7, 8.6], [], ..., [15.6], [], [7.68], [], []]\n",
      "Scale factor: [[], [0.909], [0.953, 0.972], [0.807, 0.827], [], ..., [], [0.946], [], []]\n"
     ]
    }
   ],
   "source": [
    "print(\"Electron eta:\", events.Electron.eta.compute())\n",
    "print(\"Electron pt:\", events.Electron.pt.compute())\n",
    "print(\"Scale factor:\", evaluator[\"testSF2d\"](events.Electron.eta, events.Electron.pt).compute())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building and using your own correction from a histogram\n",
    "\n",
    "To use a histogram or ratio of histograms to build your own correction, you can use `lookup_tools` to simplify the implementation. Here we create some mock data for two slightly different pt and eta spectra (say, from two different generators) and derive a correction to reweight one sample to the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/homebrew/lib/python3.11/site-packages/mplhep/utils.py:197: RuntimeWarning: All sumw are zero!  Cannot compute meaningful error bars\n",
      "  return np.abs(method_fcn(self.values, variances) - self.values)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x312239950>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGwCAYAAAC3qV8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIFElEQVR4nO3de1xVVf7/8ffhcpAjN0G5TShqeUtN8hZTGiVf0axvTs5MmaVOTlZfrdQ0tRxTKzVlLGv66q9mxKbRbPp+y/pqeU20GjTFyFuSmoaTIk0mCEc5XPbvD8aTJxEOyOGw4fV8PPajvfdaa+91tsT5sNbaa1kMwzAEAABgIj7ergAAAEBNEcAAAADTIYABAACmQwADAABMhwAGAACYDgEMAAAwHQIYAABgOn7eroCnlJeX68SJEwoODpbFYvF2dQAAgBsMw9DZs2cVGxsrH5/Lt7M02gDmxIkTiouL83Y1AABALRw/flxXXXXVZdMbbQATHBwsqeIBhISEeLk2AADAHQUFBYqLi3N+j19Oow1gLnQbhYSEEMAAAGAy1Q3/YBAvAAAwHQIYAABgOgQwAADAdBrtGBgAQNNUVlamkpISb1cDl+Hv7y9fX98rvg4BDACgUTAMQ7m5uTpz5oy3q4JqhIWFKTo6+ormaSOAAQA0CheCl8jISNlsNiYxbYAMw5DdbldeXp4kKSYmptbXIoABAJheWVmZM3iJiIjwdnVQhcDAQElSXl6eIiMja92dxCBeAIDpXRjzYrPZvFwTuOPCv9OVjFUigAEANBp0G5lDXfw7EcAAAADTIYABAACmQwADAEAdS0pK0oQJE7xdjUaNAAYAAC9KT0+XxWKp9/lrZs2apR49etTrPesSAQwAADAdAhgAAK5AUVGRRo4cqaCgIMXExOiPf/yjS/qbb76pXr16KTg4WNHR0br33nudE7kdO3ZMt9xyiySpRYsWslgsGj16tCRp3bp1uummmxQWFqaIiAjdfvvtOnLkiPO6DodD48ePV0xMjJo1a6Y2bdpo3rx5zvQzZ87o97//vVq1aqWQkBDdeuut+vLLLyVJy5cv1+zZs/Xll1/KYrHIYrFo+fLlHnxKdY8AxmTsjlLFT1ur+GlrZXeUers6ANDkTZkyRVu3btX777+vDRs2KD09Xbt373aml5SU6Nlnn9WXX36p1atX69ixY84gJS4uTv/7v/8rScrOztbJkye1ePFiSRWB0aRJk7Rr1y5t3rxZPj4++tWvfqXy8nJJ0ssvv6wPPvhAf//735Wdna0VK1YoPj7eed/f/OY3ysvL00cffaTMzExdf/31GjBggE6fPq27775bTzzxhK699lqdPHlSJ0+e1N13310/D6yOMBMvJFUERl1mrpckHZiTIpuVHw0AqE5hYaH+8pe/6G9/+5sGDBggSXrjjTd01VVXOfM88MADzv127drp5ZdfVu/evVVYWKigoCCFh4dLkiIjIxUWFubMO2zYMJd7LVu2TK1atdKBAwfUtWtX5eTk6JprrtFNN90ki8WiNm3aOPN++umn+vzzz5WXl6eAgABJUmpqqlavXq3/+Z//0dixYxUUFCQ/Pz9FR0fX+XOpD7TAAABQS0eOHJHD4VDfvn2d58LDw9WxY0fncWZmpu644w61bt1awcHBuvnmmyVJOTk5VV770KFDGj58uNq1a6eQkBBn68qFcqNHj1ZWVpY6duyoxx57TBs2bHCW/fLLL1VYWKiIiAgFBQU5t6NHj7p0Q5kZf2YDAOAhRUVFSklJUUpKilasWKFWrVopJydHKSkpcjgcVZa944471KZNG73++uuKjY1VeXm5unbt6ix3/fXX6+jRo/roo4+0adMm/fa3v1VycrL+53/+R4WFhYqJiVF6evol1724lcfMCGAAAKil9u3by9/fXzt27FDr1q0lST/++KO+/vpr3XzzzTp48KB++OEHzZ8/X3FxcZKkXbt2uVzDarVKqliQ8oIffvhB2dnZev3119WvXz9JFd1CPxcSEqK7775bd999t379619r0KBBOn36tK6//nrl5ubKz8/PZVzMz+978T3NhgAGAIBaCgoK0pgxYzRlyhRFREQoMjJSTz/9tHx8KkZotG7dWlarVa+88ooefvhh7du3T88++6zLNdq0aSOLxaI1a9botttuU2BgoFq0aKGIiAi99tpriomJUU5OjqZNm+ZSbtGiRYqJiVFCQoJ8fHz0zjvvKDo6WmFhYUpOTlZiYqKGDh2qBQsWqEOHDjpx4oTWrl2rX/3qV+rVq5fi4+N19OhRZWVl6aqrrlJwcLBzvIwZMAamATAMQ3ZHqdvbBTUpYxiGFz8hADReCxcuVL9+/XTHHXcoOTlZN910k3r27ClJatWqlZYvX6533nlHXbp00fz585WamupS/he/+IVmz56tadOmKSoqSuPHj5ePj49WrVqlzMxMde3aVRMnTtTChQtdygUHB2vBggXq1auXevfurWPHjunDDz+Uj4+PLBaLPvzwQ/Xv31+/+93v1KFDB91zzz369ttvFRUVJalikPCgQYN0yy23qFWrVnrrrbfq54HVEYvRSL/ZCgoKFBoaqvz8fIWEhHi7OlW6+A0gT6nuzSLeQgJgZufPn9fRo0fVtm1bNWvWzNvVQTWq+vdy9/ubb6kmwu6oup/z5y077gr092X5egBAvatxALNt2zYtXLhQmZmZOnnypN577z0NHTrUmX65L7MFCxZoypQpkqT4+Hh9++23Lunz5s1z6d/bs2ePxo0bp507d6pVq1Z69NFH9eSTT9a0uqaza0aybFbfy6bbHaXq9dzmf+cdUE2rSpl6PbdJkpz/dceF67uD1hoAgDfU+JunqKhI1113nR544AHdddddl6SfPHnS5fijjz7SmDFjLpmQZ86cOXrwwQedx8HBwc79goICDRw4UMnJyVq6dKn27t2rBx54QGFhYRo7dmxNq2wqNquv2wGBzepH8AAAaJJq/O03ePBgDR48+LLpP5/R7/3339ctt9yidu3auZy/sCZEZVasWCGHw6Fly5bJarXq2muvVVZWlhYtWnTZAKa4uFjFxcXO44KCAnc/UqMV6O+rA3NS3Mpb25YdAAC8waNvIZ06dUpr167VmDFjLkmbP3++IiIilJCQoIULF6q09KdxFxkZGerfv7/z3XhJSklJUXZ2tn788cdK7zVv3jyFhoY6twvv2zdlFovF2UrjznZB9Xkv38UFAEB98GgA88Ybbyg4OPiSrqbHHntMq1at0pYtW/TQQw9p7ty5LuNbcnNzna95XXDhODc3t9J7TZ8+Xfn5+c7t+PHjdfxpAABAQ+HRARTLli3TiBEjLnlFatKkSc797t27y2q16qGHHtK8efNqPYlOQECAqSbgAQAAteexFphPPvlE2dnZ+v3vf19t3r59+6q0tFTHjh2TVDGO5tSpUy55LhybddVMAIA52B2lip+2VvHT1tZoWgnUL4+1wPzlL39Rz549dd1111WbNysrSz4+PoqMjJQkJSYm6umnn1ZJSYn8/f0lSRs3blTHjh3VokULT1XZFGxWPx2bP8Tb1QAAwKtq3AJTWFiorKwsZWVlSZJzHYWLlwUvKCjQO++8U2nrS0ZGhl566SV9+eWX+uabb7RixQpNnDhR9913nzM4uffee2W1WjVmzBjt379fb7/9thYvXuzS9QQAANy3f/9+DRs2TPHx8bJYLHrppZe8XaUrUuMWmF27dumWW25xHl8IKkaNGqXly5dLklatWiXDMDR8+PBLygcEBGjVqlWaNWuWiouL1bZtW02cONElOAkNDdWGDRs0btw49ezZUy1bttTMmTMb/Rww3kTLDgA0bna7Xe3atdNvfvMbTZw40dvVuWI1boFJSkqSYRiXbBeCF0kaO3as7Ha7QkNDLyl//fXXa/v27Tpz5ozOnTunAwcOaPr06ZcMwO3evbs++eQTnT9/Xv/85z81derUmn86AABM4OzZsxoxYoSaN2+umJgYvfjii0pKStKECRMkVcx1NnnyZP3iF79Q8+bN1bdvX6WnpzvLL1++XGFhYVq/fr06d+6soKAgDRo0yGVy2d69e2vhwoW65557GsVLL0zjCgBolAzD0LkS99eBu+DcRWvH/VBYLLv10jzVzYJe03XiJk2apM8++0wffPCBoqKiNHPmTO3evVs9evSQJI0fP14HDhzQqlWrFBsbq/fee0+DBg3S3r17dc0111R8FrtdqampevPNN+Xj46P77rtPkydP1ooVK9yuh5kQwAAAGqVzJWXqMnP9FV2j34L0WpWryTpxZ8+e1RtvvKGVK1dqwIABkqS0tDTFxsZKknJycpSWlqacnBznucmTJ2vdunVKS0vT3LlzJUklJSVaunSp2rdvL6ki6JkzZ06t6m8GBDAAAHjRN998o5KSEvXp08d5LjQ0VB07dpQk7d27V2VlZerQoYNLueLiYkVERDiPbTabM3iRpJiYGOXl5Xm49t5DAAMAaJTcWQ/ucl1IF1pePnkySYGVLJ/iThdSXSksLJSvr68yMzPl6+t63aCgIOf+hWlHLrBYLDIMo87q0dAQwAAAGqUL68FVpbL0i4OaiKAAt7uCaqtdu3by9/fXzp071bp1a0lSfn6+vv76a/Xv318JCQkqKytTXl6e+vXr59G6mAkBDAAAXhQcHKxRo0ZpypQpCg8PV2RkpJ555hn5+PjIYrGoQ4cOGjFihEaOHKk//vGPSkhI0Pfff6/Nmzere/fuGjLEvSkwHA6HDhw44Nz/7rvvlJWVpaCgIF199dWe/Ige4dHFHAEAQPUWLVqkxMRE3X777UpOTtaNN96ozp07O9cSTEtL08iRI/XEE0+oY8eOGjp0qEuLjTtOnDihhIQEJSQk6OTJk0pNTVVCQoJbS/40RLTAAADgZcHBwS6vOxcVFWn27NnOCVz9/f01e/ZszZ49u9Lyo0eP1ujRo13ODR061GUMTHx8fKMaE0MAAwCAl33xxRc6ePCg+vTpo/z8fOfrz3feeaeXa9ZwEcAAAHARby2tkpqaquzsbFmtVvXs2VOffPKJWrZsWe/1MAsCGAAAvCwhIUGZmZneroapMIgXAACYDgEMAAAwHQIYAABgOgQwAADAdAhgAACA6RDAAABwMUeRNCu0YnMUebs2uAwCGAAAYDoEMAAANAGvv/66+vXrpxYtWqhFixZKTk7W559/7u1q1RoBDAAATUB6erqGDx+uLVu2KCMjQ3FxcRo4cKC+++47b1etVghgAADwsrNnz2rEiBFq3ry5YmJi9OKLLyopKUkTJkyQJBUXF2vy5Mn6xS9+oebNm6tv375KT093ll++fLnCwsK0fv16de7cWUFBQRo0aJBOnjzpzLNixQr913/9l3r06KFOnTrpz3/+s8rLy7V58+Z6/rR1g6UEAACNk2FIJfaq8zgqSb/4XOG/JGsleay2qq/rb5Mslurr+G+TJk3SZ599pg8++EBRUVGaOXOmdu/erR49ekiSxo8frwMHDmjVqlWKjY3Ve++9p0GDBmnv3r265pprJEl2u12pqal688035ePjo/vuu0+TJ092WeX6Yna7XSUlJQoPD3e7ng0JAQwAoHEqsUtzY6/sGi93r125p05I1uZuZT179qzeeOMNrVy5UgMGDJAkpaWlKTa2ou45OTlKS0tTTk6O89zkyZO1bt06paWlae7cuZKkkpISLV26VO3bt5dUEfRcWNW6MlOnTlVsbKySk5Nr9xm9jAAGAAAv+uabb1RSUqI+ffo4z4WGhqpjx46SpL1796qsrEwdOnRwKVdcXKyIiAjnsc1mcwYvkhQTE6O8vLxK7zl//nytWrVK6enpatasWV1+nHpDAAMAaJz8bRUtIVW5XBfShZaXx/ZU3l3kThdSHSksLJSvr68yMzPl6+vrkhYUFPTTLf39XdIsFosMw7jkeqmpqZo/f742bdqk7t1r2cLUABDAAAAaJ4ul+m6cytIvnrwuqKXbXUG11a5dO/n7+2vnzp1q3bq1JCk/P19ff/21+vfvr4SEBJWVlSkvL0/9+vW7onstWLBAzz//vNavX69evXrVRfW9hgAGAAAvCg4O1qhRozRlyhSFh4crMjJSzzzzjHx8fGSxWNShQweNGDFCI0eO1B//+EclJCTo+++/1+bNm9W9e3cNGTLErfu88MILmjlzplauXKn4+Hjl5uZKqmjFubglxyx4jRoAAC9btGiREhMTdfvttys5OVk33nijOnfu7ByfkpaWppEjR+qJJ55Qx44dNXToUJcWG3csWbJEDodDv/71rxUTE+PcUlNTPfWxPIoWGAAAvCw4ONjldeeioiLNnj1bY8eOlVQxvmX27NmaPXt2peVHjx6t0aNHu5wbOnSoyxiYY8eO1Xm9vYkABgAAL/viiy908OBB9enTR/n5+c7Xn++8804v16zhIoABAOBi1ubSrPx6v21qaqqys7NltVrVs2dPffLJJ2rZsmW918MsCGAAAPCyhIQEZWZmersapsIgXgAAYDoEMAAAwHQIYOBRdkep4qetVfy0tbI7Sr1dHQBAI0EAgwqOImlWaMV28SyUAAA0QAQwAADAdGocwGzbtk133HGHYmNjZbFYtHr1apf00aNHy2KxuGyDBg1yyXP69GmNGDFCISEhCgsL05gxY1RYWOiSZ8+ePerXr5+aNWumuLg4LViwoOafDgCAGrKX2NXtjW7q9kY32UsqWewRDUKNA5iioiJdd911evXVVy+bZ9CgQTp58qRze+utt1zSR4wYof3792vjxo1as2aNtm3b5pxtUJIKCgo0cOBAtWnTRpmZmVq4cKFmzZql1157rabVBQAAjVCN54EZPHiwBg8eXGWegIAARUdHV5r21Vdfad26ddq5c6dzJcxXXnlFt912m1JTUxUbG6sVK1bI4XBo2bJlslqtuvbaa5WVlaVFixa5BDoAAMA9r7/+uv76179q3759kqSePXtq7ty56tOnT53eJykpST169NBLL71Up9f9OY+MgUlPT1dkZKQ6duyoRx55RD/88IMzLSMjQ2FhYS7LeCcnJ8vHx0c7duxw5unfv7+sVqszT0pKirKzs/Xjjz9Wes/i4mIVFBS4bPA8u6NMdkdpldtPeavOd/F28fodAIArl56eruHDh2vLli3KyMhQXFycBg4cqO++++6yZRwORz3WsGbqfCbeQYMG6a677lLbtm115MgRPfXUUxo8eLAyMjLk6+ur3NxcRUZGulbCz0/h4eHOpb1zc3PVtm1blzxRUVHOtBYtWlxy33nz5l12kasmyzAkd/tvHfbK9yvNW6pAnZck9XpuoySLW7fo9dxm9+oi6cCcFNmsTBQNoGk4e/asHn74Ya1evVohISF68skn9f777ztbMoqLi/X000/rrbfe0pkzZ9S1a1e98MILSkpKkiQtX75cEyZM0Ntvv60JEybo+PHjuummm5SWlqaYmBhJclksUpL+/Oc/63//93+1efNmjRw5UpIUHx+vMWPG6NChQ1q9erXuuusuFRYWKjo6Wn/6058kSRMmTNDixYv11VdfqVOnTnI4HGrRooXef/99/e1vf9PWrVu1detWLV68WJJ09OhRxcfH1/kzq/NviHvuuce5361bN3Xv3l3t27dXenq6BgwYUNe3c5o+fbomTZrkPC4oKFBcXJzH7mcKJXZpbmzNy6VeXWWyTdJXFSu8q/P5ZTqnZjW/BwB4mGEYOld6rso8laVffO70+dOV5gn0C6zyuoF+gbJY3PvjTpImTZqkzz77TB988IGioqI0c+ZM7d69Wz169JAkjR8/XgcOHNCqVasUGxur9957T4MGDdLevXt1zTXXSJLsdrtSU1P15ptvysfHR/fdd58mT558SeBygd1uV0lJicLDw13Op6amaubMmXrmmWckSR9++KH+3//7f870rVu3qmXLlkpPT1enTp20c+dOlZSU6Je//KV69+6tr7/+Wl27dnUuSNmqVSu3n0NNePxP3Hbt2qlly5Y6fPiwBgwYoOjoaOXl5bnkKS0t1enTp53jZqKjo3Xq1CmXPBeOLze2JiAgQAEBAR74BKhK5ozkioXPLsPuKHW2vOyaMaDKVhW7o0y9nttU53UE0DSdKz2nviv7XtE1Br9b9ZjPy9lx7w7Z/G1u5T179qzeeOMNrVy50vmHflpammJjK/4AzcnJUVpamnJycpznJk+erHXr1iktLU1z586VJJWUlGjp0qVq3769pIqg50IQUZmpU6cqNjZWycnJLudvvfVWPfHEE87jpKQkPf744/r+++/l5+enAwcO6A9/+IPS09P18MMPKz09Xb1795bNVvF5rVarbDbbZb+v64rHA5h//vOf+uGHH5xNWImJiTpz5owyMzPVs2dPSdLHH3+s8vJy9e3b15nn6aefVklJifz9/SVJGzduVMeOHSvtPoIbJh+WrFX8z+Sw/9TyUoO8Nquf5GZXj83qR7cQAPzMN998o5KSEpfBtKGhoerYsaMkae/evSorK1OHDh1cyhUXFysiIsJ5bLPZnMGLJMXExFzSYHDB/PnztWrVKqWnp6tZM9dW9IvHqEpS165dFR4erq1bt8pqtSohIUG33367823krVu3Oruy6lONv00KCwt1+PBh5/HRo0eVlZWl8PBwhYeHa/bs2Ro2bJiio6N15MgRPfnkk7r66quVkpIiSercubMGDRqkBx98UEuXLlVJSYnGjx+ve+65xxlZ3nvvvZo9e7bGjBmjqVOnat++fVq8eLFefPHFOvrYTZDVVmVLSa3zAkADFegXqB337qgyz+W6kC60vHx010eVdhe504VUVwoLC+Xr66vMzEz5+vq6pAUFBTn3L/zBf4HFYqn0hYjU1FTNnz9fmzZtUvfu3S9Jb97c9fe/xWJR//79lZ6eroCAACUlJal79+4qLi7Wvn379I9//EOTJ0++ko9YKzUOYHbt2qVbbrnFeXxh3MmoUaO0ZMkS7dmzR2+88YbOnDmj2NhYDRw4UM8++6xL986KFSs0fvx4DRgwQD4+Pho2bJhefvllZ3poaKg2bNigcePGqWfPnmrZsqVmzpzJK9QAALdZLJZqu3EqS7948rrwZuFudwXVVrt27eTv76+dO3eqdevWkqT8/Hx9/fXX6t+/vxISElRWVqa8vDz169fviu61YMECPf/881q/fv0lLS1Vufnmm/X6668rICBAzz//vHx8fNS/f38tXLhQxcXFuvHGG515rVarysrKrqie7qhxAJOUlFTlK67r16+v9hrh4eFauXJllXm6d++uTz75pKbVAwDAVIKDgzVq1ChNmTJF4eHhioyM1DPPPCMfHx9ZLBZ16NBBI0aM0MiRI/XHP/5RCQkJ+v7777V582Z1795dQ4YMces+L7zwgmbOnKmVK1cqPj7e+eZvUFCQS0tOZZKSkjRx4kRZrVbddNNNznOTJ09W7969XVpt4uPjtWPHDh07dkxBQUEKDw+Xj0/dz9rCWkgAAHjZokWLlJiYqNtvv13Jycm68cYb1blzZ+f4lLS0NI0cOVJPPPGEOnbsqKFDh7q02LhjyZIlcjgc+vWvf62YmBjnlpqaWm3Zbt26KSwsTD169HAGO0lJSSorK7tk/MvkyZPl6+urLl26qFWrVsrJyXH/QdQAIyrhUTarn47Nd++vAwBoqoKDg11edy4qKtLs2bOdQyf8/f01e/bsy853Nnr0aI0ePdrl3NChQ116TI4dO1ZtPS6Xx8fHR6dPn3Y516NHj0p7ZDp06KCMjIxq73WlCGBQwdpcmpXv7VoAQJP0xRdf6ODBg+rTp4/y8/Odrz/feeedXq5Zw0UAAwDARWz+Nu0dtbfe75uamqrs7GxZrVb17NlTn3zyiVq2bFnv9TALAhgAALwsISFBmZmZ3q6GqTCIFwAAmA4BDACg0WAle3Ooi38nAhgAgOldmIXWbrdXkxMNwYV/p5/PHlwTjIEBAJier6+vwsLCnGv/2Gy2Gq0GjfphGIbsdrvy8vIUFhZ2ydIINUEAAwBoFC6sfny5BQzRcISFhV3xatUEMACARsFisSgmJkaRkZEqKSnxdnVwGf7+/lfU8nIBAQwAoFHx9fWtky9INGwM4gUAAKZDAAMAAEyHAAae5SiSZoVWbI4ib9cGANBIEMAAAADTIYABAACmQwADAABMhwAGAACYDvPA4Mo4qll35OL0avOWKlDnK/ZZkA0AUAUCGFyZ1KvrLK9N0lfNKvbtJTlSQGjt6wUAaNToQgIAAKZDCwxqzt8mPXXCvbwO+08tL5MPS1bbZbPaiwpkW9ypDioIAGjsCGBQcxaLZG1e83JWW9XlHKW1rxMAoEmhCwkAAJgOAQwAADAdAhgAAGA6BDBmw+KIAAAQwAAAAPPhLSR4lrW5NCvf27UAADQytMAAAADTIYABAACmQxdSQ2AYPy1i6ChSlf8sNVocsZp0AABMigCmISix66tmD1Tsp9agXE0WUgQAoBGhCwkAAJgOLTANjP3xg7I1D7l8hhosjujC3818AACYAAFMQ+NfzYKHF6tucUQAABopupAAAIDp1DiA2bZtm+644w7FxsbKYrFo9erVzrSSkhJNnTpV3bp1U/PmzRUbG6uRI0fqxIkTLteIj4+XxWJx2ebPn++SZ8+ePerXr5+aNWumuLg4LViwoHafEAAANDo1DmCKiop03XXX6dVXX70kzW63a/fu3frDH/6g3bt3691331V2drb+8z//85K8c+bM0cmTJ53bo48+6kwrKCjQwIED1aZNG2VmZmrhwoWaNWuWXnvttZpWFwAANEI1HgMzePBgDR48uNK00NBQbdy40eXcn/70J/Xp00c5OTlq3bq183xwcLCio6Mrvc6KFSvkcDi0bNkyWa1WXXvttcrKytKiRYs0duzYSssUFxeruLjYeVxQUFDTjwYAAEzC42Ng8vPzZbFYFBYW5nJ+/vz5ioiIUEJCghYuXKjS0lJnWkZGhvr37y+r1eo8l5KSouzsbP3444+V3mfevHkKDQ11bnFxcR75PF53YW2hWfkM4AUANFkeDWDOnz+vqVOnavjw4QoJ+enV4Mcee0yrVq3Sli1b9NBDD2nu3Ll68sknnem5ubmKiopyudaF49zc3ErvNX36dOXn5zu348ePe+ATAQCAhsBjr1GXlJTot7/9rQzD0JIlS1zSJk2a5Nzv3r27rFarHnroIc2bN08BAQG1ul9AQECtywIAAHPxSAvMheDl22+/1caNG11aXyrTt29flZaW6tixY5Kk6OhonTp1yiXPhePLjZtB02J3lCp+2lrFT1sru6O0+gIAgEalzgOYC8HLoUOHtGnTJkVERFRbJisrSz4+PoqMjJQkJSYmatu2bSopKXHm2bhxozp27KgWLVrUdZUBAIDJ1LgLqbCwUIcPH3YeHz16VFlZWQoPD1dMTIx+/etfa/fu3VqzZo3KysqcY1bCw8NltVqVkZGhHTt26JZbblFwcLAyMjI0ceJE3Xfffc7g5N5779Xs2bM1ZswYTZ06Vfv27dPixYv14osv1tHHBgAAZlbjAGbXrl265ZZbnMcXxrOMGjVKs2bN0gcffCBJ6tGjh0u5LVu2KCkpSQEBAVq1apVmzZql4uJitW3bVhMnTnQZFxMaGqoNGzZo3Lhx6tmzp1q2bKmZM2de9hVqAADQtNQ4gElKSpJhGJdNrypNkq6//npt37692vt0795dn3zySU2rBwAAmgAWc0SDZHeUSVUMzr144G5NBvEG+vvKYrFcUd0AAN5HAIMGqd+CLTqnZm7l7fXcZreve2BOimxWfuwBwOxYjRqmFKjzOtbsXh1rdq8Cdd7b1QEA1DP+FEWDEejv69zPnJFc5VIJ9sJ86eWK/U+eTJItKPTyeR1l6vXcpjqrJwDA+whg0GBcPDbFZvWTqurquSjNZvWjWwgAmhi6kAAAgOkQwAAAANOh3R0Nk8NeTfo5131HFT/KjtKfBvpWM08RAMAcCGDQMKVeXWWy7eL9lztWm/erf7+RbS/JkQIuP+AXAGAOdCEBAADToQUGDYe/TXrqhHt5HfafWmkmH5aststmtRcVyLa4Ux1UEADQUBDAoOGwWKqc++WyrLaqy9VgqQEAgDnQhQQAAEyHAAYAAJgOAQwAADAdAhgAAGA6DOKFOVmbS7PyvV0LAICX0AIDAABMhwAGAACYDgEMAAAwHQIYAABgOgQwAADAdAhgAACA6RDAAAAA0yGAAQAApkMAAwAATIcABgAAmA4BDAAAMB0CGAAAYDoEMAAAwHQIYAAAgOkQwAAAANMhgAEAAKZDAAMAAEyHAAYAAJgOAQwAADAdAhgAAGA6NQ5gtm3bpjvuuEOxsbGyWCxavXq1S7phGJo5c6ZiYmIUGBio5ORkHTp0yCXP6dOnNWLECIWEhCgsLExjxoxRYWGhS549e/aoX79+atasmeLi4rRgwYKafzoAANAo1TiAKSoq0nXXXadXX3210vQFCxbo5Zdf1tKlS7Vjxw41b95cKSkpOn/+vDPPiBEjtH//fm3cuFFr1qzRtm3bNHbsWGd6QUGBBg4cqDZt2igzM1MLFy7UrFmz9Nprr9XiIwIAgMbGr6YFBg8erMGDB1eaZhiGXnrpJc2YMUN33nmnJOmvf/2roqKitHr1at1zzz366quvtG7dOu3cuVO9evWSJL3yyiu67bbblJqaqtjYWK1YsUIOh0PLli2T1WrVtddeq6ysLC1atMgl0AEAAE1TnY6BOXr0qHJzc5WcnOw8Fxoaqr59+yojI0OSlJGRobCwMGfwIknJycny8fHRjh07nHn69+8vq9XqzJOSkqLs7Gz9+OOPld67uLhYBQUFLhsAAGic6jSAyc3NlSRFRUW5nI+KinKm5ebmKjIy0iXdz89P4eHhLnkqu8bF9/i5efPmKTQ01LnFxcVd+QcCAAANUqN5C2n69OnKz893bsePH/d2lQAAgIfUaQATHR0tSTp16pTL+VOnTjnToqOjlZeX55JeWlqq06dPu+Sp7BoX3+PnAgICFBIS4rIBAIDGqU4DmLZt2yo6OlqbN292nisoKNCOHTuUmJgoSUpMTNSZM2eUmZnpzPPxxx+rvLxcffv2debZtm2bSkpKnHk2btyojh07qkWLFnVZZQAAYEI1DmAKCwuVlZWlrKwsSRUDd7OyspSTkyOLxaIJEyboueee0wcffKC9e/dq5MiRio2N1dChQyVJnTt31qBBg/Tggw/q888/12effabx48frnnvuUWxsrCTp3nvvldVq1ZgxY7R//369/fbbWrx4sSZNmlRnHxwAAJhXjV+j3rVrl2655Rbn8YWgYtSoUVq+fLmefPJJFRUVaezYsTpz5oxuuukmrVu3Ts2aNXOWWbFihcaPH68BAwbIx8dHw4YN08svv+xMDw0N1YYNGzRu3Dj17NlTLVu21MyZM3mFGp7nKJLmVgTSeuqEZG3u3foAACplMQzD8HYlPKGgoEChoaHKz89v8ONh7IX5sqW2rtifnCNbUKiXa9S41Oj5EsAAgFe5+/3daN5CAgAATQcBDAAAMB0CGAAAYDo1HsQLmFqJXXJU8WPvsFe+Xx1/m2Sx1L5eAIAaIYBBk2Jb3Mn9zKlXu5+XAb8AUK/oQvIQu6NU8dPWKn7aWtkdpd6uDgAAjQotMGj8/G3qfH6ZJClzRrJs1mq6kC60vEw+LFlt7uUFANQrAhg0fhaLzunfEylam0tVBTAXs9roFgKABoouJAAAYDq0wNSAYRg6V1LmVt6Lx71UNwbG7ihTFR0VAADgZwhgauBcSZm6zFxf43K9nttcZXqgzuurZlVmQX2xNpdm5Xu7FgCAatCFBAAATIcWmFraNSNZNqvvZdPtjlJny8uuGQOqefOlSEqt2A30v/w1AQBABQKYWrJZfasOSlzy+lWT96c0C7O5AgBQLbqQAACA6RDAAAAA0yGAAQAApsMYGA+xWf10bP4Qb1cDNWR3lDpflT8wJ8XtcU4AgPpFCwwAADAdAhgAAGA6BDAAAMB06OBHk2J3VL2WVU3WsJKj1LmGlWEYYgYfAKg/BDBoUno9t6kGed1fw+pcSZlsAVdSMwBATdCFBAAATIcWGDR6gf6+OjAnxa28NVnDyl5YIL3sZiUcRdLc2Ir9p05UrHoNAKg1Ahg0ehaLpVbzuVS7hlUVi3kCADyLLiQAAGA6BDAAAMB06EICLlLrJSBK7JKjiv+dHPbK96vjb5MsvKANAD9HAAPUAdviTu5nTr3a/bwM+AWAStGFBAAATIcWGKC2/G3qfH6ZJClzRnLVbyw57D+1vEw+LFlt7uUFAFSKAAaoLYtF5/TvqXitzSV3X9W22ugWAoArRBcSAAAwHQIYAABgOnQhAfXB2lyale/tWgBAo0ELDAAAMB0CGAAAYDp1HsDEx8fLYrFcso0bN06SlJSUdEnaww8/7HKNnJwcDRkyRDabTZGRkZoyZYpKS0vruqoAAMCk6nwMzM6dO1VWVuY83rdvn/7jP/5Dv/nNb5znHnzwQc2ZM8d5bLP9NCdGWVmZhgwZoujoaP3jH//QyZMnNXLkSPn7+2vu3Ll1XV0AAGBCdR7AtGrVyuV4/vz5at++vW6++WbnOZvNpujo6ErLb9iwQQcOHNCmTZsUFRWlHj166Nlnn9XUqVM1a9YsWa3WSssVFxeruLjYeVxQUFAHnwYAADREHh0D43A49Le//U0PPPCALBctSLdixQq1bNlSXbt21fTp02W3/7S4XUZGhrp166aoqCjnuZSUFBUUFGj//v2Xvde8efMUGhrq3OLi4jzzoQAAgNd59DXq1atX68yZMxo9erTz3L333qs2bdooNjZWe/bs0dSpU5Wdna13331XkpSbm+sSvEhyHufm5l72XtOnT9ekSZOcxwUFBQQxAAA0Uh4NYP7yl79o8ODBio2NdZ4bO3asc79bt26KiYnRgAEDdOTIEbVv377W9woICFBAQMAV1RcAAJiDx7qQvv32W23atEm///3vq8zXt29fSdLhw4clSdHR0Tp16pRLngvHlxs3AwAAmhaPBTBpaWmKjIzUkCFDqsyXlZUlSYqJiZEkJSYmau/evcrLy3Pm2bhxo0JCQtSlSxdPVRcAAJiIR7qQysvLlZaWplGjRsnP76dbHDlyRCtXrtRtt92miIgI7dmzRxMnTlT//v3VvXt3SdLAgQPVpUsX3X///VqwYIFyc3M1Y8YMjRs3ji4iAAAgyUMBzKZNm5STk6MHHnjA5bzVatWmTZv00ksvqaioSHFxcRo2bJhmzJjhzOPr66s1a9bokUceUWJiopo3b65Ro0a5zBsDAACaNo8EMAMHDpRhGJecj4uL09atW6st36ZNG3344YeeqBoAAGgEWAsJAACYDgEMAAAwHQIYAABgOgQwAADAdAhgAACA6RDAAAAA0yGAAQAApkMAAwAATIcAxlMcRdKs0IrNUeTt2gAA0KgQwAAAANMhgAEAAKbjkbWQGi3DUKDOV+w7ilTl43PYK9+vLi8aJbujVF1mrpckHZiTIpuV//UA4ErwW7QmSuz6qtm/V9hOrUG51Ks9Uh0AAJoqupAAAIDp0AJTS/bHD8rWPOTyGRz2n1peJh+WrDb3LuzvZj40KHZHWTXppZXuV8pRqgs/BYZhyHKFdQOAxogAprb8bZK1uXt5rTXIC1Pq9dymGuTdXGV6oM7rq2YV++dKymQLuJKaAUDjRBcSAAAwHVpggFoK9PfVgTkpbuW1O0qdLS+7Zgyo8i0ke2GB9HKdVBEAGi0CGE+xNpdm5Xu7FvAgi8VSq9ehbVa/qstZfa+gVgDQNNCFBAAATIcABgAAmA5dSEA9sFn9dGz+kLq/sKNImhtbsf/UCd52A9Bk0AIDAABMhwAGAACYDl1IQENWYpccdbRo6MX8bZKFOX4BmBcBDNCA2RZ3cj9zTRYNZbwMAJOjCwkAAJgOLTBAQ+NvU+fzyyRJmTOSq570riaLhl6cFwBMjgAGaGgsFp3Tv1dztDaX3J3tl0VDATQhdCEBAADToQUGMDPW3ALQRNECAwAATIcABgAAmA4BDAAAMB0CGAAAYDoEMAAAwHQIYAAAgOnUeQAza9YsWSwWl61Tp5/Wczl//rzGjRuniIgIBQUFadiwYTp16pTLNXJycjRkyBDZbDZFRkZqypQpKi0treuqAgAAk/LIPDDXXnutNm3a9NNN/H66zcSJE7V27Vq98847Cg0N1fjx43XXXXfps88+kySVlZVpyJAhio6O1j/+8Q+dPHlSI0eOlL+/v+bOneuJ6gIAAJPxSADj5+en6OjoS87n5+frL3/5i1auXKlbb71VkpSWlqbOnTtr+/btuuGGG7RhwwYdOHBAmzZtUlRUlHr06KFnn31WU6dO1axZs2S1Wj1RZQAAYCIeGQNz6NAhxcbGql27dhoxYoRycnIkSZmZmSopKVFycrIzb6dOndS6dWtlZGRIkjIyMtStWzdFRUU586SkpKigoED79++/7D2Li4tVUFDgsgEAgMapzgOYvn37avny5Vq3bp2WLFmio0ePql+/fjp79qxyc3NltVoVFhbmUiYqKkq5ubmSpNzcXJfg5UL6hbTLmTdvnkJDQ51bXFxc3X4wAADQYNR5F9LgwYOd+927d1ffvn3Vpk0b/f3vf1dgYGBd385p+vTpmjRpkvO4oKCAIAYAgEbK469Rh4WFqUOHDjp8+LCio6PlcDh05swZlzynTp1yjpmJjo6+5K2kC8eVjau5ICAgQCEhIS4bAABonDwewBQWFurIkSOKiYlRz5495e/vr82bNzvTs7OzlZOTo8TERElSYmKi9u7dq7y8PGeejRs3KiQkRF26dPF0dQEAgAnUeRfS5MmTdccdd6hNmzY6ceKEnnnmGfn6+mr48OEKDQ3VmDFjNGnSJIWHhyskJESPPvqoEhMTdcMNN0iSBg4cqC5duuj+++/XggULlJubqxkzZmjcuHEKCAio6+oCpmZ3lKrLzPWSpANzUmSzeuTFQgBocOr8t90///lPDR8+XD/88INatWqlm266Sdu3b1erVq0kSS+++KJ8fHw0bNgwFRcXKyUlRf/93//tLO/r66s1a9bokUceUWJiopo3b65Ro0Zpzpw5dV1VAABgUnUewKxatarK9GbNmunVV1/Vq6++etk8bdq00YcffljXVQMAAI0E7c1AA2Z3lFWTXlrpfqUcpbL9e9cwDFmusG4A4E0EMEAD1uu5TdVncubdXGV6oM7rq2YV++dKymRjSBkAE2M1agAAYDq0wAANTKC/rw7MSXErr91R6mx52TVjQJVvIdkLC6SX3bioo0iaG1ux/9QJydrcrboAQH0igAEaGIvFUqvXoW1Wv6rLWX2voFYA0LAQwABNUYldclzmf3+HvfL96vjbJAtDgwHUDwIYwMRsVj8dmz+k5uUWd3IvY+rV7l+U7iYA9YhBvAAAwHRogQGaCn+bOp9fJknKnJF8+fEyDvtPLS+TD0tWW+X5fp4XAOoRAQzQVFgsOqd/TwRjbS65M1DYaqNbCECDRBcSAAAwHVpgALiyNpdm5Xu7FgBQJVpgAACA6RDAAKgfjiJpVmjF5ijydm0AmBwBDAAAMB3GwACoG9XN2ssMvwDqEAEMgLpRk/lgmOEXwBWiCwkAAJgOLTAAas/fVtFC4g5m+AVQhwhgANSexVK77h1m+AVwhQhgANQPJsgDUIcYAwMAAEyHAAYAAJgOAQwAADAdAhgAAGA6BDAAAMB0CGAAAIDpEMAAcGF3lCp+2lrFT1sru6PU29UBgEoRwAAAANNhIjugCbI7yqpIK610vzqB/r6ysGo0gHpCAAM0Qb2e2+Rmvs1uX/PAnBTZrPxKAVA/6EICAACmw59LQBMR6O+rA3NSqs1nd5Q6W152zRhQZauK3VHmdmsOANQlAhigibBYLDXu4rFZ/egWAtAg8ZsJgAub1U/H5g/xdjUAoEqMgQFgbo4iaVZoxeYo8nZtANQTAhgAAGA6dCEBaNgcdvfTq8t7MX+bxLw1gGnVeQAzb948vfvuuzp48KACAwP1y1/+Ui+88II6duzozJOUlKStW7e6lHvooYe0dOlS53FOTo4eeeQRbdmyRUFBQRo1apTmzZsnPz9iLqBJSb3aM3mfOiFZm9e8PgAahDqPBrZu3apx48apd+/eKi0t1VNPPaWBAwfqwIEDat78p18WDz74oObMmeM8ttlszv2ysjINGTJE0dHR+sc//qGTJ09q5MiR8vf319y5c+u6ygAAwGTqPIBZt26dy/Hy5csVGRmpzMxM9e/f33neZrMpOjq60mts2LBBBw4c0KZNmxQVFaUePXro2Wef1dSpUzVr1ixZrdZLyhQXF6u4uNh5XFBQUEefCEC987dVtJC4w2H/qeVl8mHJanMvLwBT8/gg3vz8fElSeHi4y/kVK1aoZcuW6tq1q6ZPny67/ae+64yMDHXr1k1RUVHOcykpKSooKND+/fsrvc+8efMUGhrq3OLi4jzwaQDUC4ulonvHre2igMVqcz8vAFPz6ICS8vJyTZgwQTfeeKO6du3qPH/vvfeqTZs2io2N1Z49ezR16lRlZ2fr3XfflSTl5ua6BC+SnMe5ubmV3mv69OmaNGmS87igoIAgBmgKrM2lWfnergWAeubRAGbcuHHat2+fPv30U5fzY8eOde5369ZNMTExGjBggI4cOaL27dvX6l4BAQEKCAi4ovoCAABz8FgX0vjx47VmzRpt2bJFV111VZV5+/btK0k6fPiwJCk6OlqnTp1yyXPh+HLjZgAAQNNR5wGMYRgaP3683nvvPX388cdq27ZttWWysrIkSTExMZKkxMRE7d27V3l5ec48GzduVEhIiLp06VLXVQYAACZT511I48aN08qVK/X+++8rODjYOWYlNDRUgYGBOnLkiFauXKnbbrtNERER2rNnjyZOnKj+/fure/fukqSBAweqS5cuuv/++7VgwQLl5uZqxowZGjduHN1EAACg7ltglixZovz8fCUlJSkmJsa5vf3225Ikq9WqTZs2aeDAgerUqZOeeOIJDRs2TP/3f//nvIavr6/WrFkjX19fJSYm6r777tPIkSNd5o0BAABNV523wBiGUWV6XFzcJbPwVqZNmzb68MMP66paALzM7ihVl5nrJUkH5qTIZmVWbQC1x28QAHXC7iirJr200v3qBPr7ysKaRQB+hgAGQJ3o9dymGuTd7HZeWmsAVMbjM/ECgCk5iqRZoRWbo8jbtQHwM/xZA6DWAv19dWBOilt57Y5SZ8vLrhkDqmxVsTvKatSiUysOu/vp1eW9mL+tYikEAB5FAAOg1iwWi9vdOzarn47NH+LhGtVATRZ1rEnep05ULG8AwKPoQgIAAKZDCwyApsPfVtFC4g6H/aeWl8mHq17J+uK8AOoFAQyApsNicb97h1WugQaNLiQAAGA6BDAAAMB0CGAAAIDpEMAAMDW7o1Tx09YqftraGi1RAMDcGMQLoEFrdGssOYqkubEV+8wZA9QaAQyABs10aywxwy9QLwhgAKAuMcMvUC8IYAA0OKZeYwlAvSCAAdDgmG6NJWb4BeodAQwAXClm+AXqHa9RAwAA0yGAAQAApkMXEoAmyZ35ZdwdHHwxr80vAzQxBDAAmqRGOb+Mu4ODL8b8MjApAhgAMAPmlwFcEMAAaDKYXwZoPAhgADQZzC9TSV7ApAhgAKChqo/5ZRhbA5MigAGApoyxNTAp5oEBAACmQwsMANQhU8wvU5OxNTVxcXcTXVONwvHjx7VmzRo9+OCD8vNz72e1NmVqw2IYhuGxq3tRQUGBQkNDlZ+fr5CQkDq5pr0wX7bU1hX7k3NkCwqtk+sCMDe7o1RdZq736D12zUiWzepbZR28Hhg5iqS5sXVzrcuha6pejRo1SuvWrdOhQ4fc/i6tTZmLufv9TQsMAJiA6Sbeg+nt3r1bf/3rX7VkyRK3A5HalKktfmoB4ArVZH6ZmjDd/DJ0TTUahmHoiSeeUOfOnfX73//eY2WuBAEMAFyhmswvUxOmC4xq8tp3bfHWVL1Yu3at0tPTtWbNGrfHsdSmzJUggAGABspTgdHFqht03CA4SuVmO8oV3KOalp3aMGmrTmlpqaZMmaJbb71Vt912m8fKXCkCGABowszRRWUoUMskSZkzkusuqLu4u8kTMxObtFXn9ddfV3Z2tlauXOn2AO/alLlSBDAAgAbOonNqVrFrbS4x6NhjCgoK9Mwzz+j+++9XQkKCx8rUBX4KAKCJ8dTYGk+5eMxOnXZ5GVZpco6kOnydvCYDjhsaf5teeOEFnT17Vs8//7zbxWpTpi406ADm1Vdf1cKFC5Wbm6vrrrtOr7zyivr06ePtagGAqdXH2BpP8VSXV3Xz7LjvovE6Jlsw8/j9n2vRokV64okndNVVV7lX5vjxGpepKw32J/jtt9/WpEmTtHTpUvXt21cvvfSSUlJSlJ2drcjISG9XDwDQiNRVYBSo8/qqWZ1cqt7NeGa2QkJCNHXqVPfLzJhR4zJ1pcHOxNu3b1/17t1bf/rTnyRJ5eXliouL06OPPqpp06ZVW96TM/HaLRade/ygApt7dpIeAEDF/CLnS8vr/Lp2R5n6vbDFjQpIkvXfBw6pyp4mQ4Eqduv+hiGdV8X3SDMVeOWFpUAV69NmE3TwRKl6vlakJUuW6OGHH3ar7O7du9WzZ88alXGHqWfidTgcyszM1PTp053nfHx8lJycrIyMjErLFBcXq7j4px+a/PyKZeULCgrqrF72wgKVFhu6oc0vpFW31tl1AQDeYWvjvXtbJM+/Hl6N85J6GjFquexrXXP11frtb3/r1vemYRh6/PHH1aFDB7fLuOvCtaprX2mQAcy//vUvlZWVKSoqyuV8VFSUDh48WGmZefPmafbs2Zecj4uL80ANv/LANQEA8KbDioiIqHGp2pRxx9mzZxUaevk1BxtkAFMb06dP16RJk5zH5eXlOn36tCIiIur0nfSCggLFxcXp+PHjHl/noSni+XoWz9dzeLaexfP1jNLSUvXp00dHjhxRTk5OlQHDxWVuuOEGxcTE6IMPPqjzeV8Mw9DZs2cVG1v1wqANMoBp2bKlfH19derUKZfzp06dUnR0dKVlAgICFBAQ4HIuLCzMU1VUSEgI/xN5EM/Xs3i+nsOz9Syeb91asmSJvvnmG0lSaGioW892yZIlOnz4sN5++223Ap7acOe6Ph658xWyWq3q2bOnNm/+aUXV8vJybd68WYmJiV6sGQAAjcOFCeiGDx9e4zIjR46s10nrKtMgAxhJmjRpkl5//XW98cYb+uqrr/TII4+oqKhIv/vd77xdNQAATO/CBHQzZsyocZnnnnvOgzVzT4PsQpKku+++W99//71mzpyp3Nxc9ejRQ+vWrbtkYG99CwgI0DPPPHNJdxXqBs/Xs3i+nsOz9Syeb926eAK6du3aufVsvTlpXWUa7DwwAADAM0aOHKn169fr8OHDCg4O9lgZT2qwLTAAAKDu7d69W2+++aaWLFnidiBSmzKeRgsMAABNhGEYuvXWW3Xq1Cnt2bNHfn7Vt2PUpkx9aBi1AAAAHrd27Vqlp6drzZo1bgcitSlTH2iBAQCgCSgpKVH37t0VGxurTZs2uTUBXW3K1JcG+xp1Q/Xqq68qPj5ezZo1U9++ffX55597u0qmM2/ePPXu3VvBwcGKjIzU0KFDlZ2d7ZLn/PnzGjdunCIiIhQUFKRhw4ZdMrEhqjd//nxZLBZNmDDBeY5ne2W+++473XfffYqIiFBgYKC6deumXbt2OdMNw9DMmTMVExOjwMBAJScn69ChQ16ssXmUlZXpD3/4g9q2bavAwEC1b99ezz77rMuaODxf923btk133HGHYmNjZbFY9Oijjyo7O1upqamyWCxuPcvFixfr4MGD2r59u1q0aKExY8aosLDQS5/oZwy4bdWqVYbVajWWLVtm7N+/33jwwQeNsLAw49SpU96umqmkpKQYaWlpxr59+4ysrCzjtttuM1q3bm0UFhY68zz88MNGXFycsXnzZmPXrl3GDTfcYPzyl7/0Yq3N5/PPPzfi4+ON7t27G48//rjzPM+29k6fPm20adPGGD16tLFjxw7jm2++MdavX28cPnzYmWf+/PlGaGiosXr1auPLL780/vM//9No27atce7cOS/W3Byef/55IyIiwlizZo1x9OhR45133jGCgoKMxYsXO/PwfN334YcfGk8//bTx7rvvGpKMkJAQY9SoUc706p5lfn6+YbVajRYtWhjbt283PvnkE+Pqq682hg8f7qVP5IoApgb69OljjBs3znlcVlZmxMbGGvPmzfNircwvLy/PkGRs3brVMAzDOHPmjOHv72+88847zjxfffWVIcnIyMjwVjVN5ezZs8Y111xjbNy40bj55pudAQzP9spMnTrVuOmmmy6bXl5ebkRHRxsLFy50njtz5owREBBgvPXWW/VRRVMbMmSI8cADD7icu+uuu4wRI0YYhsHzvRKSjIceesg4fvy4YRjuPcs9e/YYkow1a9Y483z00UeGxWIxvvvuu/r9AJWgC8lNDodDmZmZSk5Odp7z8fFRcnKyMjIyvFgz88vPz5ckhYeHS5IyMzNVUlLi8qw7deqk1q1b86zdNG7cOA0ZMsTlGUo82yv1wQcfqFevXvrNb36jyMhIJSQk6PXXX3emHz16VLm5uS7PNzQ0VH379uX5uuGXv/ylNm/erK+//lqS9OWXX+rTTz/V4MGDJfF8r9SgQYOcE9C58yx37typsLAwDRkyxJknOTlZPj4+2rFjR/1WvhINZzhxA/evf/1LZWVll8wEHBUVpYMHD3qpVuZXXl6uCRMm6MYbb1TXrl0lSbm5ubJarZcsxhkVFaXc3Fwv1NJcVq1apd27d2vnzp2XpPFsr8w333yjJUuWaNKkSXrqqae0c+dOPfbYY7JarRo1apTzGVb2e4LnW71p06apoKBAnTp1kq+vr8rKyvT8889rxIgRksTzrUPuPMvc3FxFRka6pPv5+Sk8PLxBPG8CGHjVuHHjtG/fPn366aferkqjcPz4cT3++OPauHGjmjVr5u3qNDrl5eXq1auX5s6dK0lKSEjQvn37tHTpUo0aNcrLtTO/v//971qxYoVWrlypa6+9VllZWZowYYJiY2N5vrgEXUhuatmypXx9fS95W+PUqVOKjo72Uq3Mbfz48VqzZo22bNnisq5GdHS0HA6Hzpw545KfZ129zMxM5eXl6frrr5efn5/8/Py0detWvfzyy/Lz81NUVBTP9grExMSoS5cuLuc6d+6snJwcSXI+Q35P1M6UKVM0bdo03XPPPerWrZvuv/9+TZw4UfPmzZPE861L7jzL6Oho5eXluaSXlpbq9OnTDeJ5E8C4yWq1qmfPntq8ebPzXHl5uTZv3qzExEQv1sx8DMPQ+PHj9d577+njjz9W27ZtXdJ79uwpf39/l2ednZ2tnJwcnnU1BgwYoL179yorK8u59erVSyNGjHDu82xr78Ybb7zklf+vv/5abdq0kSS1bdtW0dHRLs+3oKBAO3bs4Pm6wW63y8fH9WvJ19dX5eXlkni+dcmdZ5mYmKgzZ84oMzPTmefjjz9WeXm5+vbtW+91voS3RxGbyapVq4yAgABj+fLlxoEDB4yxY8caYWFhRm5urrerZiqPPPKIERoaaqSnpxsnT550bna73Znn4YcfNlq3bm18/PHHxq5du4zExEQjMTHRi7U2r4vfQjIMnu2V+Pzzzw0/Pz/j+eefNw4dOmSsWLHCsNlsxt/+9jdnnvnz5xthYWHG+++/b+zZs8e48847ec3XTaNGjTJ+8YtfOF+jfvfdd42WLVsaTz75pDMPz9d9Z8+eNb744gvjiy++MCQZixYtMr744gvj22+/NQzDvWc5aNAgIyEhwdixY4fx6aefGtdccw2vUZvVK6+8YrRu3dqwWq1Gnz59jO3bt3u7SqYjqdItLS3NmefcuXPGf/3XfxktWrQwbDab8atf/co4efKk9yptYj8PYHi2V+b//u//jK5duxoBAQFGp06djNdee80lvby83PjDH/5gREVFGQEBAcaAAQOM7OxsL9XWXAoKCozHH3/caN26tdGsWTOjXbt2xtNPP20UFxc78/B83bdly5ZKf9demAvGnWf5ww8/GMOHDzeCgoKMkJAQ43e/+51x9uxZL3yaS7GUAAAAMB3GwAAAANMhgAEAAKZDAAMAAEyHAAYAAJgOAQwAADAdAhgAAGA6BDAAAMB0CGAAAIDpEMAAAADTIYABYGrx8fF66aWXvF0NAPWMAAYAAJgOayEBaNCSkpLUtWtXSdKbb74pf39/PfLII5ozZ45uueUWbd261SU/v9KApoEWGAAN3htvvCE/Pz99/vnnWrx4sRYtWqQ///nPevfdd3XVVVdpzpw5OnnypE6ePOntqgKoJ37ergAAVCcuLk4vvviiLBaLOnbsqL179+rFF1/Ugw8+KF9fXwUHBys6Otrb1QRQj2iBAdDg3XDDDbJYLM7jxMREHTp0SGVlZV6sFQBvIoABAACmQwADoMHbsWOHy/H27dt1zTXXyNfXV1arlZYYoAkigAHQ4OXk5GjSpEnKzs7WW2+9pVdeeUWPP/64pIp5YLZt26bvvvtO//rXv7xcUwD1hdeoATRoSUlJuvbaa1VeXq6VK1fK19dXjzzyiJ577jlZLBZt375dDz30kLKzs1VcXMxr1EATQQADoEFLSkpSjx49mG0XgAu6kAAAgOkQwAAAANOhCwkAAJgOLTAAAMB0CGAAAIDpEMAAAADTIYABAACmQwADAABMhwAGAACYDgEMAAAwHQIYAABgOv8fZ8Ryk6YuN2QAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import hist\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "dists = (\n",
    "    hist.Hist.new\n",
    "    .StrCat([\"gen1\", \"gen2\", \"gen2rwt\"], name=\"dataset\")\n",
    "    .Reg(20, 0, 100, name=\"pt\")\n",
    "    .Reg(4, -3, 3, name=\"eta\")\n",
    "    .Weight()\n",
    "    .fill(\n",
    "        dataset=\"gen1\",\n",
    "        pt=np.random.exponential(scale=10.0, size=10000) + np.random.exponential(scale=10.0, size=10000),\n",
    "        eta=np.random.normal(scale=1, size=10000)\n",
    "    )\n",
    "    .fill(\n",
    "        dataset=\"gen2\",\n",
    "        pt=np.random.exponential(scale=10.0, size=10000) + np.random.exponential(scale=15.0, size=10000),\n",
    "        eta=np.random.normal(scale=1.1, size=10000)\n",
    "    )\n",
    ")\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "dists[:, :, sum].plot1d(ax=ax)\n",
    "ax.legend(title=\"dataset\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we derive a correction as a function of $p_T$ and $\\eta$ to `gen2` such that it agrees with `gen1`. We'll set it to 1 anywhere we run out of statistics for the correction, to avoid divide by zero issues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<coffea.lookup_tools.dense_lookup.dense_lookup object at 0x319bfba50> 2 dimensional histogram with axes:\n",
      "\t1: [  0.   5.  10.  15.  20.  25.  30.  35.  40.  45.  50.  55.  60.  65.\n",
      "  70.  75.  80.  85.  90.  95. 100.]\n",
      "\t2: [-3.  -1.5  0.   1.5  3. ]\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "ColormeshArtists(pcolormesh=<matplotlib.collections.QuadMesh object at 0x319d24150>, cbar=<matplotlib.colorbar.Colorbar object at 0x319d18090>, text=[])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAG2CAYAAADcEepCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1dElEQVR4nO3de3hU1bnH8d8kIRMoTgBDLmC4eeFSEBBKDGgFiQYOD0r1cChyBBHxSIkFoiIoEPFCvFSEWjSKImpBEKtYEeHQYOBBwy0QixeiXDQ5SAKIEAiQwMw+fyBTpwTcCWvYmfD9PM96yuzZe62V1SG+vO/ae1yWZVkCAADAeRfm9AQAAAAuVARiAAAADiEQAwAAcAiBGAAAgEMIxAAAABxCIAYAAOAQAjEAAACHEIgBAAA4hEAMAADAIQRiAAAADgmZQOzFF1/UlVdeKY/HI4/Ho+TkZH300UdnvWbRokVq06aNoqKi1KFDBy1duvQ8zRYAAOCXhUwgdskll+jJJ59UXl6eNm7cqOuvv14333yzvvjii0rP//TTTzV48GCNGDFCmzdv1oABAzRgwAB9/vnn53nmAAAAlXOF8pd+N2rUSM8884xGjBhx2nuDBg1SWVmZlixZ4j929dVXq1OnTsrKyjqf0wQAAKhUhNMTqA6v16tFixaprKxMycnJlZ6Tm5ur9PT0gGOpqalavHjxWfsuLy9XeXm5/7XP59P+/ft18cUXy+VynfPcAQCoKSzL0qFDh9SkSROFhZkpkh07dkwVFRVG+jpXkZGRioqKcnoaZxVSgdiWLVuUnJysY8eOqX79+nrvvffUrl27Ss8tLi5WXFxcwLG4uDgVFxefdYzMzExNnTrV2JwBAKjpioqKdMkll5xzP8eOHVPL5vVVvMdrYFbnLj4+Xjt37qzRwVhIBWKtW7dWfn6+Dh48qHfeeUfDhg3TqlWrzhiMVcfEiRMDMmkHDx5Us2bN1PSxSQoz/H9k3SaHjfYXbMcrwoPSb90NvwpKv153ULqVNzI4/cavD9K/IIOUyI3avi84HUv6MSkhKP2eiArOYkRvPxqUfiN37AlKvyeaNQ5Kv8digvOXbvfvgvN3w/IF5/PgCgudHT++o+UquvdpXXTRRUb6q6ioUPEer77LayHPRc5uQy895FPzLt+qoqKCQMyUyMhIXXbZZZKkLl26aMOGDZo5c6Zeeuml086Nj49XSUlJwLGSkhLFx8efdQy32y23+/RfJmFRUQqra/b/yPB6x432F2zeiOB8XMLdQfoLEqRATEEKxCIigvRLK0iBWERYsBZYiqgTnM+EFRmcxYiICM5/eCPCgvVhC876RtQJzmcirF5w/m4QiP2L6a039S9yqf5Fzm7n8QXrl59hIXPXZGV8Pl/Afq6fS05OVnZ2dsCxFStWnHFPGQAAMMNr+WpECwUhkxGbOHGi+vbtq2bNmunQoUOaP3++cnJytHz5cknS0KFD1bRpU2VmZkqSxowZo+uuu07PPvus+vXrpwULFmjjxo16+eWXnfwxAACo9Xyy5JOzmUGnx7crZAKxPXv2aOjQodq9e7eio6N15ZVXavny5brhhhskSYWFhQF3fHTv3l3z58/XpEmT9NBDD+nyyy/X4sWL1b59e6d+BAAAgAAhE4i9+uqrZ30/JyfntGMDBw7UwIEDgzQjAABQGZ98crow6PwM7AmZQAwAAIQGr2XJ6/Dz4p0e366Q3qwPAAAQysiIAQAAo9isbx+BGAAAMMonS14CMVsoTQIAADiEjBgAADCK0qR9BGIAAMAo7pq0j9IkAACAQ8iIAQAAo3w/NafnEAoIxAAAgFHeGnDXpNPj20UgBgAAjPJaJ5vTcwgF7BEDAABwCBkxAABgFHvE7CMQAwAARvnkklcux+cQCihNAgAAOISMGAAAMMpnnWxOzyEUEIgBAACjvDWgNOn0+HZRmgQAAHAIGTEAAGAUGTH7CMQAAIBRPssln+XwXZMOj28XpUkAAACHkBEDAABGUZq0j0AMAAAY5VWYvA4X3byOjm4fgRgAADDKqgF7xCz2iAEAAOBsyIgBAACj2CNmH4EYAAAwymuFyWs5vEcsRL7iiNIkAACAQ8iIAQAAo3xyyedwrsen0EiJEYgBAACj2CNmH6VJAAAAh5ARAwAARtWMzfqUJgEAwAXo5B4xh7/0m9IkAAAAzoaMGAAAMMpXA75rkrsmAQDABYk9YvYRiAEAAKN8CuM5YjaxRwwAAMAhZMQAAIBRXsslr+XwA10dHt8uMmIAAMAo70+b9Z1uVbF69Wr1799fTZo0kcvl0uLFi21f+8knnygiIkKdOnWq2kKJQAwAAEBlZWXq2LGjZs2aVaXrDhw4oKFDh6p3797VGpfSJAAAMMpnhcnn8F2TvireNdm3b1/17du3yuPcc889uu222xQeHl6lLNopZMQAAIBRTpckf16aLC0tDWjl5eXGfs7XXntNO3bsUEZGRrX7IBADAAC1VmJioqKjo/0tMzPTSL/ffPONJkyYoL/+9a+KiKh+gZHSJAAAMMon5+9a9P30v0VFRfJ4PP7jbrf7nPv2er267bbbNHXqVF1xxRXn1BeBGAAAMKpmPND15PgejycgEDPh0KFD2rhxozZv3qy0tLST4/l8sixLERER+t///V9df/31tvoiEAMAAKgCj8ejLVu2BBx74YUXtHLlSr3zzjtq2bKl7b4IxAAAgFE147smqzb+4cOHtW3bNv/rnTt3Kj8/X40aNVKzZs00ceJE7dq1S2+88YbCwsLUvn37gOtjY2MVFRV12vFfQiAGAACM8skln5zeI1a18Tdu3KhevXr5X6enp0uShg0bprlz52r37t0qLCw0OkeJQAwAABgWihmxnj17yjrLs8fmzp171usfeeQRPfLII1UaU+LxFQAAAI4hIwYAAIyqznc9BmMOoYBADAAAGOWzXPI5/Rwxh8e3KzTCRQAAgFqIjBgAADDKVwNKk04/UNYuAjEAAGCUzwqTz+G7Jp0e367QmCUAAEAtREYMAAAY5ZVLXocf6Or0+HYRiAEAAKMoTdoXGrMEAACohciIAQAAo7xyvjTodXR0+wjEAACAUZQm7SMQAwAARoXil347JTRmCQAAUAuREQMAAEZZcsnn8B4xi8dXAACACxGlSftCY5YAAAC1EBkxAABglM9yyWc5Wxp0eny7CMQAAIBRXoXJ63DRzenx7QqNWQIAANRCZMQAAIBRlCbtIxADAABG+RQmn8NFN6fHtys0ZgkAAFALkREDAABGeS2XvA6XBp0e3y4CMQAAYBR7xOwjEAMAAEZZVph8Dj/Z3uLJ+gAAADgbMmIAAMAor1zyOvyl206PbxeBGAAAMMpnOb9Hy2c5OrxtlCYBAAAcQkYMAAAY5asBm/WdHt+u0JjlT1avXq3+/furSZMmcrlcWrx48VnPz8nJkcvlOq0VFxefnwkDAHAB8slVI1ooCKlArKysTB07dtSsWbOqdF1BQYF2797tb7GxsUGaIQAAgH0hVZrs27ev+vbtW+XrYmNj1aBBA/MTAgAAp+HJ+vaFVEasujp16qSEhATdcMMN+uSTT5yeDgAAtdqpPWJOt1AQUhmxqkpISFBWVpa6du2q8vJyvfLKK+rZs6fWrVunq666qtJrysvLVV5e7n9dWlp6vqYLAAAuMLU6EGvdurVat27tf929e3dt375dzz33nN58881Kr8nMzNTUqVNPO/5YyiLVuyjc6Px6RO022t8pj5dcH5R+vyltHJR+dyQFJ33s/dEdlH5dJ4Iz3523mP18nXLRN8HpN/pX8UHpV5LCy4PzAKDIQ96g9CtXcD4TX2UkBqXfK+YcC0q/x+sFZx0u/1P5L59UDa7vgvM72NumeVD61af5xrs8YR3Xd8Z7/WmzvtPPEWOzfs3UrVs3bdu27YzvT5w4UQcPHvS3oqKi8zg7AABCn1UD7pi0QiQQq9UZscrk5+crISHhjO+73W653cHJpAAAcCHwWTUgI8ZmffMOHz6s/Px85efnS5J27typ/Px8FRYWSjqZzRo6dKj//BkzZuj999/Xtm3b9Pnnn2vs2LFauXKlRo8e7cT0AQBADVXVZ5W+++67uuGGG9S4cWN5PB4lJydr+fLlVR43pAKxjRs3qnPnzurcubMkKT09XZ07d9aUKVMkSbt37/YHZZJUUVGh++67Tx06dNB1112nzz77TP/4xz/Uu3dvR+YPAMCFwOm7Jatz12RVn1W6evVq3XDDDVq6dKny8vLUq1cv9e/fX5s3b67SuCFVmuzZs6cs68ybeOfOnRvwevz48Ro/fnyQZwUAAH4uFEuTVX1W6YwZMwJeT5s2Te+//74++OADf8LIjpDKiAEAANREPp9Phw4dUqNGjap0XUhlxAAAQM1XE77r8dT4//480GDdlPenP/1Jhw8f1n/9139V6ToyYgAAwKhTpUmnmyQlJiYqOjra3zIzM43/vPPnz9fUqVP19ttvV/n7rMmIAQCAWquoqEgej8f/2nQ2bMGCBbrrrru0aNEipaSkVPl6AjEAAGBUTdqs7/F4AgIxk9566y3deeedWrBggfr161etPgjEAACAUTUpELPr8OHDAd+8c+pZpY0aNVKzZs00ceJE7dq1S2+88Yakk+XIYcOGaebMmUpKSlJxcbEkqW7duoqOjrY9LnvEAADABa+qzyp9+eWXdeLECY0ePVoJCQn+NmbMmCqNS0YMAAAYFYoZsao+qzQnJ6caszodgRgAADDKkhx/fMWZQ6qahUAMAAAYFYoZMaewRwwAAMAhZMQAAIBRZMTsIxADAABGEYjZR2kSAADAIWTEAACAUWTE7CMQAwAARlmWS5bDgZDT49tFaRIAAMAhZMQAAIBRPrkcf6Cr0+PbRSAGAACMYo+YfZQmAQAAHEJGDAAAGMVmffsIxAAAgFGUJu0jEAMAAEaREbOPPWIAAAAOISMGAACMsmpAaTJUMmIEYgAAwChLkmU5P4dQQGkSAADAIWTEAACAUT655OLJ+rYQiAEAAKO4a9I+SpMAAAAOISMGAACM8lkuuXigqy0EYgAAwCjLqgF3TYbIbZOUJgEAABxCRgwAABjFZn37CMQAAIBRBGL2EYgBAACj2KxvH3vEAAAAHEJGDAAAGMVdk/YRiAEAAKNOBmJO7xFzdHjbKE0CAAA4hIwYAAAwirsm7SMQAwAARlk/NafnEAooTQIAADiEjBgAADCK0qR9BGIAAMAsapO2EYgBAACzakBGTE6PbxN7xAAAABxCRgwAABjFk/XtIyMGAACMOrVZ3+lWFatXr1b//v3VpEkTuVwuLV68+BevycnJ0VVXXSW3263LLrtMc+fOrfJaEYgBAIALXllZmTp27KhZs2bZOn/nzp3q16+fevXqpfz8fI0dO1Z33XWXli9fXqVxKU0CAACzLJfzm+WrOH7fvn3Vt29f2+dnZWWpZcuWevbZZyVJbdu21Zo1a/Tcc88pNTXVdj8EYgAAwKiatEestLQ04Ljb7Zbb7T7n/nNzc5WSkhJwLDU1VWPHjq1SP5QmAQBArZWYmKjo6Gh/y8zMNNJvcXGx4uLiAo7FxcWptLRUR48etd0PGTEAAGBWDXqga1FRkTwej/+wiWyYSQRiAADAqJr0FUcejycgEDMlPj5eJSUlAcdKSkrk8XhUt25d2/1QmgQAAKii5ORkZWdnBxxbsWKFkpOTq9QPgRgAADDPcrhV0eHDh5Wfn6/8/HxJJx9PkZ+fr8LCQknSxIkTNXToUP/599xzj3bs2KHx48dr69ateuGFF/T2229r3LhxVRqX0iQAADCqJpUm7dq4caN69erlf52eni5JGjZsmObOnavdu3f7gzJJatmypT788EONGzdOM2fO1CWXXKJXXnmlSo+ukAjEAACAaTVos75dPXv2lHWWZ25U9tT8nj17avPmzVWcWCBKkwAAAA4hIwYAAAxz/dScnkPNRyAGAADMCsHSpFMoTQIAADiEjBgAADCLjJhtBGIAAMAsy3WyOT2HEEBpEgAAwCFkxAAAgFGWdbI5PYdQQCAGAADMYo+YbZQmAQAAHEJGDAAAmMVmfdsIxAAAgFEu62Rzeg6hgEAMAACYxR4x29gjBgAA4BAyYgAAwCz2iNlGIAYAAMyiNGkbpUkAAACHkBEDAABmkRGzjUAMAACYRSBm2zkFYl9++aUKCwtVUVERcPymm246p0kBAABcCKoViO3YsUO/+93vtGXLFrlcLlk/fbOmy3XyDgWv12tuhgAAILRw16Rt1dqsP2bMGLVs2VJ79uxRvXr19MUXX2j16tXq2rWrcnJyDE8RAACEklNP1ne6hYJqZcRyc3O1cuVKxcTEKCwsTGFhYbrmmmuUmZmpP/7xj9q8ebPpeQIAANQ61cqIeb1eXXTRRZKkmJgYff/995Kk5s2bq6CgwNzsKjFr1iy1aNFCUVFRSkpK0vr16896/qJFi9SmTRtFRUWpQ4cOWrp0aVDnBwDABc+qIS0EVCsQa9++vT777DNJUlJSkp5++ml98sknevTRR9WqVSujE/y5hQsXKj09XRkZGdq0aZM6duyo1NRU7dmzp9LzP/30Uw0ePFgjRozQ5s2bNWDAAA0YMECff/550OYIAABgV7UCsUmTJsnn80mSHn30Ue3cuVPXXnutli5dqpkzZxqd4M9Nnz5dI0eO1PDhw9WuXTtlZWWpXr16mjNnTqXnz5w5U3369NEDDzygtm3b6rHHHtNVV12lv/zlL0GbIwAAFzqXnN8fFhpb9au5Ryw1NdX/58suu0xbt27V/v371bBhQ/+dk6ZVVFQoLy9PEydO9B8LCwtTSkqKcnNzK70mNzdX6enpp8198eLFZxynvLxc5eXl/telpaXnNnEAAIAzqFYgduedd2rmzJn+fWKS1KhRI5WVlenee+89Y4bqXOzbt09er1dxcXEBx+Pi4rR169ZKrykuLq70/OLi4jOOk5mZqalTp552/OZfHZLnV+HVmPmZPbrvN0b7O6WorGFQ+t1bVj8o/SZcfDAo/RYdahyUfusVmf0cnOKNCkq3OlEvOP0eahacdZCkunt9Qen3WMPgPMO60Q8Vv3xSNbSZdSgo/Z6IDs6HzQoPUg4iSP/AP9CnTVD6Lelb/ssnVcMljZOM93ni+DHpw/eN98vjK+yrVmny9ddf19GjR087fvToUb3xxhvnPCknTZw4UQcPHvS3oqIip6cEAEBocXqTfght1q/SPw9LS0tlWZYsy9KhQ4cUFfWvf1V5vV4tXbpUsbGxxicpnbw7Mzw8XCUlJQHHS0pKFB8fX+k18fHxVTpfktxut9xu97lPGAAA4BdUKSPWoEEDNWrUSC6XS1dccYUaNmzobzExMbrzzjs1evTooEw0MjJSXbp0UXZ2tv+Yz+dTdna2kpOTK70mOTk54HxJWrFixRnPBwAABjidCautGbGPP/5YlmXp+uuv19/+9jc1atTI/15kZKSaN2+uJk2aGJ/kKenp6Ro2bJi6du2qbt26acaMGSorK9Pw4cMlSUOHDlXTpk2VmZkp6eQ3AFx33XV69tln1a9fPy1YsEAbN27Uyy+/HLQ5AgBwoasJT7Z3eny7qhSIXXfddZKknTt3qrCwUC+99JK2b9+ud955R02bNtWbb76pli1b6pprrgnKZAcNGqS9e/dqypQpKi4uVqdOnbRs2TL/hvzCwkKFhf0ryde9e3fNnz9fkyZN0kMPPaTLL79cixcvVvv27YMyPwAAgKqo1i1EGzdu1O23364hQ4Zo8+bN/sc9HDx4UNOmTQvq0+vT0tKUlpZW6XuVfc/lwIEDNXDgwKDNBwAA/JuaUBp0enybqnXX5OOPP66srCzNnj1bderU8R/v0aOHNm3aZGxyAAAgBDm9N6wmBII2VSsQKygo0G9/+9vTjkdHR+vAgQPnOicAAIALQrUCsfj4eG3btu2042vWrAnqd00CAICaz/GvN6oBNwvYVa1AbOTIkRozZozWrVsnl8ul77//XvPmzdP999+vUaNGmZ4jAAAIJaeerO90CwHV2qw/YcIE+Xw+9e7dW0eOHNFvf/tbud1u3X///br33ntNzxEAAISSmrBHy+nxbapWIOZyufTwww/rgQce0LZt23T48GG1a9dO9esH57sIAQAAaqNqlSZPiYyMVLt27dStWzeCMAAAIMn5vWHnskds1qxZatGihaKiopSUlKT169ef9fwZM2aodevWqlu3rhITEzVu3DgdO3bM9njnFIgBAACcxunHVlSzNLpw4UKlp6crIyNDmzZtUseOHZWamqo9e/ZUev78+fM1YcIEZWRk6KuvvtKrr76qhQsX6qGHHrI9JoEYAACApOnTp2vkyJEaPny42rVrp6ysLNWrV09z5syp9PxPP/1UPXr00G233aYWLVroxhtv1ODBg38xi/ZzBGIAAMCsmlCW/CkjVlpaGtBOfRvQv6uoqFBeXp5SUlL8x8LCwpSSkqLc3NxKr+nevbvy8vL8gdeOHTu0dOlS/cd//IftparWZn0AAIAzqkF3TSYmJgYczsjI0COPPHLa6fv27ZPX6/V/f/UpcXFx2rp1a6VD3Hbbbdq3b5+uueYaWZalEydO6J577qlSaZJADAAA1FpFRUXyeDz+126321jfOTk5mjZtml544QUlJSVp27ZtGjNmjB577DFNnjzZVh8EYgAAwKwalBHzeDwBgdiZxMTEKDw8XCUlJQHHS0pKFB8fX+k1kydP1u2336677rpLktShQweVlZXp7rvv1sMPP6ywsF/eAcYeMQAAYJTT+8Oq8/iKyMhIdenSRdnZ2f5jPp9P2dnZSk5OrvSaI0eOnBZshYeHS5Isy94EyIgBAABISk9P17Bhw9S1a1d169ZNM2bMUFlZmYYPHy5JGjp0qJo2barMzExJUv/+/TV9+nR17tzZX5qcPHmy+vfv7w/IfgmBGAAAgKRBgwZp7969mjJlioqLi9WpUyctW7bMv4G/sLAwIAM2adIkuVwuTZo0Sbt27VLjxo3Vv39/PfHEE7bHJBADAABm1aA9YlWVlpamtLS0St/LyckJeB0REaGMjAxlZGRUbzARiAEAAMPO5SuGTM4hFLBZHwAAwCFkxAAAgHkhkpFyGoEYAAAwK4T3iJ1vlCYBAAAcQkYMAAAYxWZ9+wjEAACAWZQmbaM0CQAA4BAyYgAAwChKk/YRiAEAALMoTdpGaRIAAMAhZMQAAIBZZMRsIxADAABGsUfMPgIxAABgFhkx29gjBgAA4BAyYgAAwCwyYrYRiAEAAKPYI2YfpUkAAACHkBEDAABmUZq0jUAMAAAYRWnSPkqTAAAADiEjBgAAzKI0aRuBGAAAMItAzDZKkwAAAA4hIwYAAIxy/dScnkMoIBADAABmUZq0jUAMAAAYxeMr7GOPGAAAgEPIiAEAALMoTdpGIAYAAMwLkUDIaZQmAQAAHEJGDAAAGMVmffsIxAAAgFnsEbON0iQAAIBDyIgBAACjKE3aRyAGAADMojRpG6VJAAAAh5ARAwAARlGatI9ADAAAmEVp0jYCMQAAYBaBmG3sEQMAAPjJrFmz1KJFC0VFRSkpKUnr168/6/kHDhzQ6NGjlZCQILfbrSuuuEJLly61PR4ZMQAAYFSo7hFbuHCh0tPTlZWVpaSkJM2YMUOpqakqKChQbGzsaedXVFTohhtuUGxsrN555x01bdpU3333nRo0aGB7TAIxAABgVoiWJqdPn66RI0dq+PDhkqSsrCx9+OGHmjNnjiZMmHDa+XPmzNH+/fv16aefqk6dOpKkFi1aVGlMSpMAAKDWKi0tDWjl5eWVnldRUaG8vDylpKT4j4WFhSklJUW5ubmVXvP3v/9dycnJGj16tOLi4tS+fXtNmzZNXq/X9vwIxAAAgFEuy6oRTZISExMVHR3tb5mZmZXOed++ffJ6vYqLiws4HhcXp+Li4kqv2bFjh9555x15vV4tXbpUkydP1rPPPqvHH3/c9lpRmgQAAGbVoNJkUVGRPB6P/7Db7TY2hM/nU2xsrF5++WWFh4erS5cu2rVrl5555hllZGTY6oNADAAA1FoejycgEDuTmJgYhYeHq6SkJOB4SUmJ4uPjK70mISFBderUUXh4uP9Y27ZtVVxcrIqKCkVGRv7iuJQmAQCAUafumnS6VUVkZKS6dOmi7Oxs/zGfz6fs7GwlJydXek2PHj20bds2+Xw+/7Gvv/5aCQkJtoIwiUAMAACYZtWQVkXp6emaPXu2Xn/9dX311VcaNWqUysrK/HdRDh06VBMnTvSfP2rUKO3fv19jxozR119/rQ8//FDTpk3T6NGjbY9JaRIAAEDSoEGDtHfvXk2ZMkXFxcXq1KmTli1b5t/AX1hYqLCwf+WwEhMTtXz5co0bN05XXnmlmjZtqjFjxujBBx+0PSaBGAAAMCpUH+gqSWlpaUpLS6v0vZycnNOOJScna+3atdUbTARiAADAtBp012RNRyAGAACMCuWM2PnGZn0AAACHkBEDAABmUZq0jUAMAAAYFyqlQadRmgQAAHAIGTEAAGCWZZ1sTs8hBBCIAQAAo7hr0j5KkwAAAA4hIwYAAMzirknbCMQAAIBRLt/J5vQcQgGlSQAAAIeQEQMAAGZRmrQtZDJiTzzxhLp376569eqpQYMGtq6544475HK5AlqfPn2CO1EAAC5wp+6adLqFgpDJiFVUVGjgwIFKTk7Wq6++avu6Pn366LXXXvO/drvdwZgeAAA4heeI2RYygdjUqVMlSXPnzq3SdW63W/Hx8UGYEQAAwLkJmdJkdeXk5Cg2NlatW7fWqFGj9MMPP5z1/PLycpWWlgY0AABgn9MlSUqTNUSfPn10yy23qGXLltq+fbseeugh9e3bV7m5uQoPD6/0mszMTH/27ecGbrtBEb8yW9as8FY+h3O1/esmQelXkcG5F7g0qm5Q+o36Pjgf74YF3qD0W9o8OJ+H8uTDQek3Zlbwfn0cTgzOFoLY/y0MSr+HOzcNSr9RQelVqrPrx6D0u/eehkHpt+Fnwfk7t6ffsaD06wpKr9Ku/6ow3qfvSIX0ofFu2axfBY5mxCZMmHDaZvp/b1u3bq12/7///e910003qUOHDhowYICWLFmiDRs2KCcn54zXTJw4UQcPHvS3oqKiao8PAABwNo5mxO677z7dcccdZz2nVatWxsZr1aqVYmJitG3bNvXu3bvSc9xuNxv6AQA4BzWhNOj0+HY5Gog1btxYjRs3Pm/j/d///Z9++OEHJSQknLcxAQC44HDXpG0hs1m/sLBQ+fn5KiwslNfrVX5+vvLz83X48L/2v7Rp00bvvfeeJOnw4cN64IEHtHbtWn377bfKzs7WzTffrMsuu0ypqalO/RgAAAB+IbNZf8qUKXr99df9rzt37ixJ+vjjj9WzZ09JUkFBgQ4ePChJCg8P1z//+U+9/vrrOnDggJo0aaIbb7xRjz32GKVHAACCiNKkfSETiM2dO/cXnyFm/SwNWbduXS1fvjzIswIAAKfhrknbQqY0CQAAUNuETEYMAACEBkqT9hGIAQAAs3zWyeb0HEIAgRgAADCLPWK2sUcMAADAIWTEAACAUS45v0crWN/5aRqBGAAAMIsn69tGaRIAAMAhZMQAAIBRPL7CPgIxAABgFndN2kZpEgAAwCFkxAAAgFEuy5LL4c3yTo9vF4EYAAAwy/dTc3oOIYDSJAAAgEPIiAEAAKMoTdpHIAYAAMzirknbKE0CAACzTj1Z3+lWDbNmzVKLFi0UFRWlpKQkrV+/3tZ1CxYskMvl0oABA6o0HoEYAACApIULFyo9PV0ZGRnatGmTOnbsqNTUVO3Zs+es13377be6//77de2111Z5TAIxAABg1Kkn6zvdqmr69OkaOXKkhg8frnbt2ikrK0v16tXTnDlzzniN1+vVkCFDNHXqVLVq1arKYxKIAQAAs5wuSf6sNFlaWhrQysvLK51yRUWF8vLylJKS4j8WFhamlJQU5ebmnvFHffTRRxUbG6sRI0ZUa6kIxAAAQK2VmJio6Ohof8vMzKz0vH379snr9SouLi7geFxcnIqLiyu9Zs2aNXr11Vc1e/bsas+PuyYBAIBRLt/J5vQcJKmoqEgej8d/3O12G+n/0KFDuv322zV79mzFxMRUux8CMQAAYNY53LVodA6SPB5PQCB2JjExMQoPD1dJSUnA8ZKSEsXHx592/vbt2/Xtt9+qf//+/mM+38noLyIiQgUFBbr00kt/cVxKkwAA4IIXGRmpLl26KDs723/M5/MpOztbycnJp53fpk0bbdmyRfn5+f520003qVevXsrPz1diYqKtccmIAQAAs0L0ga7p6ekaNmyYunbtqm7dumnGjBkqKyvT8OHDJUlDhw5V06ZNlZmZqaioKLVv3z7g+gYNGkjSacfPhkAMAAAYFapfcTRo0CDt3btXU6ZMUXFxsTp16qRly5b5N/AXFhYqLMxsMZFADAAA4CdpaWlKS0ur9L2cnJyzXjt37twqj0cgBgAAzKpBm/VrOgIxAABgliXJ4cdXOL5HzSYCMQAAYFSo7hFzAo+vAAAAcAgZMQAAYJYl5/dohUZCjEAMAAAYxmZ92yhNAgAAOISMGAAAMMsnyVUD5hACCMQAAIBR3DVpH6VJAAAAh5ARAwAAZrFZ3zYCMQAAYBaBmG2UJgEAABxCRgwAAJhFRsw2AjEAAGAWj6+wjUAMAAAYxeMr7GOPGAAAgEPIiAEAALPYI2YbgRgAADDLZ0kuhwMhX2gEYpQmAQAAHEJGDAAAmEVp0jYCMQAAYFgNCMTk9Pj2UJoEAABwCBkxAABgFqVJ2wjEAACAWT5LjpcGuWsSAAAAZ0NGDAAAmGX5Tjan5xACCMQAAIBZ7BGzjUAMAACYxR4x29gjBgAA4BAyYgAAwCxKk7YRiAEAALMsOR8IhUYcRmkSAADAKWTEAACAWZQmbSMQAwAAZvl8khx+jpcvNJ4jRmkSAADAIWTEAACAWZQmbSMQAwAAZhGI2UZpEgAAwCFkxAAAgFl8xZFtZMQAAIBRluWrEa06Zs2apRYtWigqKkpJSUlav379Gc+dPXu2rr32WjVs2FANGzZUSkrKWc+vDIEYAAAwy7JOZqScbNXYI7Zw4UKlp6crIyNDmzZtUseOHZWamqo9e/ZUen5OTo4GDx6sjz/+WLm5uUpMTNSNN96oXbt22R6TQAwAAEDS9OnTNXLkSA0fPlzt2rVTVlaW6tWrpzlz5lR6/rx58/SHP/xBnTp1Ups2bfTKK6/I5/MpOzvb9pjsEQMAAGZZNWCP2E8ZsdLS0oDDbrdbbrf7tNMrKiqUl5eniRMn+o+FhYUpJSVFubm5toY8cuSIjh8/rkaNGtmeJhkxAABgls9XM5qkxMRERUdH+1tmZmalU963b5+8Xq/i4uICjsfFxam4uNjWj/3ggw+qSZMmSklJsb1UZMQAAECtVVRUJI/H439dWTbMhCeffFILFixQTk6OoqKibF9HIAYAAMyqQaVJj8cTEIidSUxMjMLDw1VSUhJwvKSkRPHx8We99k9/+pOefPJJ/eMf/9CVV15ZpWlSmgQAAEZZPl+NaFURGRmpLl26BGy0P7XxPjk5+YzXPf3003rssce0bNkyde3atcprRUYMAABAUnp6uoYNG6auXbuqW7dumjFjhsrKyjR8+HBJ0tChQ9W0aVP/PrOnnnpKU6ZM0fz589WiRQv/XrL69eurfv36tsYkEAMAAGbVoNJkVQwaNEh79+7VlClTVFxcrE6dOmnZsmX+DfyFhYUKC/tXMfHFF19URUWF/vM//zOgn4yMDD3yyCO2xiQQAwAAZvksyRV6gZgkpaWlKS0trdL3cnJyAl5/++231Rrj59gjBgAA4BAyYgAAwCzLklS973o0O4eaj0AMAAAYZfksWQ6XJi0CMQAAcEGyfHI+I+bw+DaFzB6xm266Sc2aNVNUVJQSEhJ0++236/vvvz/rNceOHdPo0aN18cUXq379+rr11ltPe1AbAACAU0ImEOvVq5fefvttFRQU6G9/+5u2b99+2u2i/27cuHH64IMPtGjRIq1atUrff/+9brnllvM0YwAALkyWz6oRLRSETGly3Lhx/j83b95cEyZM0IABA3T8+HHVqVPntPMPHjyoV199VfPnz9f1118vSXrttdfUtm1brV27VldfffV5mzsAABcUSpO2hUwg9nP79+/XvHnz1L1790qDMEnKy8vT8ePHA74BvU2bNmrWrJlyc3PPGIiVl5ervLzc//rgwYOSpBNHKgz+BCed8IYb71OSfEePBaVfeYP0ofadCEq33mOuoPR74rg3KP16y4PzefAeCc7n4cSJ4P368FYE51+yJ3zlv3xSdfo9Hqw1Ds58FaR18AXrs+YN1nyD8zsilPiOnlxb0xvbT+i4489zPaHjzk7ALiuEjB8/3qpXr54lybr66qutffv2nfHcefPmWZGRkacd/81vfmONHz/+jNdlZGScehwwjUaj0WgXRNu+fbuR/04fPXrUio+Pd/znOdXi4+Oto0ePGvnZgsVlWc7d3zlhwgQ99dRTZz3nq6++Ups2bSRJ+/bt0/79+/Xdd99p6tSpio6O1pIlS+Rynf6vmvnz52v48OEB2S1J6tatm3r16nXGcf89I3bgwAE1b95chYWFio6OruqPiF9QWlqqxMREFRUVyePxOD2dWof1DT7WOLhY3+A6ePCgmjVrph9//FENGjQw0uexY8dUUWG+ilQdkZGRioqKcnoaZ+VoafK+++7THXfccdZzWrVq5f9zTEyMYmJidMUVV6ht27ZKTEzU2rVrK/1W9Pj4eFVUVOjAgQMBH66SkhLFx8efcTy32y23233a8ejoaH4JBJHH42F9g4j1DT7WOLhY3+D6+fcnnquoqKgaH/zUJI4GYo0bN1bjxo2rda3Pd3K/0r9nvE7p0qWL6tSpo+zsbN16662SpIKCAhUWFlYauAEAAJxvIbFZf926ddqwYYOuueYaNWzYUNu3b9fkyZN16aWX+oOqXbt2qXfv3nrjjTfUrVs3RUdHa8SIEUpPT1ejRo3k8Xh07733Kjk5mTsmAQBAjRASgVi9evX07rvvKiMjQ2VlZUpISFCfPn00adIkfxnx+PHjKigo0JEjR/zXPffccwoLC9Ott96q8vJypaam6oUXXqjS2G63WxkZGZWWK3HuWN/gYn2DjzUOLtY3uFhf5zm6WR8AAOBCFjJP1gcAAKhtCMQAAAAcQiAGAADgEAIxAAAAhxCIncWsWbPUokULRUVFKSkpSevXr3d6SiEpMzNTv/nNb3TRRRcpNjZWAwYMUEFBQcA5x44d0+jRo3XxxRerfv36uvXWW1VSUuLQjEPbk08+KZfLpbFjx/qPsb7nbteuXfrv//5vXXzxxapbt646dOigjRs3+t+3LEtTpkxRQkKC6tatq5SUFH3zzTcOzjh0eL1eTZ48WS1btlTdunV16aWX6rHHHgv4/kPWt2pWr16t/v37q0mTJnK5XFq8eHHA+3bWc//+/RoyZIg8Ho8aNGigESNG6PDhw+fxp7gwEIidwcKFC5Wenq6MjAxt2rRJHTt2VGpqqvbs2eP01ELOqlWrNHr0aK1du1YrVqzQ8ePHdeONN6qsrMx/zrhx4/TBBx9o0aJFWrVqlb7//nvdcsstDs46NG3YsEEvvfSSrrzyyoDjrO+5+fHHH9WjRw/VqVNHH330kb788ks9++yzatiwof+cp59+Wn/+85+VlZWldevW6Ve/+pVSU1N17Fhwvgi7Nnnqqaf04osv6i9/+Yu++uorPfXUU3r66af1/PPP+89hfaumrKxMHTt21KxZsyp93856DhkyRF988YVWrFihJUuWaPXq1br77rvP149w4XDuay5rtm7dulmjR4/2v/Z6vVaTJk2szMxMB2dVO+zZs8eSZK1atcqyLMs6cOCAVadOHWvRokX+c7766itLkpWbm+vUNEPOoUOHrMsvv9xasWKFdd1111ljxoyxLIv1NeHBBx+0rrnmmjO+7/P5rPj4eOuZZ57xHztw4IDldrutt95663xMMaT169fPuvPOOwOO3XLLLdaQIUMsy2J9z5Uk67333vO/trOeX375pSXJ2rBhg/+cjz76yHK5XNauXbvO29wvBGTEKlFRUaG8vDylpKT4j4WFhSklJUW5ubkOzqx2OHjwoCSpUaNGkqS8vDwdP348YL3btGmjZs2asd5VMHr0aPXr1y9gHSXW14S///3v6tq1qwYOHKjY2Fh17txZs2fP9r+/c+dOFRcXB6xxdHS0kpKSWGMbunfvruzsbH399deSpM8++0xr1qxR3759JbG+ptlZz9zcXDVo0EBdu3b1n5OSkqKwsDCtW7fuvM+5NguJJ+ufb/v27ZPX61VcXFzA8bi4OG3dutWhWdUOPp9PY8eOVY8ePdS+fXtJUnFxsSIjIwO+nF06ud7FxcUOzDL0LFiwQJs2bdKGDRtOe4/1PXc7duzQiy++qPT0dD300EPasGGD/vjHPyoyMlLDhg3zr2NlvzNY4182YcIElZaWqk2bNgoPD5fX69UTTzyhIUOGSBLra5id9SwuLlZsbGzA+xEREWrUqBFrbhiBGM6r0aNH6/PPP9eaNWucnkqtUVRUpDFjxmjFihWKiopyejq1ks/nU9euXTVt2jRJUufOnfX5558rKytLw4YNc3h2oe/tt9/WvHnzNH/+fP36179Wfn6+xo4dqyZNmrC+qPUoTVYiJiZG4eHhp91VVlJSovj4eIdmFfrS0tK0ZMkSffzxx7rkkkv8x+Pj41VRUaEDBw4EnM9625OXl6c9e/boqquuUkREhCIiIrRq1Sr9+c9/VkREhOLi4ljfc5SQkKB27doFHGvbtq0KCwslyb+O/M6ongceeEATJkzQ73//e3Xo0EG33367xo0bp8zMTEmsr2l21jM+Pv60m9NOnDih/fv3s+aGEYhVIjIyUl26dFF2drb/mM/nU3Z2tpKTkx2cWWiyLEtpaWl67733tHLlSrVs2TLg/S5duqhOnToB611QUKDCwkLW24bevXtry5Ytys/P97euXbtqyJAh/j+zvuemR48epz1y5euvv1bz5s0lSS1btlR8fHzAGpeWlmrdunWssQ1HjhxRWFjgf47Cw8Pl8/kksb6m2VnP5ORkHThwQHl5ef5zVq5cKZ/Pp6SkpPM+51rN6bsFaqoFCxZYbrfbmjt3rvXll19ad999t9WgQQOruLjY6amFnFGjRlnR0dFWTk6OtXv3bn87cuSI/5x77rnHatasmbVy5Upr48aNVnJyspWcnOzgrEPbz++atCzW91ytX7/eioiIsJ544gnrm2++sebNm2fVq1fP+utf/+o/58knn7QaNGhgvf/++9Y///lP6+abb7ZatmxpHT161MGZh4Zhw4ZZTZs2tZYsWWLt3LnTevfdd62YmBhr/Pjx/nNY36o5dOiQtXnzZmvz5s2WJGv69OnW5s2bre+++86yLHvr2adPH6tz587WunXrrDVr1liXX365NXjwYKd+pFqLQOwsnn/+eatZs2ZWZGSk1a1bN2vt2rVOTykkSaq0vfbaa/5zjh49av3hD3+wGjZsaNWrV8/63e9+Z+3evdu5SYe4fw/EWN9z98EHH1jt27e33G631aZNG+vll18OeN/n81mTJ0+24uLiLLfbbfXu3dsqKChwaLahpbS01BozZozVrFkzKyoqymrVqpX18MMPW+Xl5f5zWN+q+fjjjyv9vTts2DDLsuyt5w8//GANHjzYql+/vuXxeKzhw4dbhw4dcuCnqd1clvWzRxcDAADgvGGPGAAAgEMIxAAAABxCIAYAAOAQAjEAAACHEIgBAAA4hEAMAADAIQRiAAAADiEQAwAAcAiBGIAaq0WLFpoxY4bT0wCAoCEQAwAAcAhfcQTAMT179lT79u0lSW+++abq1KmjUaNG6dFHH1WvXr20atWqgPP5dQWgtiEjBsBRr7/+uiIiIrR+/XrNnDlT06dP1yuvvKJ3331Xl1xyiR599FHt3r1bu3fvdnqqAGBchNMTAHBhS0xM1HPPPSeXy6XWrVtry5Yteu655zRy5EiFh4froosuUnx8vNPTBICgICMGwFFXX321XC6X/3VycrK++eYbeb1eB2cFAOcHgRgAAIBDCMQAOGrdunUBr9euXavLL79c4eHhioyMJDMGoFYjEAPgqMLCQqWnp6ugoEBvvfWWnn/+eY0ZM0bSyeeIrV69Wrt27dK+ffscnikAmMfjKwA4pmfPnvr1r38tn8+n+fPnKzw8XKNGjdLjjz8ul8ultWvX6n/+539UUFCg8vJyHl8BoNYhEAPgmJ49e6pTp048PR/ABYvSJAAAgEMIxAAAABxCaRIAAMAhZMQAAAAcQiAGAADgEAIxAAAAhxCIAQAAOIRADAAAwCEEYgAAAA4hEAMAAHAIgRgAAIBDCMQAAAAc8v8iG0DLDOqMOQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from coffea.lookup_tools.dense_lookup import dense_lookup\n",
    "\n",
    "num = dists[\"gen1\", :, :].values()\n",
    "den = dists[\"gen2\", :, :].values()\n",
    "sf = np.where(\n",
    "    (num > 0) & (den > 0),\n",
    "    num / np.maximum(den, 1) * den.sum() / num.sum(),\n",
    "    1.0,\n",
    ")\n",
    "\n",
    "corr = dense_lookup(sf, [ax.edges for ax in dists.axes[1:]])\n",
    "print(corr)\n",
    "\n",
    "# a quick way to plot the scale factor is to steal the axis definitions from the input histograms:\n",
    "sfhist = hist.Hist(*dists.axes[1:], data=sf)\n",
    "sfhist.plot2d()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we generate some new mock data as if it was drawn from `gen2` and reweight it with our `corr` to match `gen1`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x319d25290>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGwCAYAAAC3qV8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOJElEQVR4nO3de1xU1d4/8M+GmYEZ7qDcCkUt75B4IyoNkyOiebLsppZYptXBSlEzyxS8gUqa9XT0yZPo6eixen5qHjVvJGqFqBiJmpSm4UkupcltRmaA/fuD2DIJwwzMMGz4vF+vebX3Xt+995qtwde11l5LEEVRBBEREZGMONi7AkRERESWYgJDREREssMEhoiIiGSHCQwRERHJDhMYIiIikh0mMERERCQ7TGCIiIhIdhT2roCtVFdX4+rVq3Bzc4MgCPauDhEREZlBFEWUlpYiMDAQDg4Nt7O02QTm6tWrCAoKsnc1iIiIqAmuXLmCO++8s8HyNpvAuLm5Aah5AO7u7nauDREREZmjpKQEQUFB0u/xhrTZBKa228jd3Z0JDBERkcw0NvyDg3iJiIhIdpjAEBERkewwgSEiIiLZabNjYIiIqH2qqqqCwWCwdzWoAUqlEo6Ojs2+DhMYIiJqE0RRREFBAW7cuGHvqlAjPD094e/v36x52pjAEBFRm1CbvPj6+kKj0XAS01ZIFEVotVoUFRUBAAICApp8LSYwREQke1VVVVLy4uPjY+/qkAlqtRoAUFRUBF9f3yZ3J3EQLxERyV7tmBeNRmPnmpA5av+cmjNWiQkMERG1Gew2kgdr/DkxgSEiIiLZYQJDREREssMEhoiIyMoiIyMxY8YMe1ejTWMCQ0REZEfp6ekQBKHF569JSEhAv379WvSe1sQEhoiIiGSHCQwREVEzlJeXY9KkSXB1dUVAQADeeecdo/KPP/4YAwcOhJubG/z9/TFhwgRpIrfLly9j2LBhAAAvLy8IgoDJkycDAPbu3YsHHngAnp6e8PHxwcMPP4yLFy9K19Xr9Zg+fToCAgLg7OyMzp07IykpSSq/ceMGXnjhBXTs2BHu7u546KGH8N133wEANm7ciMTERHz33XcQBAGCIGDjxo02fErWxwRGZrQGLUI2hSBkUwi0Bq29q0NE1O7NmTMHhw8fxueff479+/cjPT0dp06dksoNBgMWL16M7777Djt27MDly5elJCUoKAj/7//9PwBAbm4u8vPzsWbNGgA1iVF8fDxOnjyJtLQ0ODg44NFHH0V1dTUA4L333sPOnTvx6aefIjc3F5s3b0ZwcLB03yeeeAJFRUX44osvkJWVhf79+2P48OG4fv06nnrqKcyaNQt9+vRBfn4+8vPz8dRTT7XMA7MSzsRLRETURGVlZfjoo4/wr3/9C8OHDwcAbNq0CXfeeacU8/zzz0vbXbt2xXvvvYdBgwahrKwMrq6u8Pb2BgD4+vrC09NTih03bpzRvTZs2ICOHTvi3Llz6Nu3L/Ly8nD33XfjgQcegCAI6Ny5sxT71Vdf4fjx4ygqKoKTkxMAICUlBTt27MD//d//Ydq0aXB1dYVCoYC/v7/Vn0tLYAsMERFRE128eBF6vR7h4eHSMW9vb/To0UPaz8rKwpgxY9CpUye4ubnhwQcfBADk5eWZvPaPP/6I8ePHo2vXrnB3d5daV2rPmzx5MrKzs9GjRw+8+uqr2L9/v3Tud999h7KyMvj4+MDV1VX6XLp0yagbSs7YAkNERGQj5eXliI6ORnR0NDZv3oyOHTsiLy8P0dHR0Ov1Js8dM2YMOnfujPXr1yMwMBDV1dXo27evdF7//v1x6dIlfPHFFzh48CCefPJJREVF4f/+7/9QVlaGgIAApKen33bduq08csYEhoiIqIm6desGpVKJzMxMdOrUCQDw+++/44cffsCDDz6I8+fP49q1a0hOTkZQUBAA4OTJk0bXUKlUAGoWpKx17do15ObmYv369RgyZAiAmm6hP3N3d8dTTz2Fp556Co8//jhGjhyJ69evo3///igoKIBCoTAaF/Pn+9a9p9wwgZEZnaHKaFujtGNliIjaOVdXV0yZMgVz5syBj48PfH198dZbb8HBoWaERqdOnaBSqfD+++/jpZdewpkzZ7B48WKja3Tu3BmCIGDXrl0YNWoU1Go1vLy84OPjgw8//BABAQHIy8vDG2+8YXTeqlWrEBAQgLCwMDg4OOCzzz6Dv78/PD09ERUVhYiICIwdOxYrVqxA9+7dcfXqVezevRuPPvooBg4ciODgYFy6dAnZ2dm488474ebmJo2XkQOOgSEAwDVtqfR20zVtqb2rQ0QkGytXrsSQIUMwZswYREVF4YEHHsCAAQMAAB07dsTGjRvx2WefoXfv3khOTkZKSorR+XfccQcSExPxxhtvwM/PD9OnT4eDgwO2bt2KrKws9O3bFzNnzsTKlSuNznNzc8OKFSswcOBADBo0CJcvX8aePXvg4OAAQRCwZ88eDB06FM899xy6d++Op59+Gj///DP8/PwA1AwSHjlyJIYNG4aOHTvi3//+d8s8MCsRRFEU7V0JWygpKYGHhweKi4vh7u5u7+pYzTVtKSI/uw8AkP7EN/DRuLXq6xIRtYSbN2/i0qVL6NKlC5ydne1dHWqEqT8vc39/swupFRBFEbpKnVmxWoPOaFttcDTrPLVCzWXmiYiozWAC0wroKnUI3xLeeCAAsVoB4Y+Ov5ht0RAcKs06L3NCJjRKTVOrSERE1KpwDAwRERHJDltgWoG6w5C+eDQNaoW6wdjftWV4dPcIAMCOMXvgpXFtMFZXqUPM9uG33YOIiEjumMC0Ajcrq6XtB5K+BkRVw8GCHm49azaj3jlmduzNymq4mAglIiKSEyYw7cR1bZnJ8t/rlP/eSGxdXs4u0nwHRERELcXiBObIkSNYuXIlsrKykJ+fj+3bt2Ps2LFSeUNvuqxYsQJz5swBAAQHB+Pnn382Kk9KSjKapOf06dOIi4vDiRMn0LFjR7zyyit4/fXXLa2u7BydOww+JrqFrmlLMerz2thIk687X9OWSbGP/dHt1JC6g4PH/meU2YOD+co1ERHZg8UJTHl5Oe655x48//zzeOyxx24rz8/PN9r/4osvMGXKlNtW1Vy0aBGmTp0q7bu53folWFJSghEjRiAqKgrr1q1DTk4Onn/+eXh6emLatGmWVllWNCpHaFQN/7HoKhV1YhWNxJr3ijUREd2i1Vei94J9AIBzi6JN/pwl+7H4TyUmJgYxMTENlv95We7PP/8cw4YNQ9euXY2Ou7m5NbiE9+bNm6HX67FhwwaoVCr06dMH2dnZWLVqVYMJTEVFBSoqKqT9kpISc79Sm+Xl7IL0J74xK9aSwcHXtWWNtugQERHZkk0HLxQWFmL37t2YMmXKbWXJycnw8fFBWFgYVq5cicrKW10WGRkZGDp0qLTAFQBER0cjNzcXv//+e733SkpKgoeHh/SpXTSrPXNwcICPxs2sT92ExUvjajLW20RyQ0RErdPZs2cxbtw4BAcHQxAEvPvuu/auUrPYNIHZtGkT3NzcbutqevXVV7F161YcOnQIL774IpYtW2Y0vqWgoEBaq6FW7X5BQUG995o3bx6Ki4ulz5UrV6z8bYiIiORLq9Wia9euSE5ObrAHRE5smsBs2LABEydOvG2dg/j4eERGRiI0NBQvvfQS3nnnHbz//vtGXUCWcnJygru7u9GnLVIr1Cj9Phml3yebnC+GiIjko7S0FBMnToSLiwsCAgKwevVqREZGYsaMGQBqhknMnj0bd9xxB1xcXBAeHo709HTp/I0bN8LT0xP79u1Dr1694OrqipEjRxqNSx00aBBWrlyJp59+WlarTjfEZiOTjh49itzcXHzyySeNxoaHh6OyshKXL19Gjx494O/vj8LCQqOY2v22kDUSEZHtiaIInaHKZIxWf/sblzr9rXOulVVAq7o9prGBvWqlo0Xrz8XHx+Prr7/Gzp074efnhwULFuDUqVPo168fAGD69Ok4d+4ctm7disDAQGzfvh0jR45ETk4O7r777prvotUiJSUFH3/8MRwcHPDMM89g9uzZ2Lx5s9n1kBObJTAfffQRBgwYgHvuuafR2OzsbDg4OMDX1xcAEBERgbfeegsGgwFKpRIAcODAAfTo0QNeXl62qrIsaFQKXE4ebe9qEBG1ejpDlfQ2UVMNWZHepPMseXuptLQUmzZtwpYtWzB8eM3s6ampqQgMDAQA5OXlITU1FXl5edKx2bNnY+/evUhNTcWyZcsAAAaDAevWrUO3bt0A1CQ9ixYtalL95cDiBKasrAwXLlyQ9i9duoTs7Gx4e3ujU6dOAGreAPrss8/wzjvv3HZ+RkYGMjMzMWzYMLi5uSEjIwMzZ87EM888IyUnEyZMQGJiIqZMmYK5c+fizJkzWLNmDVavXt3U70lERNQq/fTTTzAYDBg8eLB0zMPDAz169AAA5OTkoKqqCt27dzc6r6KiAj4+PtK+RqORkhcACAgIQFFRkY1rbz8WJzAnT57EsGHDpP34+HgAQGxsLDZu3AgA2Lp1K0RRxPjx428738nJCVu3bkVCQgIqKirQpUsXzJw5U7oOUPMHt3//fsTFxWHAgAHo0KEDFixY0ObngLEntdKx3m0iIrlSKx1xblG0yZiGupBqW16Ovh4Jter2n4nmdCFZS1lZGRwdHZGVlQVHR+Prurreeiu0tseiliAIbXodPIsTmMjIyEYfyLRp0xpMNvr3749jx441ep/Q0FAcPXrU0upRE2mUGuTE5ti7GkREViMIQqOJRn3ldZMaH1cnm09k17VrVyiVSpw4cULqySguLsYPP/yAoUOHIiwsDFVVVSgqKsKQIUNsWhc54fSCREREduTm5obY2FjMmTMH3t7e8PX1xcKFC+Hg4ABBENC9e3dMnDgRkyZNwjvvvIOwsDD8+uuvSEtLQ2hoKEaPNm9cpF6vx7lz56TtX375BdnZ2XB1dcVdd91ly69oE1yFj4iIyM5WrVqFiIgIPPzww4iKisL999+PXr16SdOQpKamYtKkSZg1axZ69OiBsWPHGrXYmOPq1asICwtDWFgY8vPzkZKSgrCwMLzwwgu2+lo2xRYYIiIiO3NzczN63bm8vByJiYnScAylUonExEQkJibWe/7kyZMxefJko2Njx441GvIRHBzcpsbEMIEhIiKys2+//Rbnz5/H4MGDUVxcLL3+/Mgjj9i5Zq0XExgiIqI67DXfVkpKCnJzc6FSqTBgwAAcPXoUHTp0aPF6yAUTGCIiIjsLCwtDVlaWvashKxzES0RERLLDBIaIiIhkhwkMERERyQ4TGCIiIpIdJjBEREQkO0xgiIiI6tKXAwkeNR99ub1rQw1gAkNERESywwSGiIioHVi/fj2GDBkCLy8veHl5ISoqCsePH7d3tZqMCQwREVE7kJ6ejvHjx+PQoUPIyMhAUFAQRowYgV9++cXeVWsSJjBERER2VlpaiokTJ8LFxQUBAQFYvXo1IiMjMWPGDABARUUFZs+ejTvuuAMuLi4IDw9Henq6dP7GjRvh6emJffv2oVevXnB1dcXIkSORn58vxWzevBl/+9vf0K9fP/Ts2RP/+Mc/UF1djbS0tBb+ttbBpQSIiKhtEkXAoDUdo6+nvO6xst8AVT0xKo3p6yo1gCA0Xsc/xMfH4+uvv8bOnTvh5+eHBQsW4NSpU+jXrx8AYPr06Th37hy2bt2KwMBAbN++HSNHjkROTg7uvvtuAIBWq0VKSgo+/vhjODg44JlnnsHs2bONVrmuS6vVwmAwwNvb2+x6tiZMYIiIqG0yaIFlgc27xnuhTTvvzauAysWs0NLSUmzatAlbtmzB8OHDAQCpqakIDKype15eHlJTU5GXlycdmz17Nvbu3YvU1FQsW7YMAGAwGLBu3Tp069YNQE3SU7uqdX3mzp2LwMBAREVFNe072hkTGCIiIjv66aefYDAYMHjwYOmYh4cHevToAQDIyclBVVUVunfvbnReRUUFfHx8pH2NRiMlLwAQEBCAoqKieu+ZnJyMrVu3Ij09Hc7Oztb8Oi2GCQwREbVNSk1NS4gpDXUh1ba8vHq6/u4ic7qQrKSsrAyOjo7IysqCo6OjUZmrq+utWyqVRmWCIEAUxduul5KSguTkZBw8eBChoU1sYWoFmMAQEVHbJAiNd+PUV1538jrXDmZ3BTVV165doVQqceLECXTq1AkAUFxcjB9++AFDhw5FWFgYqqqqUFRUhCFDhjTrXitWrMDSpUuxb98+DBw40BrVtxsmMERERHbk5uaG2NhYzJkzB97e3vD19cXChQvh4OAAQRDQvXt3TJw4EZMmTcI777yDsLAw/Prrr0hLS0NoaChGjx5t1n2WL1+OBQsWYMuWLQgODkZBQQGAmlacui05csHXqImIiOxs1apViIiIwMMPP4yoqCjcf//96NWrlzQ+JTU1FZMmTcKsWbPQo0cPjB071qjFxhxr166FXq/H448/joCAAOmTkpJiq69lU2yBISIisjM3Nzej153Ly8uRmJiIadOmAagZ35KYmIjExMR6z588eTImT55sdGzs2LFGY2AuX75s9XrbExMYIiIiO/v2229x/vx5DB48GMXFxdLrz4888oida9Z6MYEhIiKqS+UCJBS3+G1TUlKQm5sLlUqFAQMG4OjRo+jQoUOL10MumMAQERHZWVhYGLKysuxdDVnhIF4iIiKSHSYwZFNagxYhm0IQsikE2sbWJCEiIjITExgiIiKSHSYwREREJDtMYIiIiEh2mMAQERHVwbF78sAEhoiIiGSHCQw1i65SB61B2+DnmrZMir2mLTMZW/dT3xLwRETUdOvXr8eQIUPg5eUFLy8vREVF4fjx41a/T2RkJGbMmGH16/4ZJ7KjGvpyYFlgzfabV81ePj5m+3CT5WK1AsIfaXLMtmgIDpVmXTdzQiY0So1ZsURE1Lj09HSMHz8e9913H5ydnbF8+XKMGDECZ8+exR133FHvOXq9HiqVqoVrah6LW2COHDmCMWPGIDAwEIIgYMeOHUblkydPhiAIRp+RI0caxVy/fh0TJ06Eu7s7PD09MWXKFJSVlRnFnD59GkOGDIGzszOCgoKwYsUKy78dERGRDJSWlmLixIlwcXFBQEAAVq9ebdSSUVFRgdmzZ+OOO+6Ai4sLwsPDkZ6eLp2/ceNGeHp6Yt++fejVqxdcXV0xcuRI5OfnSzGbN2/G3/72N/Tr1w89e/bEP/7xD1RXVyMtLU2KCQ4OxuLFizFp0iS4u7tj2rRpePzxxzF9+nQpZsaMGRAEAefPnwdQk+S4uLjg4MGDmDx5Mg4fPow1a9ZIOYCtFpG0uAWmvLwc99xzD55//nk89thj9caMHDkSqamp0r6Tk5NR+cSJE5Gfn48DBw7AYDDgueeew7Rp07BlyxYAQElJCUaMGIGoqCisW7cOOTk5eP755+Hp6SmtzEn24+zojNLzNQuNZb0dBbXSscHYa9oyjNpR00rzxWP74KNxbTBWV6lD5KeRVq0rEbVfoihCV6kzGVNfed1j129erzdGrVCbvK5aoYYgCGbWFIiPj8fXX3+NnTt3ws/PDwsWLMCpU6fQr18/AMD06dNx7tw5bN26FYGBgdi+fTtGjhyJnJwc3H333QAArVaLlJQUfPzxx3BwcMAzzzyD2bNnG61yXZdWq4XBYIC3t7fR8ZSUFCxYsAALFy4EAOzZswf/+7//K5UfPnwYHTp0QHp6Onr27IkTJ07AYDDgvvvuw6BBg/DDDz+gb9++0oKUHTt2NPs5WMLiBCYmJgYxMTEmY5ycnODv719v2ffff4+9e/fixIkTGDhwIADg/fffx6hRo5CSkoLAwEBs3rwZer0eGzZsgEqlQp8+fZCdnY1Vq1YxgWkFBEEAxJomRbVCDY2y4b9GOmWVtK1RqtktREQtRlepQ/iW8GZdI2ab6d93DbGkG7y0tBSbNm3Cli1bMHx4zT/4UlNTERhY062fl5eH1NRU5OXlScdmz56NvXv3IjU1FcuWLQMAGAwGrFu3Dt26dQNQk/TUJhH1mTt3LgIDAxEVFWV0/KGHHsKsWbOk/cjISLz22mv49ddfoVAocO7cObz99ttIT0/HSy+9hPT0dAwaNAgaTc33ValU0Gg0DeYB1mKTMTDp6enw9fWFl5cXHnroISxZsgQ+Pj4AgIyMDHh6ekrJCwBERUXBwcEBmZmZePTRR5GRkYGhQ4ca9btFR0dj+fLl+P333+Hl5XXbPSsqKlBRUSHtl5SU2OKr0Z9o9VWNlFcabasVDY+B0VXeuhYH8RJRe/HTTz/BYDBg8ODB0jEPDw/06NEDAJCTk4Oqqip0797d6LyKigrpdysAaDQaKXkBgICAABQVFdV7z+TkZGzduhXp6elwdnY2Kqv7+xkA+vbtC29vbxw+fBgqlQphYWF4+OGH8cEHHwCoaZGJjIy0/Is3k9UTmJEjR+Kxxx5Dly5dcPHiRbz55puIiYlBRkYGHB0dUVBQAF9fX+NKKBTw9vZGQUEBAKCgoABdunQxivHz85PK6ktgkpKSkJiYaO2vI2+iCJg7h4FeW/92vbGVUOMmAGDgkgMATDSTCnq49azZHLI8XWq5aSz2ZmU1XFrnuDEikgm1Qo3MCZkmYxrqQqptefnisS/q7S4ypwvJWsrKyuDo6IisrCw4Ohp32bu63uqWVyqVRmWCINT7j8GUlBQkJyfj4MGDCA0Nva3cxcX4JQ5BEDB06FCkp6fDyckJkZGRCA0NRUVFBc6cOYNvvvkGs2fPbs5XbBKrJzBPP/20tB0SEoLQ0FB069YN6enpUtOYLcybNw/x8fHSfklJCYKCgmx2P1kwaG+9WWSJlLtMFmsAfP9Hwt7r5gbo4GwynojIHgRBaLQbp77yupPXeTt727zru2vXrlAqlThx4gQ6deoEACguLsYPP/yAoUOHIiwsDFVVVSgqKsKQIUOada8VK1Zg6dKl2Ldv320tLaY8+OCDWL9+PZycnLB06VI4ODhg6NChWLlyJSoqKnD//fdLsSqVClVVplvnrcHmr1F37doVHTp0wIULFzB8+HD4+/vf1qRVWVmJ69evS/1l/v7+KCwsNIqp3W+oT83Jyem2wcJke1nzo0y+cn1NW4pRn9dsH50bCR+Nm4nYMimWiKi9cHNzQ2xsLObMmQNvb2/4+vpi4cKFcHBwgCAI6N69OyZOnIhJkybhnXfeQVhYGH799VekpaUhNDQUo0ePNus+y5cvx4IFC7BlyxYEBwdLvR6urq5GLTn1iYyMxMyZM6FSqfDAAw9Ix2bPno1BgwYZtdoEBwcjMzMTly9fhqurK7y9veHgYP1p52yewPz3v//FtWvXEBAQAACIiIjAjRs3kJWVhQEDBgAAvvzyS1RXVyM8PFyKeeutt2AwGKQmsQMHDqBHjx71dh+RGWZfAFQm/hWh195qebEgVqNSACoTg3grb5VpVIqa+AZjG36biYioLVu1ahVeeuklPPzww3B3d8frr7+OK1euSONTUlNTsWTJEsyaNQu//PILOnTogHvvvRcPP/yw2fdYu3Yt9Ho9Hn/8caPjCxcuREJCgslzQ0JC4Onpie7du0vJTmRkJKqqqm4b/zJ79mzExsaid+/e0Ol0uHTpEoKDg82up7ksTmDKyspw4cIFaf/SpUvIzs6Gt7c3vL29kZiYiHHjxsHf3x8XL17E66+/jrvuugvR0dEAgF69emHkyJGYOnUq1q1bB4PBgOnTp+Ppp5+WRldPmDABiYmJmDJlCubOnYszZ85gzZo1WL16tZW+djuk0pg9OZ1FsY1QK9Qo/T5Z2iYiotu5ubkZve5cXl6OxMRE6c1bpVKJxMTEBsd6Tp48GZMnTzY6NnbsWKMxMObMx9JQjIODA65fv250rF+/fvWOsenevTsyMjIavVdzWZzAnDx5EsOGDZP2a8edxMbGYu3atTh9+jQ2bdqEGzduIDAwECNGjMDixYuNunc2b96M6dOnY/jw4XBwcMC4cePw3nvvSeUeHh7Yv38/4uLiMGDAAHTo0AELFizgK9RERNQmffvttzh//jwGDx6M4uJi6fXnRx55xM41a70sTmAiIyNNvuK6b9++Rq/h7e0tTVrXkNDQUBw9etTS6hERETWLRqlBTmxOi983JSUFubm5UKlUGDBgAI4ePYoOHTq0eD3kgmshERER2VlYWBiysrLsXQ1Z4WrUREREJDtsgaEaKhcgodjql9WoFLicbN4rfkREzcVZvOXBGn9ObIEhIiLZq51yQ6s1c/ZxsqvaP6c/zx5sCbbAEBGR7Dk6OsLT01OaKFWj0Vi0GjS1DFEUodVqUVRUBE9Pz9uWRrAEExgiImoTamdqb2gBQ2o9PD09m71aNRMYIiJqEwRBQEBAAHx9fWEwGOxdHWqAUqlsVstLLSYwRETUpjg6OlrlFyS1bhzES0RERLLDBIaIiIhkhwkMERERyQ4TGCIiIpIdJjBEREQkO0xgiIiISHaYwJBt6cuBBI+aj77c3rUhIqI2ggkMERERyQ4TGCIiIpIdJjBEREQkO0xgiIiISHa4FhI1j15rfnmjsXUG+Ypi0+tERERtHhMYap6Uu6wWqxYEIDioZsegBeDe9HoREVGbxi4kIiIikh22wJDllBrgzavmxeq1t1peZl8AVJoGQ3XFhcCeR6xQQSIiauuYwJDlBAFQuVh+nkpj+jxlw8kNERFRXexCIiIiItlhAkNERESywwSGiIiIZIcJjNxwcUQiIiIO4iUbU7kACcX2rgUREbUxbIEhIiIi2WECQ0RERLLDBIaIiIhkh2NgWoO6CxfqywGFibzSosURGyknIiKSKSYwrYHhVqKhXtPT/JWYLVlIkYiIqA1hFxLJktagRcimEIRsCoHWwJYmIqL2hi0wrYzub6eg8fBrOMCCxRGNcJ0hIiJqQ5jAtDbKRhY8rKuxxRGJiIjaKIu7kI4cOYIxY8YgMDAQgiBgx44dUpnBYMDcuXMREhICFxcXBAYGYtKkSbh69arRNYKDgyEIgtEnOTnZKOb06dMYMmQInJ2dERQUhBUrVjTtGxIREVGbY3ECU15ejnvuuQcffPDBbWVarRanTp3C22+/jVOnTmHbtm3Izc3FX//619tiFy1ahPz8fOnzyiuvSGUlJSUYMWIEOnfujKysLKxcuRIJCQn48MMPLa0uERERtUEWdyHFxMQgJiam3jIPDw8cOHDA6Nj//M//YPDgwcjLy0OnTp2k425ubvD396/3Ops3b4Zer8eGDRugUqnQp08fZGdnY9WqVZg2bZqlVSYZ0lXeNDk4V1epq3e7MWqFGoIgNKtuRERkfzYfA1NcXAxBEODp6Wl0PDk5GYsXL0anTp0wYcIEzJw5EwpFTXUyMjIwdOhQqFQqKT46OhrLly/H77//Di8vr9vuU1FRgYqKCmm/pKTENl+IWkTM7jFmx0Z+Gml2bOaETGg4oJmISPZsmsDcvHkTc+fOxfjx4+Hu7i4df/XVV9G/f394e3vjm2++wbx585Cfn49Vq1YBAAoKCtClSxeja/n5+Ull9SUwSUlJSExMtOG3aSW4OCIREZHtEhiDwYAnn3wSoihi7dq1RmXx8fHSdmhoKFQqFV588UUkJSXBycmpSfebN2+e0XVLSkoQFBTUtMqTXTg7OiPz8hUAgO6181C7uDcYq6vUSS0v6U+mQ61QmxVLRERtg00SmNrk5eeff8aXX35p1PpSn/DwcFRWVuLy5cvo0aMH/P39UVhYaBRTu9/QuBknJ6cmJz/UOgiCAE3tLMQKtcmuHp2h6taOqGK3EBFRO2P1mXhrk5cff/wRBw8ehI+PT6PnZGdnw8HBAb6+vgCAiIgIHDlyBAaDQYo5cOAAevToUW/3EREREbUvFrfAlJWV4cKFC9L+pUuXkJ2dDW9vbwQEBODxxx/HqVOnsGvXLlRVVaGgoAAA4O3tDZVKhYyMDGRmZmLYsGFwc3NDRkYGZs6ciWeeeUZKTiZMmIDExERMmTIFc+fOxZkzZ7BmzRqsXr3aSl+biIiI5MziBObkyZMYNmyYtF877iQ2NhYJCQnYuXMnAKBfv35G5x06dAiRkZFwcnLC1q1bkZCQgIqKCnTp0gUzZ840Gr/i4eGB/fv3Iy4uDgMGDECHDh2wYMECvkJNREREAJqQwERGRkI0sVqyqTIA6N+/P44dO9bofUJDQ3H06FFLq0dERETtANdColZJq68C9JUNlovVSpR+nyxta03E6ipvDfhtLMEmIiJ5YAJDrdKQFYegg7NZsQOXpJkOEPRw61mzebOyGi4q0+FERNT6Wf0tJCIiIiJbYwsMtRpqpaO0nTU/qmbW4QZoy4rR4b2uAIDfXv0JGlePBmOvacsw6nPr1ZOIiOyPCQy1GnUXWdSoFIDKxF/POmUalaImvgG6SscGy4iISJ7YhURERESywwSGiIiIZIddSNQ66bWNlOuMt/Um/irry29t8zVqIqI2gQkMtU4pd5kuF1VA7ZCZNSGAoG8wVC0IQPAfK5MbtABMLy5KREStH7uQSJY0dRIWjYnkhYiI2ia2wFDrodQAb141L1avvdVKM/sCoNI0GKorLgT2PGKFChIRUWvBBIZaD0EwOfdLg1Qa0+cpG05uiIhIntiFRERERLLDBIaIiIhkhwkMERERyQ4TGCIiIpIdDuIleVK5AAnF9q4FERHZCVtgiIiISHaYwBAREZHsMIEhIiIi2WECQ0RERLLDBIaIiIhkhwkMERERyQ4TGCIiIpIdJjBEREQkO0xgiIiISHaYwBAREZHsMIEhIiIi2WECQ0RERLLDBIaIiIhkhwkMERERyQ4TGCIiIpIdJjBEREQkO0xgiIiISHaYwBAREZHsMIEhIiIi2WECQ0RERLJjcQJz5MgRjBkzBoGBgRAEATt27DAqF0URCxYsQEBAANRqNaKiovDjjz8axVy/fh0TJ06Eu7s7PD09MWXKFJSVlRnFnD59GkOGDIGzszOCgoKwYsUKy78dERERtUkWJzDl5eW455578MEHH9RbvmLFCrz33ntYt24dMjMz4eLigujoaNy8eVOKmThxIs6ePYsDBw5g165dOHLkCKZNmyaVl5SUYMSIEejcuTOysrKwcuVKJCQk4MMPP2zCVyQiIqK2RmHpCTExMYiJiam3TBRFvPvuu5g/fz4eeeQRAMA///lP+Pn5YceOHXj66afx/fffY+/evThx4gQGDhwIAHj//fcxatQopKSkIDAwEJs3b4Zer8eGDRugUqnQp08fZGdnY9WqVUaJDhEREbVPVh0Dc+nSJRQUFCAqKko65uHhgfDwcGRkZAAAMjIy4OnpKSUvABAVFQUHBwdkZmZKMUOHDoVKpZJioqOjkZubi99//73ee1dUVKCkpMToQ0RERG2TVROYgoICAICfn5/RcT8/P6msoKAAvr6+RuUKhQLe3t5GMfVdo+49/iwpKQkeHh7SJygoqPlfiIiIiFqlNvMW0rx581BcXCx9rly5Yu8qERERkY1YNYHx9/cHABQWFhodLywslMr8/f1RVFRkVF5ZWYnr168bxdR3jbr3+DMnJye4u7sbfYiIiKhtsmoC06VLF/j7+yMtLU06VlJSgszMTERERAAAIiIicOPGDWRlZUkxX375JaqrqxEeHi7FHDlyBAaDQYo5cOAAevToAS8vL2tWmYiIiGTI4gSmrKwM2dnZyM7OBlAzcDc7Oxt5eXkQBAEzZszAkiVLsHPnTuTk5GDSpEkIDAzE2LFjAQC9evXCyJEjMXXqVBw/fhxff/01pk+fjqeffhqBgYEAgAkTJkClUmHKlCk4e/YsPvnkE6xZswbx8fFW++JEREQkXxa/Rn3y5EkMGzZM2q9NKmJjY7Fx40a8/vrrKC8vx7Rp03Djxg088MAD2Lt3L5ydnaVzNm/ejOnTp2P48OFwcHDAuHHj8N5770nlHh4e2L9/P+Li4jBgwAB06NABCxYs4CvUZHv6cmBZTSKNN68CKhf71oeIiOplcQITGRkJURQbLBcEAYsWLcKiRYsajPH29saWLVtM3ic0NBRHjx61tHpERETUDrSZt5CIiIio/WACQ0RERLJjcRcSkawZtDXjXBqi19a/3RilBhCEpteLiIgswgSG2hX13/sDJsZwGUm5y/wLc8AvEVGLYheSjWj1lQh+YzeC39gNrb7S3tUhIiJqU9gCQ22fUiNt6l47D43GreFYvfZWy8vsC4BKY14sERG1KCYw1PbVHZuicjG/q0elYbcQEVErxS4kIiIikh22wFhAFEXoDFVmxdYd99LYGBit3rxrEhERUQ0mMBbQGarQe8E+i88buCTNZLlaKIGiZ1NrRURE1P4wgSGqS+UCJBTbuxZERNQIJjBNdHJ+FDQqxwbLr2lLMerzoQCAPY8cgY+JN1902t8Q+fkyAICzsuFrEhERUQ0mME2kUTlCo2r48ekqFXViFSZjUSdW4GyuREREjeJbSERERCQ7TGCI6tAatAjZFIKQTSHQGixYC4mIiFoUExgiIiKSHY6BsYAoioCgBwDoKnWAYGLA7R9xtdum/jWvq9RZq4rUTHXn+dEZqqBR2rEyRETUICYwFrhZdRNuPRcAACI/W2D2eZGfRtqoRmQpXaUOWkPDiafWoDPaVpuIRaVOWqZALYrg8GsiopbDBIbalZjtw02Wi9UKCH90rMZsi4bg0MhK4sFBAIDMqpvQwNUaVSQiIjMwgWmiLx5Ng7e64V9Yukqd1PKS/mQ61Ap1wxczaIGVNasaqx2drVlNIiKiNokJTBOpFWpolBrrxIpizQcwXjmZrMLZ0Rml5xcBAI7OHWZyAsLftWV4dPcIAMCOMXvgpTGRpGqvIWb3GAB/jI8iIqIWwwSG2jxBEABRBQAYkvx1I8F6uP2xLlXUO8ek8+pTdw2rm4YquFijskREZBa+Rk1ERESywxYYavPUSkecWxRtVmzNGlY120fnRppcw+rajUKM+mKZNapIREQWYgJjIxqlBjmxOfauBqGmC8nkWlRG3FD6fTIAwGe8m+n1rkyMpbmNvhxYFliz/ebVmlWviYioydiFRERERLLDBIaIiIhkh11IRNZg0NZ0EzVEr61/uzFKDV+tJyKqBxMYojo0KgUuJ4+2+Dz13/vfmsunMSl3mX9hjpchIqoXu5CIiIhIdtgCQ9RUdWZX1r12HhoTr1xDr73V8jL7AqAyMTNz3VgiIqoXExiipqo7NkXlYn5Xj0rDbiEiomZiFxIRERHJDhMYIiIikh0mMERERCQ7HAND1BJULkBCsb1rQUTUZrAFhoiIiGTH6glMcHAwBEG47RMXFwcAiIyMvK3spZdeMrpGXl4eRo8eDY1GA19fX8yZMweVlZXWrioRERHJlNW7kE6cOIGqqipp/8yZM/jLX/6CJ554Qjo2depULFq0SNrXaG7NiVFVVYXRo0fD398f33zzDfLz8zFp0iQolUosW7bM2tUlIiIiGbJ6AtOxY0ej/eTkZHTr1g0PPvigdEyj0cDf37/e8/fv349z587h4MGD8PPzQ79+/bB48WLMnTsXCQkJUKlU9Z5XUVGBiooKab+kpMQK34aIiIhaI5uOgdHr9fjXv/6F559/HkKdSb82b96MDh06oG/fvpg3bx602luL22VkZCAkJAR+fn7SsejoaJSUlODs2bMN3ispKQkeHh7SJygoyDZfioiIiOzOpm8h7dixAzdu3MDkyZOlYxMmTEDnzp0RGBiI06dPY+7cucjNzcW2bdsAAAUFBUbJCwBpv6CgoMF7zZs3D/Hx8dJ+SUkJkxgiIqI2yqYJzEcffYSYmBgEBgZKx6ZNmyZth4SEICAgAMOHD8fFixfRrVu3Jt/LyckJTk5OzaovERERyYPNupB+/vlnHDx4EC+88ILJuPDwcADAhQsXAAD+/v4oLCw0iqndb2jcDBEREbUvNktgUlNT4evri9GjR5uMy87OBgAEBAQAACIiIpCTk4OioiIp5sCBA3B3d0fv3r1tVV0iIiKSEZt0IVVXVyM1NRWxsbFQKG7d4uLFi9iyZQtGjRoFHx8fnD59GjNnzsTQoUMRGhoKABgxYgR69+6NZ599FitWrEBBQQHmz5+PuLg4dhERERERABslMAcPHkReXh6ef/55o+MqlQoHDx7Eu+++i/LycgQFBWHcuHGYP3++FOPo6Ihdu3bh5ZdfRkREBFxcXBAbG2s0bwwRERG1bzZJYEaMGAFRFG87HhQUhMOHDzd6fufOnbFnzx5bVI2IiIjaAK6FRERERLLDBIaIiIhkhwkMERERyQ4TGCIiIpIdJjBEREQkO0xgiIiISHaYwBAREZHsMIEhIiIi2WECYyv6ciDBo+ajL7d3bYiIiNoUJjBEREQkO0xgiIiISHZsshZSm1V3fSd9OaAwkf/ptfVvNxZLbZJWX4neC/YBAM4tioZGxf/1iIiagz9FLWG4lWio1/Q0TmhMSbnLRhUiIiJqn9iFRERERLLDFpgm0v3tFDQefg0H6LW3Wl5mXwBUGvMurDQzjloVrb4KakWlifLKerfrpa9E7d8CURQhWKF+RERtDROYplJqAJWLebEqC2JJloYsPwSIKrNiBy5JM1muxk1871yzrTNUQePU3NoREbU97EIiIiIi2WELjK2oXICEYnvXgmzIuc5baFlvR0GtUDcYq9VXSi0vJ+cPN/kWkrasBHjPevUkImqLmMAQNZEg3BqdolY6QqNs+H8nXaUObr3eqDnP4RtoVM4NX1jlaLU6EhG1VexCIiIiItlhCwyRFegqdaYDBL3RttbQ8OSFukod8EfrjmjuXENERO0MExgiK4j8NNK6scFBAID0qpsw+f6avhxYFliz/eZVvu1GRO0Gu5CIiIhIdtgCQ9REaoUamRMyzYrVVeqklpf0J9NNvrF0vbgIMbvH1OwYtDWtLA2xZM2tupQaqZuKiEiOmMAQNZEgCNA0YeZktUJt8jyd4tYbSuq/97fNmlvsbiIimWMXEhEREckOW2CIWoBGqUFObI55wXVaZ3SvnYdG49ZwrCVrbtWNJSKSOSYwRK1N3bEpKheuuUVEVA92IREREZHsMIEhIiIi2WEXEpGccdFQImqn2AJDREREssMEhoiIiGSHCQwRERHJDhMYIiIikh0mMERERCQ7TGCIiIhIdqyewCQkJEAQBKNPz549pfKbN28iLi4OPj4+cHV1xbhx41BYWGh0jby8PIwePRoajQa+vr6YM2cOKisrrV1VIiIikimbzAPTp08fHDx48NZNFLduM3PmTOzevRufffYZPDw8MH36dDz22GP4+uuvAQBVVVUYPXo0/P398c033yA/Px+TJk2CUqnEsmXLbFFdIiIikhmbJDAKhQL+/v63HS8uLsZHH32ELVu24KGHHgIApKamolevXjh27Bjuvfde7N+/H+fOncPBgwfh5+eHfv36YfHixZg7dy4SEhKgUqnqvWdFRQUqKiqk/ZKSElt8NSIiImoFbDIG5scff0RgYCC6du2KiRMnIi8vDwCQlZUFg8GAqKgoKbZnz57o1KkTMjIyAAAZGRkICQmBn5+fFBMdHY2SkhKcPXu2wXsmJSXBw8ND+gQFBdniqxEREVErYPUEJjw8HBs3bsTevXuxdu1aXLp0CUOGDEFpaSkKCgqgUqng6elpdI6fnx8KCgoAAAUFBUbJS215bVlD5s2bh+LiYulz5coV634xolZIa9AiZFMIQjaFQGvQ2rs6REQtxupdSDExMdJ2aGgowsPD0blzZ3z66adQq9XWvp3EyckJTk5ONrs+ERERtR42X8zR09MT3bt3x4ULF/CXv/wFer0eN27cMGqFKSwslMbM+Pv74/jx40bXqH1Lqb5xNURtma5SB63B0WR5fdv1qtQBggAAUIsiBKvUkIjIPmyewJSVleHixYt49tlnMWDAACiVSqSlpWHcuHEAgNzcXOTl5SEiIgIAEBERgaVLl6KoqAi+vr4AgAMHDsDd3R29e/e2dXWJWpWY7cPNjo38NLLxoOCasWGZVTehgWsTa0VEZH9WT2Bmz56NMWPGoHPnzrh69SoWLlwIR0dHjB8/Hh4eHpgyZQri4+Ph7e0Nd3d3vPLKK4iIiMC9994LABgxYgR69+6NZ599FitWrEBBQQHmz5+PuLg4dhER/YlYrYDgUHnbNhFRW2f1BOa///0vxo8fj2vXrqFjx4544IEHcOzYMXTs2BEAsHr1ajg4OGDcuHGoqKhAdHQ0/v73v0vnOzo6YteuXXj55ZcREREBFxcXxMbGYtGiRdauKlGr5OzojNLzNX/fs96OglrZcBfSNW0ZRu2oaaX54rF98NE03Kqi011D5PZR1q0sEZGdWD2B2bp1q8lyZ2dnfPDBB/jggw8ajOncuTP27Nlj7aoRyYIgCIBYM9+RWK0CxIYTGIhK422x/nmSAADVt2JFUWxuNYmI7MrmY2CIqOkGLjloOkDQw+2PlTqGLE83mcCohRIo/oi9aaiCi3WqSERkF0xgiORMVKH0+2R714KIqMUxgSFqZdRKR5xbFG1WrFZfiYFL0gAAJ+cPh0bV8P/S124UYtQXZqwnpi8HlgXWbL95FVCxrYaIWh8mMEStjCAIJhORhmhUCpPn6VQmxtIQEckMExii9sigrWlpqY9eW/92Y5QaaaI8IiJbYwJDJGMalQKXk0dbfJ767/0Bc95ESrnL/Iuyu4mIWpBNVqMmIiIisiW2wBC1F0qNtKl77Tw0Grf64/TaWy0vsy8AKk39cX+OJSJqQUxgiNqLuuNTVC7mdfeoNOwWIqJWiV1IREREJDtsgSEiYyoXIKHY3rUgIjKJLTBEREQkO0xgiKhl6MuBBI+aT0Nz0BARmYkJDBEREckOx8AQkXU0NmsvZ/glIitiAkNE1mHJfDCc4ZeImoldSERERCQ7bIEhoqZTampaSMzBGX6JyIqYwBBR0wlC07p3OMMvETUTExgiahmcII+IrIhjYIiIiEh2mMAQERGR7DCBISIjWoMWIZtCELIpBFqDBfO1EBG1II6BIWqHdJU6aA2ODZbVt90YtUINgRPOEVELYQJD1A7FbB9uVlzkp5FmXzNzQiY0ShOvRhMRWRG7kIiIiEh22AJD1E44Ozqj9PwiAEDW21FQK+vvQrqmLcOoHTUtNHvGpsFH49rgNXWVOotaaYiIrIUJDFE7IQgCIKoA1IxX0Sjr/99fp6yStjVKNbuFiKhVYgJD1A5p9VUNlonVSpR+nyxta/WVDcbqKm9dRxRF61WQiKgRTGCI2qGBSw6aGZdmOkDQw61nzebNymq4qJpZMSIiM3EQLxEREckOW2CI2gm10hHnFkU3GqfVV0otLyfnD4dG1fCPiWvaMoz63GpVJCIyGxMYonZCEASTyUh9NCqFyXN0lfW/yUREZGtMYIjIiEalwOXk0fauBhGRSRwDQ0Typi8HEjxqPvpye9eGiFoIW2CIqHXTN7KgZN3yxmLrUmoArt1EJFtWT2CSkpKwbds2nD9/Hmq1Gvfddx+WL1+OHj16SDGRkZE4fPiw0Xkvvvgi1q1bJ+3n5eXh5ZdfxqFDh+Dq6orY2FgkJSVBoWDORdSupNxlm9g3rwIqF8vrQ0StgtWzgcOHDyMuLg6DBg1CZWUl3nzzTYwYMQLnzp2Di8utHxZTp07FokWLpH2N5tZsn1VVVRg9ejT8/f3xzTffID8/H5MmTYJSqcSyZcusXWUiIiKSGasnMHv37jXa37hxI3x9fZGVlYWhQ4dKxzUaDfz9/eu9xv79+3Hu3DkcPHgQfn5+6NevHxYvXoy5c+ciISEBKhVnyyJq05SamhYSc+i1t1peZl8AVCaWPqgbS0SyZvNBvMXFxQAAb29vo+ObN29Ghw4d0LdvX8ybNw9a7a2+64yMDISEhMDPz086Fh0djZKSEpw9e7be+1RUVKCkpMToQ0QyJQg13TtmfeokLCqN+bFEJGs2HVBSXV2NGTNm4P7770ffvn2l4xMmTEDnzp0RGBiI06dPY+7cucjNzcW2bdsAAAUFBUbJCwBpv6CgoN57JSUlITEx0UbfhIiIiFoTmyYwcXFxOHPmDL766iuj49OmTZO2Q0JCEBAQgOHDh+PixYvo1q1bk+41b948xMfHS/slJSUICgpqWsWJSD5ULkBCsb1rQUQtzGZdSNOnT8euXbtw6NAh3HnnnSZjw8PDAQAXLlwAAPj7+6OwsNAopna/oXEzTk5OcHd3N/oQERFR22T1BEYURUyfPh3bt2/Hl19+iS5dujR6TnZ2NgAgICAAABAREYGcnBwUFRVJMQcOHIC7uzt69+5t7SoTERGRzFi9CykuLg5btmzB559/Djc3N2nMioeHB9RqNS5evIgtW7Zg1KhR8PHxwenTpzFz5kwMHToUoaGhAIARI0agd+/eePbZZ7FixQoUFBRg/vz5iIuLg5OTk7WrTERERDJj9RaYtWvXori4GJGRkQgICJA+n3zyCQBApVLh4MGDGDFiBHr27IlZs2Zh3Lhx+M9//iNdw9HREbt27YKjoyMiIiLwzDPPYNKkSUbzxhAREVH7ZfUWGFEUTZYHBQXdNgtvfTp37ow9e/ZYq1pEZGdafSV6L9gHADi3KNrilbGJiOriTxAisgqtvgpqRaWJ8sp6txujVjpC4JpFRPQnTGCIyCqGLD8EiObNkj1wSZrZ12VrDRHVx+Yz8RIRyZK+HEjwqPnoy+1dGyL6E/6zhoiazFlx699AWW9HQa1QNxir1VdKLS8n5w832aqi1Vdh4JKD1qtoffRa88sbi61LqalZCoGIbIoJDBE1Wd2xKWqlIzRK836kaFQK+3cLWbKooyWxb16tmR2YiGyKCQwRtQiNSoHLyaPtXQ0iaiOYwBCRVegqdY2WR34aCQBIfzLdZHeTrrIKEPQAGp+awSJKTU0LiTn02lstL7MvmF7Jum4sEbUIJjBEZBW1yYm1Yt161vz3ZlUUXKBsWqX+TBCa1r2j0rBbiKiVYQJDRFQfrnJN1KoxgSGiJlMr1MickGlWrCVdSNd1ZYjZPtwaVSSiNooJDBE1mSAI0ChNjA2pQ6PUICc2x6xYnaGqOdUionaAE9kRERGR7LAFhohaNXPWWDJ3gry6uMYSkbwxgSGiVm3oyn0Qq02tsaSHW68lAICBS+cDMG89prMJD8PFyUpvNxlVx4wZfs19PbsuzvBLZIQJDBG1aq7dl5gsF6tv/Rhz7Z4MwcG8la6t+np2XZzhl6hFcAwMEbU6dddYIiKqD1tgiKjV0Sg1Zr+ebQmbvZ7NGX6JWhwTGCJqdSx5PdsSNns925IZfjlBHpFVMIEhonZJV6mD1uBo9euqFWq+3UTUApjAEFG7ZKuZfjMnZNqk9YiIjHGkHBFRPeq+3VR3u9HzrLl6NhE1iC0wRNRuODs6o/T8IjOj9XDrWfMKd9kPb8DU/DKCg1563ftmZTVczJuKxjKcX4bICBMYImo3BEEARHOzCxVKv082K1KsbnqdzMb5ZYiMMIEhonZDrXTEuUXRZsVaskTBNW0ZRn1ulSoSkZmYwBBRuyEIgtlrJWlUClxOHm1WrK7S+m8zAeD8MkQmMIEhImqtWmJ+GY6tIZliAkNE1J5xbA3JFF+jJiIiItlhCwwRkRVp9VVQK8xbEdsSaqWj9Wb4tWRsjSXqdjexa6pNuHLlCnbt2oWpU6dCoTAvZWjKOU3BBIaIyIqGrtwHsdrEq9oicGtOGT1g5u/iE2/GwMXJmj+ynQBYOTGqi11TbcL8+fOxd+9eTJw4Ee7u7jY7pymYwBARWVHthHYNEasVEBwqb9tuzOBkvenEqM7Ee6Xn58PUxHt1nU14GC5OSrNiqX05deoU/vnPf2Lt2rVmJyJNOaepmMAQETWTs8L2wwnNSYxuxSabnRjdrIqCC6yUwLBrqs0QRRGzZs1Cr1698MILL9jsnOZgAkNE1EwapQaZEzKtfl2tQYthnw2z+nVtxpLXvpuKXVMtYvfu3UhPT8euXbvMHsfSlHOagwkMEVEzCYJgkxWo1Qq1TRKj67oyaTXu69oyq1/f6vTlUP/ROqIWRXOHDVETVVZWYs6cOXjooYcwatQom53TXExgiIhaKVslRjpDlbT92O4RJmObOmanMRZfNzgIAJD+yAH4aNysUgezu6baWbfU+vXrkZubiy1btpg9wLsp5zQXExgiIpINnYMDtNb6BekgSAmGOuUu81p2GumWEgHoaus350JNEmNlaoXaZklCSUkJFi5ciGeffRZhYWE2O8camMAQEbUzXs4uSH/iG3tXw2zXtWVSS1Ft15fV1Lbs/PxfqEWx/hhRRG0aogVMtqroBAGRne+s2bHR+KXMJw5Bo1Bb/8JKDZYvX47S0lIsXbrU7NOaco41tOoE5oMPPsDKlStRUFCAe+65B++//z4GDx5s72oREcmag4OD9bphWpmmdnlJSUcLcq6uxk0Hh9u2G7XyLqChZKsZrjx7HKtWrcKsWbNw553mPY8rV65YfI61tNoE5pNPPkF8fDzWrVuH8PBwvPvuu4iOjkZubi58fX3tXT0iImohzo7OKD2/yLxgCyYKFBz0jb6eDjQ9KTLZqgM0uWVH11j3kQXXreuthYlwd3fH3LlzzYoHaiats/Qca2m1CcyqVaswdepUPPfccwCAdevWYffu3diwYQPeeOMNO9eOiIhaikalwLnEMVa/bnlFJQYtM2PCP4tmTxahRgUA4AFRCVPBIoCbcAYAOOOm6csKeiiwDIBtWot0l3W4uPnfWLt2LdzczGudqztpnbnnWFOrTGD0ej2ysrIwb9486ZiDgwOioqKQkZFR7zkVFRWoqKiQ9ouLa5aVLykpsVq9SktKUaWrkraVDtYfnEVERC1DJYo4NutBe1fDLNe1ZXh8T1XjgU0giiLyN+ej211d8eSTT5r1e1MURbz22mvo3r272eeYq/ZaYiPdZK0ygfntt99QVVUFPz8/o+N+fn44f/58veckJSUhMTHxtuNBQUE2qWMXdLfJdYmIiOzhIn6Cj4+Pxec15RxzlJaWwsPDo8HyVpnANMW8efMQHx8v7VdXV+P69evw8fGx6utmJSUlCAoKwpUrV2y+zkN7xOdrW3y+tsNna1t8vrZRWVmJwYMH4+LFi8jLyzOZMNQ9595770VAQAB27txp9Ve6RVFEaWkpAgMDTca1ygSmQ4cOcHR0RGFhodHxwsJC+Pv713uOk5MTnJycjI55enraqopwd3fn/0Q2xOdrW3y+tsNna1t8vta1du1a/PTTTwAADw8Ps57t2rVrceHCBXzyySdmJTxNYc51bb8CWROoVCoMGDAAaWlp0rHq6mqkpaUhIiLCjjUjIiJqG2onoBs/frzF50yaNKlFJ62rT6tMYAAgPj4e69evx6ZNm/D999/j5ZdfRnl5ufRWEhERETVd7QR08+fPt/icJUsaf/3c1lplFxIAPPXUU/j111+xYMECFBQUoF+/fti7d+9tA3tbmpOTExYuXHhbdxVZB5+vbfH52g6frW3x+VpX3QnounbtatazteekdfURxMbeUyIiIqI2ZdKkSdi3bx8uXLhg9hwuTTnHllptCwwRERFZ36lTp/Dxxx9bPGmdpefYGltgiIiI2glRFPHQQw+hsLAQp0+fhkLReDtGU85pCa2jFkRERGRzu3fvRnp6Onbt2mV2ItKUc1oCW2CIiIjaAYPBgNDQUAQGBuLgwYNmTUDXlHNaSqt9jbq1+uCDDxAcHAxnZ2eEh4fj+PHj9q6S7CQlJWHQoEFwc3ODr68vxo4di9zcXKOYmzdvIi4uDj4+PnB1dcW4ceNum9iQGpecnAxBEDBjxgzpGJ9t8/zyyy945pln4OPjA7VajZCQEJw8eVIqF0URCxYsQEBAANRqNaKiovDjjz/ascbyUVVVhbfffhtdunSBWq1Gt27dsHjxYqM1cfh8zXfkyBGMGTMGgYGBEAQBr7zyCnJzc5GSkgJBEMx6lmvWrMH58+dx7NgxeHl5YcqUKSgrK7PTN/oTkcy2detWUaVSiRs2bBDPnj0rTp06VfT09BQLCwvtXTVZiY6OFlNTU8UzZ86I2dnZ4qhRo8ROnTqJZWVlUsxLL70kBgUFiWlpaeLJkyfFe++9V7zvvvvsWGv5OX78uBgcHCyGhoaKr732mnScz7bprl+/Lnbu3FmcPHmymJmZKf7000/ivn37xAsXLkgxycnJooeHh7hjxw7xu+++E//617+KXbp0EXU6nR1rLg9Lly4VfXx8xF27domXLl0SP/vsM9HV1VVcs2aNFMPna749e/aIb731lrht2zYRgOju7i7GxsZK5Y09y+LiYlGlUoleXl7isWPHxKNHj4p33XWXOH78eDt9I2NMYCwwePBgMS4uTtqvqqoSAwMDxaSkJDvWSv6KiopEAOLhw4dFURTFGzduiEqlUvzss8+kmO+//14EIGZkZNirmrJSWloq3n333eKBAwfEBx98UEpg+GybZ+7cueIDDzzQYHl1dbXo7+8vrly5Ujp248YN0cnJSfz3v//dElWUtdGjR4vPP/+80bHHHntMnDhxoiiKfL7NAUB88cUXxStXroiiaN6zPH36tAhA3LVrlxTzxRdfiIIgiL/88kvLfoF6sAvJTHq9HllZWYiKipKOOTg4ICoqChkZGXasmfwVFxcDALy9vQEAWVlZMBgMRs+6Z8+e6NSpE5+1meLi4jB69GijZwjw2TbXzp07MXDgQDzxxBPw9fVFWFgY1q9fL5VfunQJBQUFRs/Xw8MD4eHhfL5muO+++5CWloYffvgBAPDdd9/hq6++QkxMDAA+3+YaOXKkNAGdOc/yxIkT8PT0xOjRo6WYqKgoODg4IDMzs2UrX4/WM5y4lfvtt99QVVV120zAfn5+OH/+vJ1qJX/V1dWYMWMG7r//fvTt2xcAUFBQAJVKddtinH5+figoKLBDLeVl69atOHXqFE6cOHFbGZ9t8/z0009Yu3Yt4uPj8eabb+LEiRN49dVXoVKpEBsbKz3D+n5O8Pk27o033kBJSQl69uwJR0dHVFVVYenSpZg4cSIA8PlakTnPsqCgAL6+vkblCoUC3t7ereJ5M4Ehu4qLi8OZM2fw1Vdf2bsqbcKVK1fw2muv4cCBA3B2drZ3ddqc6upqDBw4EMuWLQMAhIWF4cyZM1i3bh1iY2PtXDv5+/TTT7F582Zs2bIFffr0QXZ2NmbMmIHAwEA+X7oNu5DM1KFDBzg6Ot72tkZhYSH8/f3tVCt5mz59Onbt2oVDhw4Zravh7+8PvV6PGzduGMXzWTcuKysLRUVF6N+/PxQKBRQKBQ4fPoz33nsPCoUCfn5+fLbNEBAQgN69exsd69WrF/Ly8gBAeob8OdE0c+bMwRtvvIGnn34aISEhePbZZzFz5kwkJSUB4PO1JnOepb+/P4qKiozKKysrcf369VbxvJnAmEmlUmHAgAFIS0uTjlVXVyMtLQ0RERF2rJn8iKKI6dOnY/v27fjyyy/RpUsXo/IBAwZAqVQaPevc3Fzk5eXxWTdi+PDhyMnJQXZ2tvQZOHAgJk6cKG3z2Tbd/ffff9sr/z/88AM6d+4MAOjSpQv8/f2Nnm9JSQkyMzP5fM2g1Wrh4GD8a8nR0RHV1dUA+HytyZxnGRERgRs3biArK0uK+fLLL1FdXY3w8PAWr/Nt7D2KWE62bt0qOjk5iRs3bhTPnTsnTps2TfT09BQLCgrsXTVZefnll0UPDw8xPT1dzM/Plz5arVaKeemll8ROnTqJX375pXjy5EkxIiJCjIiIsGOt5avuW0iiyGfbHMePHxcVCoW4dOlS8ccffxQ3b94sajQa8V//+pcUk5ycLHp6eoqff/65ePr0afGRRx7ha75mio2NFe+44w7pNept27aJHTp0EF9//XUphs/XfKWlpeK3334rfvvttyIAcdWqVeK3334r/vzzz6IomvcsR44cKYaFhYmZmZniV199Jd599918jVqu3n//fbFTp06iSqUSBw8eLB47dszeVZIdAPV+UlNTpRidTif+7W9/E728vESNRiM++uijYn5+vv0qLWN/TmD4bJvnP//5j9i3b1/RyclJ7Nmzp/jhhx8alVdXV4tvv/226OfnJzo5OYnDhw8Xc3Nz7VRbeSkpKRFfe+01sVOnTqKzs7PYtWtX8a233hIrKiqkGD5f8x06dKjen7W1c8GY8yyvXbsmjh8/XnR1dRXd3d3F5557TiwtLbXDt7kdlxIgIiIi2eEYGCIiIpIdJjBEREQkO0xgiIiISHaYwBAREZHsMIEhIiIi2WECQ0RERLLDBIaIiIhkhwkMERERyQ4TGCIiIpIdJjBEJGvBwcF499137V0NImphTGCIiIhIdrgWEhG1apGRkejbty8A4OOPP4ZSqcTLL7+MRYsWYdiwYTh8+LBRPH+kEbUPbIEholZv06ZNUCgUOH78ONasWYNVq1bhH//4B7Zt24Y777wTixYtQn5+PvLz8+1dVSJqIQp7V4CIqDFBQUFYvXo1BEFAjx49kJOTg9WrV2Pq1KlwdHSEm5sb/P397V1NImpBbIEholbv3nvvhSAI0n5ERAR+/PFHVFVV2bFWRGRPTGCIiIhIdpjAEFGrl5mZabR/7Ngx3H333XB0dIRKpWJLDFE7xASGiFq9vLw8xMfHIzc3F//+97/x/vvv47XXXgNQMw/MkSNH8Msvv+C3336zc02JqKXwNWoiatUiIyPRp08fVFdXY8uWLXB0dMTLL7+MJUuWQBAEHDt2DC+++CJyc3NRUVHB16iJ2gkmMETUqkVGRqJfv36cbZeIjLALiYiIiGSHCQwRERHJDruQiIiISHbYAkNERESywwSGiIiIZIcJDBEREckOExgiIiKSHSYwREREJDtMYIiIiEh2mMAQERGR7DCBISIiItn5/4D3pntXigQGAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ptvals = np.random.exponential(scale=10.0, size=10000) + np.random.exponential(scale=15.0, size=10000)\n",
    "etavals = np.random.normal(scale=1.1, size=10000)\n",
    "\n",
    "dists.fill(\n",
    "    dataset=\"gen2rwt\",\n",
    "    pt=ptvals,\n",
    "    eta=etavals,\n",
    "    weight=corr(ptvals, etavals)\n",
    ")\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "dists[:, :, sum].plot1d(ax=ax)\n",
    "ax.legend(title=\"dataset\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `corr()` can accept also jagged arrays if need be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CMS high-level tools\n",
    "\n",
    "### Applying energy scale transformations with jetmet_tools\n",
    "\n",
    "The `coffea.jetmet_tools` package provides a convenience class [JetTransformer](https://coffeateam.github.io/coffea/api/coffea.jetmet_tools.JetTransformer.html#coffea.jetmet_tools.JetTransformer) which applies specified corrections and computes uncertainties in one call. First we build the desired jet correction stack to apply. This will usually be some set of the various JEC and JER correction text files that depends on the jet cone size (AK4, AK8) and the pileup mitigation algorithm, as well as the data-taking year they are associated with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Fall17_17Nov2017_V32_MC_L2Relative_AK4PFPuppi', 'Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi']\n"
     ]
    }
   ],
   "source": [
    "from coffea.jetmet_tools import FactorizedJetCorrector, JetCorrectionUncertainty\n",
    "from coffea.jetmet_tools import JECStack, CorrectedJetsFactory\n",
    "import awkward as ak\n",
    "import numpy as np\n",
    "\n",
    "ext = extractor()\n",
    "ext.add_weight_sets([\n",
    "    \"* * data/Fall17_17Nov2017_V32_MC_L2Relative_AK4PFPuppi.jec.txt\",\n",
    "    \"* * data/Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi.junc.txt\",\n",
    "])\n",
    "ext.finalize()\n",
    "\n",
    "jec_stack_names = [\n",
    "    \"Fall17_17Nov2017_V32_MC_L2Relative_AK4PFPuppi\",\n",
    "    \"Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi\"\n",
    "]\n",
    "\n",
    "evaluator = ext.make_evaluator()\n",
    "\n",
    "jec_inputs = {name: evaluator[name] for name in jec_stack_names}\n",
    "jec_stack = JECStack(jec_inputs)\n",
    "### more possibilities are available if you send in more pieces of the JEC stack\n",
    "# mc2016_ak8_jxform = JECStack([\"more\", \"names\", \"of\", \"JEC parts\"])\n",
    "\n",
    "print(dir(evaluator))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we prepare some auxilary variables that are used to parameterize the jet energy corrections, such as jet area, mass, and event $\\rho$ (mean pileup energy density), and pass all of these into the `CorrectedJetsFactory`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting columns: {'chHEF', 'jetId', 'pt_gen', 'phi', 'bRegCorr', 'electronIdx2', 'bRegRes', 'mass', 'btagCMVA', 'electronIdx1', 'pt', 'partonFlavour', 'muEF', 'jercCHPUF', 'muonIdx1', 'genJetIdx', 'hadronFlavour', 'btagDeepB', 'btagDeepC', 'neEmEF', 'btagDeepFlavC', 'electronIdx2G', 'genJetIdxG', 'pt_raw', 'btagDeepFlavB', 'qgl', 'muonIdx1G', 'electronIdxG', 'muonSubtrFactor', 'electronIdx1G', 'rho', 'eta', 'nConstituents', 'neHEF', 'cleanmask', 'chEmEF', 'muonIdx2', 'jercCHF', 'btagCSVV2', 'puId', 'muonIdxG', 'rawFactor', 'nMuons', 'mass_raw', 'muonIdx2G', 'area', 'nElectrons'}\n",
      "new columns: {'mass_jec', 'jet_energy_uncertainty_jes', 'JES_jes', 'jet_energy_correction', 'pt_orig', 'pt_jec', 'mass_orig'}\n"
     ]
    }
   ],
   "source": [
    "name_map = jec_stack.blank_name_map\n",
    "name_map['JetPt'] = 'pt'\n",
    "name_map['JetMass'] = 'mass'\n",
    "name_map['JetEta'] = 'eta'\n",
    "name_map['JetA'] = 'area'\n",
    "\n",
    "jets = events.Jet\n",
    "    \n",
    "jets['pt_raw'] = (1 - jets['rawFactor']) * jets['pt']\n",
    "jets['mass_raw'] = (1 - jets['rawFactor']) * jets['mass']\n",
    "jets['pt_gen'] = ak.values_astype(ak.fill_none(jets.matched_gen.pt, 0), np.float32)\n",
    "jets['rho'] = ak.broadcast_arrays(events.fixedGridRhoFastjetAll, jets.pt)[0]\n",
    "name_map['ptGenJet'] = 'pt_gen'\n",
    "name_map['ptRaw'] = 'pt_raw'\n",
    "name_map['massRaw'] = 'mass_raw'\n",
    "name_map['Rho'] = 'rho'\n",
    "    \n",
    "corrector = FactorizedJetCorrector(\n",
    "    Fall17_17Nov2017_V32_MC_L2Relative_AK4PFPuppi=evaluator['Fall17_17Nov2017_V32_MC_L2Relative_AK4PFPuppi'],\n",
    ")\n",
    "uncertainties = JetCorrectionUncertainty(\n",
    "    Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi=evaluator['Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi']\n",
    ")\n",
    "\n",
    "jet_factory = CorrectedJetsFactory(name_map, jec_stack)\n",
    "corrected_jets = jet_factory.build(jets)\n",
    "\n",
    "print('starting columns:', set(ak.fields(jets)))\n",
    "print('new columns:', set(ak.fields(corrected_jets)) - set(ak.fields(jets)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we show that the corrected jets indeed have a different $p_T$ and mass than we started with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "untransformed pt ratios [[1.12, 1.09, 1.2, 1.35, 1.27], [1.03, 1.08, ..., 1, 0.918], ..., [1.13, 0.978]]\n",
      "untransformed mass ratios [[1.12, 1.09, 1.2, 1.35, 1.27], [1.03, 1.08, ..., 1, 0.918], ..., [1.13, 0.978]]\n",
      "transformed pt ratios [[1.2, 1.3, 1.46, 2.09, 2.1], [1.09, 1.29, ..., 1.22, 1.83], ..., [1.37, 1.15]]\n",
      "transformed mass ratios [[1.2, 1.3, 1.46, 2.09, 2.1], [1.09, 1.29, ..., 1.22, 1.83], ..., [1.37, 1.15]]\n",
      "JES UP pt ratio [[1.22, 1.35, 1.56, 2.34, 2.37], [1.1, 1.32, ..., 1.94], ..., [1.41, 1.17]]\n",
      "JES DOWN pt ratio [[1.19, 1.25, 1.35, 1.83, 1.83], [1.08, 1.26, ..., 1.73], ..., [1.33, 1.12]]\n"
     ]
    }
   ],
   "source": [
    "print('untransformed pt ratios', (jets.pt/jets.pt_raw).compute())\n",
    "print('untransformed mass ratios', (jets.mass/jets.mass_raw).compute())\n",
    "\n",
    "print('transformed pt ratios', (corrected_jets.pt/corrected_jets.pt_raw).compute())\n",
    "print('transformed mass ratios', (corrected_jets.mass/corrected_jets.mass_raw).compute())\n",
    "\n",
    "print('JES UP pt ratio', (corrected_jets.JES_jes.up.pt/corrected_jets.pt_raw).compute())\n",
    "print('JES DOWN pt ratio', (corrected_jets.JES_jes.down.pt/corrected_jets.pt_raw).compute())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying CMS b-tagging corrections with btag_tools\n",
    "The `coffea.btag_tools` module provides the high-level utility [BTagScaleFactor](https://coffeateam.github.io/coffea/api/coffea.btag_tools.BTagScaleFactor.html#coffea.btag_tools.BTagScaleFactor) which calculates per-jet weights for b-tagging as well as light flavor mis-tagging efficiencies. Uncertainties can be calculated as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SF: [[1.52, 1.56, 1.59, 1.6, 1.6], [0.969, 1.57, ..., 1.6, 1.6], ..., [1.6, 1.6]]\n",
      "systematic +: [[1.72, 1.77, 1.79, 1.8, 1.8], [1.01, 1.78, ..., 1.8, 1.8], ..., [1.8, 1.8]]\n",
      "systematic -: [[1.31, 1.36, 1.38, 1.4, 1.4], [0.925, 1.37, ..., 1.4, 1.4], ..., [1.4, 1.4]]\n"
     ]
    }
   ],
   "source": [
    "from coffea.btag_tools import BTagScaleFactor\n",
    "\n",
    "btag_sf = BTagScaleFactor(\"data/DeepCSV_102XSF_V1.btag.csv.gz\", \"medium\")\n",
    "\n",
    "print(\"SF:\", btag_sf.eval(\"central\", events.Jet.hadronFlavour, abs(events.Jet.eta), events.Jet.pt).compute())\n",
    "print(\"systematic +:\", btag_sf.eval(\"up\", events.Jet.hadronFlavour, abs(events.Jet.eta), events.Jet.pt).compute())\n",
    "print(\"systematic -:\", btag_sf.eval(\"down\", events.Jet.hadronFlavour, abs(events.Jet.eta), events.Jet.pt).compute())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using correctionlib\n",
    "\n",
    "For the most part, using correctionlib is straightforward. We'll show here how to convert the custom correction we derived earlier (`corr`) into a correctionlib object, and save it in the json format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"><span style=\"font-weight: bold\">CorrectionSet</span> (<span style=\"font-style: italic\">schema v2</span>)\n",
       "my custom corrections\n",
       "📂\n",
       "└── 📈 <span style=\"font-weight: bold\">gen2_to_gen1</span> (v0)\n",
       "    Reweights gen2 to agree with gen1\n",
       "    Node counts: <span style=\"font-weight: bold\">MultiBinning</span>: 1\n",
       "    ╭──────────── ▶ input ─────────────╮ ╭──────────── ▶ input ────────────╮\n",
       "    │ <span style=\"font-weight: bold\">pt</span> (real)                        │ │ <span style=\"font-weight: bold\">eta</span> (real)                      │\n",
       "    │ pt                               │ │ eta                             │\n",
       "    │ Range: [0.0, 100.0), overflow ok │ │ Range: [-3.0, 3.0), overflow ok │\n",
       "    ╰──────────────────────────────────╯ ╰─────────────────────────────────╯\n",
       "    ╭─── ◀ output ───╮\n",
       "    │ <span style=\"font-weight: bold\">out</span> (real)     │\n",
       "    │ <span style=\"font-style: italic\">No description</span> │\n",
       "    ╰────────────────╯\n",
       "</pre>\n"
      ],
      "text/plain": [
       "\u001b[1mCorrectionSet\u001b[0m (\u001b[3mschema v2\u001b[0m)\n",
       "my custom corrections\n",
       "📂\n",
       "└── 📈 \u001b[1mgen2_to_gen1\u001b[0m (v0)\n",
       "    Reweights gen2 to agree with gen1\n",
       "    Node counts: \u001b[1mMultiBinning\u001b[0m: 1\n",
       "    ╭──────────── ▶ input ─────────────╮ ╭──────────── ▶ input ────────────╮\n",
       "    │ \u001b[1mpt\u001b[0m (real)                        │ │ \u001b[1meta\u001b[0m (real)                      │\n",
       "    │ pt                               │ │ eta                             │\n",
       "    │ Range: [0.0, 100.0), overflow ok │ │ Range: [-3.0, 3.0), overflow ok │\n",
       "    ╰──────────────────────────────────╯ ╰─────────────────────────────────╯\n",
       "    ╭─── ◀ output ───╮\n",
       "    │ \u001b[1mout\u001b[0m (real)     │\n",
       "    │ \u001b[3mNo description\u001b[0m │\n",
       "    ╰────────────────╯\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import correctionlib, rich\n",
    "import correctionlib.convert\n",
    "\n",
    "# without a name, the resulting object will fail validation\n",
    "sfhist.name = \"gen2_to_gen1\"\n",
    "sfhist.label = \"out\"\n",
    "clibcorr = correctionlib.convert.from_histogram(sfhist)\n",
    "clibcorr.description = \"Reweights gen2 to agree with gen1\"\n",
    "# set overflow bins behavior (default is to raise an error when out of bounds)\n",
    "clibcorr.data.flow = \"clamp\"\n",
    "\n",
    "cset = correctionlib.schemav2.CorrectionSet(\n",
    "    schema_version=2,\n",
    "    description=\"my custom corrections\",\n",
    "    corrections=[clibcorr],\n",
    ")\n",
    "rich.print(cset)\n",
    "\n",
    "with open(\"data/mycorrections.json\", \"w\") as fout:\n",
    "    fout.write(cset.json(exclude_unset=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use this new correction in a similar way to the original `corr()` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.01495674, 1.40799728, 1.31112463, ..., 0.37951701, 1.16222439,\n",
       "       0.73213844])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ceval = cset.to_evaluator()\n",
    "\n",
    "ceval[\"gen2_to_gen1\"].evaluate(ptvals, etavals)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the time of writing, `correctionlib` does not support jagged arrays. A `correctionlib_wrapper` provided in `coffea.lookup_tools` allows for the processing of jagged array inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>[[1, 0.273, 0.722, 1.02, 1.02],\n",
       " [0.496, 0.439, 0.912, 0.952, 1.02, 0.952, 1.16, 1.02],\n",
       " [1, 0.393, 0.609, 0.516, 1],\n",
       " [0.496, 0.69, 0.952],\n",
       " [0.397, 0.347, 0.722, 0.952, 0.952],\n",
       " [0.778, 0.439, 0.732, 0.935, 0.679, 0.952, 1.02, 1.09],\n",
       " [0.331, 0.519, 0.69, 0.776],\n",
       " [0.69, 0.776, 0.679, 0.952],\n",
       " [0.679],\n",
       " [0.993, 0.668, 0.439, 0.732, 0.776, 0.722, 1.02, 1.16, 1.02],\n",
       " ...,\n",
       " [0.888, 0.935],\n",
       " [0.679, 1.02, 1.02, 1.02],\n",
       " [0.273, 0.443, 0.776, 1.13, 1.09, 0.952, 1.16, 1.16, 1.02],\n",
       " [0.749, 0.722, 0.935],\n",
       " [1.13, 1.09],\n",
       " [0.912, 1, 1.09],\n",
       " [1.01],\n",
       " [0.607, 0.551, 1, 0.679, 1, 1.09],\n",
       " [0.952, 1.16]]\n",
       "--------------------------------------------------------------\n",
       "type: 40 * var * float64</pre>"
      ],
      "text/plain": [
       "<Array [[1, 0.273, 0.722, 1.02, 1.02], ..., [...]] type='40 * var * float64'>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from coffea.lookup_tools.correctionlib_wrapper import correctionlib_wrapper\n",
    "\n",
    "wrap_c = correctionlib_wrapper(ceval[\"gen2_to_gen1\"])\n",
    "wrap_c(events.Jet.pt, events.Jet.eta).compute()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can use the awkward utilities `flatten` and `unflatten` to convert awkward arrays into numpy arrays for evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>[[1, 0.273, 0.722, 1.02, 1.02],\n",
       " [0.496, 0.439, 0.912, 0.952, 1.02, 0.952, 1.16, 1.02],\n",
       " [1, 0.393, 0.609, 0.516, 1],\n",
       " [0.496, 0.69, 0.952],\n",
       " [0.397, 0.347, 0.722, 0.952, 0.952],\n",
       " [0.778, 0.439, 0.732, 0.935, 0.679, 0.952, 1.02, 1.09],\n",
       " [0.331, 0.519, 0.69, 0.776],\n",
       " [0.69, 0.776, 0.679, 0.952],\n",
       " [0.679],\n",
       " [0.993, 0.668, 0.439, 0.732, 0.776, 0.722, 1.02, 1.16, 1.02],\n",
       " ...,\n",
       " [0.888, 0.935],\n",
       " [0.679, 1.02, 1.02, 1.02],\n",
       " [0.273, 0.443, 0.776, 1.13, 1.09, 0.952, 1.16, 1.16, 1.02],\n",
       " [0.749, 0.722, 0.935],\n",
       " [1.13, 1.09],\n",
       " [0.912, 1, 1.09],\n",
       " [1.01],\n",
       " [0.607, 0.551, 1, 0.679, 1, 1.09],\n",
       " [0.952, 1.16]]\n",
       "--------------------------------------------------------------\n",
       "type: 40 * var * float64</pre>"
      ],
      "text/plain": [
       "<Array [[1, 0.273, 0.722, 1.02, 1.02], ..., [...]] type='40 * var * float64'>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def myJetSF(jets):\n",
    "    j, nj = ak.flatten(jets), ak.num(jets)\n",
    "    sf = ceval[\"gen2_to_gen1\"].evaluate(np.array(j.pt), np.array(j.eta))\n",
    "    return ak.unflatten(sf, nj)\n",
    "\n",
    "myJetSF(events.Jet.compute())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}