{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Applying corrections to columnar data\n", "\n", "This is a rendered copy of [applying_corrections.ipynb](https://github.com/scikit-hep/coffea/blob/master/binder/applying_corrections.ipynb). You can optionally run it interactively on [binder at this link](https://mybinder.org/v2/gh/coffeateam/coffea/master?filepath=binder%2Fapplying_corrections.ipynb)\n", "\n", "Here we will show how to apply corrections to columnar data using:\n", "\n", "- the `coffea.lookup_tools` package, which was 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`); while some subtools have been deprecated in favor of correctionlib, you may still come across older formats until everyone has fully transitioned.\n", "- the [correctionlib](https://cms-nanoaod.github.io/correctionlib/) package, which provides an experiment-agnostic serializable data format for common correction functions. This format is the preferred format moving forward, as CMS has standardized on this." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Test data**:\n", "We'll use NanoEvents to construct some test data." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [], "source": [ "import awkward as ak\n", "from coffea.nanoevents import NanoEventsFactory, NanoAODSchema\n", "\n", "NanoAODSchema.warn_missing_crossrefs = False\n", "\n", "fname = \"coffea/tests/samples/nano_dy.root\"\n", "access_log = []\n", "events = NanoEventsFactory.from_root(\n", " {fname: \"Events\"},\n", " schemaclass=NanoAODSchema,\n", " metadata={\"dataset\": \"DYJets\"},\n", " mode=\"virtual\",\n", " access_log=access_log,\n", ").events()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Coffea lookup_tools\n", "\n", "The entrypoint for `coffea.lookup_tools` is the [extractor class](https://coffea-hep.readthedocs.io/en/latest/api/coffea.lookup_tools.extractor.html#coffea.lookup_tools.extractor)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [] }, "outputs": [], "source": [ "from coffea.lookup_tools import extractor" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "~/Dropbox/work/pyhep_dev/coffea/binder/data ~/Dropbox/work/pyhep_dev/coffea/binder\n", "~/Dropbox/work/pyhep_dev/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://coffea-hep.readthedocs.io/en/latest/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": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "available evaluator keys:\n", "\t testSF2d\n", "testSF2d: 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: \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": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Electron eta: [??, ??, ??, ??, ??, ??, ??, ??, ??, ..., ??, ??, ??, ??, ??, ??, ??, ??, ??]\n", "Electron pt: [??, ??, ??, ??, ??, ??, ??, ??, ??, ..., ??, ??, ??, ??, ??, ??, ??, ??, ??]\n", "Scale factor: [[], [0.909], [0.953, 0.972], [0.807, 0.827], [], ..., [], [0.946], [], []]\n" ] } ], "source": [ "print(\"Electron eta:\", events.Electron.eta)\n", "print(\"Electron pt:\", events.Electron.pt)\n", "print(\"Scale factor:\", evaluator[\"testSF2d\"](events.Electron.eta, events.Electron.pt))" ] }, { "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": { "tags": [] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/iason/micromamba/envs/awkward-coffea/lib/python3.13/site-packages/mplhep/utils.py:652: 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": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGwCAYAAAC3qV8qAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPm5JREFUeJzt3Ql4U3W6x/G3lBYo0LLJUi27si8KgnVBtFwqIMqIzAgooAyIgooIIg5iWQQGEESuFx4XwFEQ9V5EBhx2BRVkFdkUQUEc2dxogUjTltzn/XdOTLBLUhKak3w/z3NsknOSnJ5W8uv736JcLpdLAAAAbKREcZ8AAACAvwgwAADAdggwAADAdggwAADAdggwAADAdggwAADAdggwAADAdkpKmDp//rwcPXpUypcvL1FRUcV9OgAAwAc6Pd3p06clMTFRSpQoEXkBRsNLUlJScZ8GAAAogu+//16uuOKKyAswWnmxLkB8fHxxnw4AAPBBRkaGKUBYn+MRF2CsZiMNLwQYAADspbDuH3TiBQAAtkOAAQAAtkOAAQAAthO2fWAAAJEpJydHsrKyivs0kI+YmBiJjo6Wi0WAAQCEzfwhx48fl1OnThX3qaAQFSpUkOrVq1/UPG0EGABAWLDCS9WqVSUuLo5JTEM0ZDocDjl58qS5X6NGjSK/FgEGABAWzUZWeKlcuXJxnw4KUKZMGfNVQ4z+vIranEQnXgCA7Vl9XrTygtBn/Zwupq8SAQYAEDZY+y5yfk4EGAAAYDsEGAAAYDsEGAAAAqx9+/YydOhQrmsQEWAAAChGH330kekTcqnnr0lLS5OWLVuKXRFgAACA7RBgAAC4CGfPnpU+ffpIuXLlzMRszz//vNf+N954Q1q3bi3ly5c3s8/26tXLPZHb4cOH5ZZbbjG3K1asaCox/fr1M/dXrFghN954o5m1Vue2uf322+Wbb75xv67T6ZQhQ4aY9yxdurTUqlVLJk2a5N6vFZ2//vWvctlll0l8fLzceuut8sUXX5h98+fPl7Fjx5r7+p666WN2QoCBF4czW2o/tdxsehsAULARI0bI+vXr5f3335dVq1aZJqEdO3a49+tcJ+PHjzdhYcmSJSa0WCElKSlJ/u///s/c3r9/vxw7dkxmzpzpDkbDhg2Tbdu2ydq1a6VEiRLypz/9Sc6fP2/2v/jii7J06VJ55513zHMXLFggtWvXdr9vjx49TFD617/+Jdu3b5drrrlGUlJS5JdffpG//OUv8sQTT0iTJk3Me+qmj9kJM/ECAFBEZ86ckddee03efPNNEw7U66+/LldccYX7mAceeMB9u27duiZ4XHvttea5WrWpVKmS2aez0mq1xdK9e3ev95o7d66ppuzbt0+aNm0qR44ckSuvvNJUabSCohUYyyeffCJbtmwxAaZUqVLmsWnTppkA9b//+78ycOBA894lS5Y0VSE7ogIDAEARaZOONuW0bdvW/ZgGkgYNGrjva/Wja9euUrNmTdOMdPPNN5vHNYAU5MCBA9KzZ08TerQJyKquWM/r16+f7Ny507zXo48+aqo/Fq32aEDSpicNKtZ26NAhr2YoO6MCAwBAkGgzUGpqqtm0iUcrKBpA9L4Gn4Jo6NGqyiuvvCKJiYmm6UgrL9bzrrnmGhNItIlozZo18uc//1k6dOhgKiwaXrRvjDZnXcizymNnBBgAAIqoXr16EhMTI5s3bzYVFvXrr7/K119/bSotX331lfz8888yefJk099FaZ8WT7Gxse4FKS36HO3XouHlpptucjcLXSg+Pt70XdHt7rvvlttuu830cdFwo6tzaxORZ7+YC9/X8z3thgADAEARabNM//79TUdeba7Rfix/+9vfTIdbpaFGg8KsWbNk0KBBsmfPHtOh15NWWbQPy7Jly6Rz585mtWYdkaSv9/LLL5tKilZtnnrqKa/nTZ8+3ey7+uqrzfu9++67pj+LVli0EpOcnCzdunWTKVOmyFVXXSVHjx6V5cuXm47AOipKg41WcLQZSvvsaPOW1V/GDugDE+ZcLpcZTeTPZvHnOfo+ABCJpk6daqok2uSjwUE71bZq1crs0yYjHZ6s4aJx48amEqOdaT1dfvnlZkizBpRq1aqZodEaSBYtWmT6z2iz0eOPP27ex1P58uVNONEwop2CdXTTBx98YJ6rgUhvt2vXTu6//34TYO655x757rvvzHtYnYS1YqPDuPU833rrLbGTKFeYfvJkZGRIQkKCpKenmxJbpNJw0XjMyqC/z75xqRIXS0EPQPE4d+6cqSbUqVPHzIkC+/68fP38pgIDAABshz+ZI8i20R0kLja60IpN6wlr/3N8SoFVFYczR1pPWBPw8wQAoDAEmAii4cWfZh49lmYhAEAoogkJAADYDgEGAADYDgEGAADYDgEGAADYDgEGAIALRmPWfmq52Twn94TNA8yGDRvMbIO6sJTO9KdLc3vSx/LaPGcQ1OmLL9yvsxN62rVrl5nZUCe40fUjdLZBBJ+OOjo8uYvZGIEEAAibAKMra7Zo0UJeeumlPPcfO3bMa5s7d64JKDplsadx48Z5HffII494zcLXsWNHsz6ETqOs4SctLc2sCQH+OgAA+G/v3r3ms9gqIrzwwgsSUfPAdOrUyWz50YWkPL3//vtmnYW6dev+YQ2HC4+16JLjuly4hh9dBKtJkyZmsSlduGrgwIH+njIAABHP4XCYz+IePXqYtZXsLqh9YE6cOGFWvtSVOi+kTUa60qauoqkVluzs39sZN23aZBagspYYV6mpqWZpcV2mPC+ZmZmmcuO5AQBgB6dPn5bevXtL2bJlzQrTM2bMkPbt28vQoUPdn3HDhw83Cz/qMW3btpWPPvrI/XxdMFJXoV65cqU0atTIrJKtCzVqC4dFF3zUz1td1NFOq04Xy0y8r7/+uqm03HXXXV6PP/roo3LNNddIpUqVZOPGjTJq1ChzkbXCoo4fP24WePJkrZ6p+3SZ8QtNmjTJrOYJAIDStYp/y8rx+2J4dtwtaifeMjHRppnGV8OGDZNPP/1Uli5daj7vxowZIzt27JCWLVua/bpC9b59+8wK1doH9b333jMBZffu3XLllVfmnqvDYVa6fuONN8yK1Pfee68JPdqqEY6CGmC0CUgT5YUrTeoPytK8eXNTaXnwwQdNCClqKtQQ5Pm6WoHRzr8AgMik4aXxmJUX9RrW2nD+2jcu1eeBEFp90T/4Fy5cKCkpKeaxefPmmaCijhw5Yu7rV+ux4cOHy4oVK8zjEydONI9lZWXJnDlzpF69eu7Qo/1Nw1XQAszHH39smnzefvvtQo/VUpg2IR0+fFgaNGhg+sZo85Mn635+/WY0+IRDSQwAEFm+/fZbEz7atGnjfiwhIcF8HiqtsuTk5MhVV13l9bzMzEzTFcMSFxfnDi9Km6JOnjwp4SpoAea1116TVq1amRFLhdEOulruqlq1qrmfnJwsf/vb38wPNCYmxjy2evVq88PMq/ko0kqd/pQ3dcVoAIhE2oyjlRB/6b+rVuVl2+iUIk0poe8dKGfOnJHo6GgzKle/eipXrpz7tvV5adEmLP1sCVcli3IhDx486L5/6NAhE0C0P0vNmjXdzTfvvvuuPP/88394vnbQ3bx5sxmZpP1j9L72hta2Oiuc9OrVy/Rn0c6/I0eOlD179sjMmTNNp6ZwVdRSZ1HLmwAQ7vQD/GLns9LnB3tOLB0ZpOFj69at7s/R9PR0+frrr82AFh3sohUYrabo/GjI5fdPZdu2bSZ8WKx+J3379jW9oJV2MtLU17Nnzz88X5t5dL/O66LlL+2sqwHGs/+Kls5WrVolgwcPNlWcKlWqmA5NDKEGAIQb/WNeP0NHjBhhigHaGvHss8+algkNYdp0pP1J+/TpYwoDGmh+/PFHWbt2relH2qVLF5/eR6cn0Y7A1u0ffvjBFCC0ilO/fn0J+wCjw7oKK0lp0MgvbOjoo88++6zQ99EfivajiUTbRneQuNjogJc3A1nSBAAEjo7CHTRokNx+++0SHx8vTz75pHz//ffuQTDaWXfChAnyxBNPmOBRpUoVue6668zxvjp69KgJPxYdsaTbzTff7DUk2y6CWxdDkWh48TWUXIryJgAg+FUYz+HOOuu9dqWwigHaxKT385supF+/fmbz1K1bN6+Cg87AG059YvjkAwCgmH3++efy1VdfmZFI2v/FGv585513FvephSwCDAAAeSxqe6lpc45OP6Jzo2n/T+1GoU1FyBsBBgCAYqZ9U3SYNEJkLSQAAIBgIMAAAADboQnJhoqrfRYAgFBBBQYAANgOAQYAANgOAQYAAE/OsyJpCbmb3kZIIsAAAADbIcAAABABXnnlFbOadcWKFc3WoUMH2bJli9gVAQYAgAjw0UcfSc+ePeXDDz+UTZs2SVJSknTs2NEsDmlHBBgAAIrZ6dOnpXfv3lK2bFmpUaOGzJgxQ9q3by9Dhw41+zMzM2X48OFy+eWXm2Patm3rtYL0/PnzpUKFCrJy5Upp1KiRlCtXTm677TY5duyY+xhdLPLhhx+Wli1bSsOGDeXVV1+V8+fPy9q1a8WOCDAAgPCkKy9rJ1y/N8fvr6G3i/Iafq76PGzYMPn0009l6dKlsnr1arMO0o4dO9z7hwwZYqomixYtkl27dkmPHj1MQDlw4ID7GIfDYdZTeuONN2TDhg1y5MgRE3ryo8dnZWVJpUqVxI6YyA4AEJ6yHCITEy/uNabVL9rznj4qElvW5+rL66+/LgsXLpSUlBTz2Lx58yQxMffcNYjoff1qPTZ8+HBZsWKFeXzixInmMQ0jc+bMkXr16rlDj7WqdV5GjhxpXk/7wtgRAQaXhMOZLY3HrDS3941LNbMJAwBEvv32WxM+2rRp474cCQkJ0qBBA3N79+7dkpOTI1dddZXX5crMzJTKlSu778fFxbnDi9KmqJMnT+Z5iSdPnmyqOdoMVbp0aVv+GPgUAQCEp5i43EqIv7TZyKq8DD8oEhtXtPcOkDNnzkh0dLRZrVq/eipXrtzvbxkT47UvKipKXHk0ZWkzkwaYNWvWSPPmzcWuCDAAgPAUFeVzM06+NLxc7GsUom7duiZ8bN26VWrWrGkeS09Pl6+//lratWsnV199tanAaDVFh0FfjClTpshzzz1nOvu2bt1a7IwAAwBAMSpfvrz07dtXRowYYTrUVq1aVZ599lkpUaKEqaJo05GOUOrTp488//zzJtD8+OOPZvSQVlC6dPFtcd+///3vMmbMGNPXpnbt2nL8+HF3FcezkmMXjEICAKCYTZ8+XZKTk+X22283nWpvuOEGMxza6p+inXU1wDzxxBOmb0y3bt28Kja+mD17tjidTrn77rtN/xhr0yYlO6ICAwBACFRhdJ4Wy9mzZ2Xs2LEycOBAc1+bmPS+bnnp16+f2TxpyPHsA3P48GEJJwQYBITDmVPI/uw8b/uiTEy0KaMCQLj6/PPP5auvvjIjkbT/izX8+c477yzuUwtZBBgEROsJa/w41r9ZHxl2DeCS0k67aemX/KJrU87+/fslNjZWWrVqZSazq1KlyiU/D7sgwAAAUMy0Y64Ok4bvCDAoMm3a0eqIL7TZyKq8bBudUuhEdtok5U9VBwAQWQgwKDLtl1KUGXX1OczECwC4GAyjBgAAtkOAAQAAtkOAAQAAtkOAAQDAgyPLIc1eb2Y2vY3QRIABAAC2Q4DBJaGjjg5P7mI2RiABwKX3yiuvmNWsK1asaDZdc2nLli0Bf5/27dvL0KFDJdgIMAAARICPPvpIevbsKR9++KFs2rRJkpKSpGPHjvLDDz/k+xxd/DFUEWAAAChmp0+flt69e0vZsmXNCtEzZszwqmRkZmbK8OHD5fLLLzfHtG3b1gQSy/z586VChQqycuVKs4p1uXLl5LbbbpNjx465j9HFIh9++GFp2bKlNGzYUF599VU5f/68rF37+/IutWvXlvHjx5uVr+Pj481ikrp69ZAhQ9zH6DnpPGC6dpMVcvSc1qxZYxaUXL9+vcycOdMco1uwFpEkwAAAwpKuxKydcP3dfsv+zf0aersor+G5CrQvhg0bJp9++qksXbpUVq9ebdZB2rFjh3u/BgitmixatEh27dolPXr0MAHlwIED7mMcDodZT+mNN96QDRs2yJEjR0zoyY8en5WVJZUqVfJ6XF+jRYsWZoHJZ555Rm6++WavsKQBRddosh7bunWreZ3rr7/eBJfk5GQZMGCACU+6aaUnGJiJFwAQljR8tF3Y9qJeo/077Yv0vM29NktcTJzP1ZfXX39dFi5cKCkpKeaxefPmSWJiormtQUTv61frseHDh8uKFSvM4xMnTjSPaYiYM2eO1KtXzx16rFWt8zJy5EjzetoXxtOtt94qTzzxhPu+VoIee+wx+fHHH6VkyZKyb98+E2w0wAwaNMh8vfbaayUuLvf71cUo9Xb16tUlmPyuwGiq69q1q/mmtTS0ZMkSr/1aPrLKRtamKdHTL7/8YkplWp7Sklf//v3lzJkzXsdowtTORqVLlzbpbcqUKUX9HgEACFnffvutCR9t2rRxP5aQkCANGjQwt3fv3i05OTly1VVXmaYha1u/fr1888037udoaLDCi9KmqJMnT+b5npMnTzbVnPfee898znpq3bq11/2mTZuaKo2+n1aGdOHJ22+/3dxX+lVDzqXmdwXm7NmzprT0wAMPyF133ZXnMRpYNBVaSpUq5bVfw4uWlbRMpj+0+++/37SzafpUGRkZpmORpkJNk/rD0/fTsKPHAQBQmDIly5hKSFEqN1bl5aM/f2Rex19FeU5+9A/86Ohos1q1fvVUrlw59+2YmBivfVpAyKspS5uINMBon5XmzZv/Yb/2Z7nwddq1a2cqLfp5rmFFn6f9cvbs2SMbN24ssKkqZAJMp06dzFYQ/QbzKx19+eWXpuylbWZWyps1a5Z07tzZXFSt7GhHI+0UNHfuXFOKatKkiezcuVOmT5+eb4DRC6mbRUMQACDCF5z1sRmnoCBysa9RmLp165rwoZ+LNWvWNI+lp6fL119/bYKDVjy0AqPVFG2ZuBhTpkyR5557znT2vbDSUhDtB6PDsPXzXZ9fokQJc25Tp041n7033HCD+1j93NbzDbagdOLVlFa1alVT/nrooYfk559/du/TTkhaSfG8cFpp0YuxefNm9zF6YfQiWFJTU2X//v3y66+/5vmekyZNMiU3awtWpyEAAAKpfPny0rdvXxkxYoQZ4rx3717TtUI/FzWEadORtlzoyKDFixfLoUOHzPwt+rm3fPlyn9/n73//u+m7osUBHW10/Phxs13YhSMvWnXRvi96bjfeeKP7MS046Oe5Z9VGX1s/z3X00U8//WRGOtkiwGjz0T/+8Q8zLEsvlraNacXGSmN6sTTceNJOQdq+pvusY6pVq+Z1jHXfOuZCo0aNMonV2r7//vtAf2sAAASFtjDo6B3tW6J/1GtFQ4dDW/1TtFuGBhjtXKvFgW7dunlVbHwxe/Zs07qhw6K1f4y1aetHYZo1a2aKDzoE22q20gCjn+0X9n/R5iRt6mrcuLFcdtllpvOxLUYh3XPPPV7fsLaTaacircpYvauDQctaF/a1AQDALlUYrWZ49jcdO3asu9uENjHpfd3y0q9fP7N50pDj2QfGl/lY8jtGq0E6AMeThpm8+thoxUhbUoIt6PPAaNuejhc/ePCgua99Yy7sFZ2dnW0ujNVvRr+eOHHC6xjrfrCHZQEAcKnpnCtvvfWWGVWk879ok5G68847+WEUV4D597//bfrAaJlKaYns1KlTpje1Zd26daaNTGcWtI7R4do6QsmiI5a0bKbrNyCInGdF0hJyN70NABFGO+3u7rvbbMHuwJvXBHLahKQVGB2yrAUABKgJSTv7WNUUpZ2JdISQ9mHRTctb3bt3N5USTZJPPvmk1K9f33TCVdqmp/1kdJY+HSKtIUUn29GmJ2uCnl69epnX0U5MOtGODtPS2f10amUAAMKNjjTy/MMeQajAbNu2zVxo3azpj/X2mDFjTKcdnYDujjvuMG1gGkBatWplUqRn/xRt59N1GLRPjA6f1h7NL7/8snu/jiJatWqVCUf6fO20pK/PHDAAAKBIFRjtbVzQGg86trwwWqmxJq3Lj3b+1eADAICv/F2DCPb9ObGYIwDA9qxZaHWBQoQ+6+d04ezB/mAxx3CnKTfLj/+hnY68bxdGO7pFRRXwumdFJub2cZKnj4rEek9VDQAXQ7sw6Dwl1ihXXRdIJ4FDCK4Q7nCYn5P+vC5cGsEfBJhwp+HFCg7+mlbf92MJJQCKmTXNRn4LGCJ0aHi52GlRCDAAgLCgFRedskNne/echgOhRZuNLqbyYiHARJLhB0ViC5nTQJuNrMpLYcd7HgsAIUI/HAPxAYnQRoCJJBpG/Ol74u/xAABcIoxCAgAAtkMFBoFR2Iglf0c3ObOljJzLvc28DgCACxBgEBj+9IXx4VjtefNl7iry4sg6IlIq4SJODgAQbmhCAgAAtkMFBt60025auu+T1+n8L77wZ3STVl3OZkjczIb8dAAAeSLAoOh0lsuijFLyZXSTM7vIpwUACH80IQEAANshwAAAANshwAAAANshwNiRruyclpC76W0AACIMnXgReqObAAAoBBUYAABgOwQYAABgOwQYAABgO/SBCRUu1++LF5qOuSUDszCiLwsnAgBgMwSYUJHlkC9LP5B7e1qQFlEEACBM0IQEAABshwpMCHI89pXElY0P2MKIXosvAgAQBggwoUiDhq+LJPqyMCIAAGGGJiQAAGA7BBgAAGA7BBgAAGA7BBgAAGA7dOK1IxZGBABEOCowAADAdggwAADAdggwAADAdggwAADAdggwAADAdggwAAAg/APMhg0bpGvXrpKYmChRUVGyZMkS976srCwZOXKkNGvWTMqWLWuO6dOnjxw9etTrNWrXrm2e67lNnjzZ65hdu3bJTTfdJKVLl5akpCSZMmXKxXyfAAAgkgPM2bNnpUWLFvLSSy/9YZ/D4ZAdO3bIM888Y74uXrxY9u/fL3fccccfjh03bpwcO3bMvT3yyCPufRkZGdKxY0epVauWbN++XaZOnSppaWny8ssvF+V7BAAAkT6RXadOncyWl4SEBFm9erXXY//93/8tbdq0kSNHjkjNmjXdj5cvX16qV6+e5+ssWLBAnE6nzJ07V2JjY6VJkyayc+dOmT59ugwcONDfUwYAAGEm6H1g0tPTTRNRhQoVvB7XJqPKlSvL1VdfbSos2dnZ7n2bNm2Sdu3amfBiSU1NNdWcX3/9Nc/3yczMNJUbzw0AAISnoC4lcO7cOdMnpmfPnhIfH+9+/NFHH5VrrrlGKlWqJBs3bpRRo0aZZiStsKjjx49LnTp1vF6rWrVq7n0VK1b8w3tNmjRJxo4dG8xvBwAAhHuA0Q69f/7zn8Xlcsns2bO99g0bNsx9u3nz5qbS8uCDD5oQUqpUqSK9n4Ygz9fVCox2/gUAAOGnZDDDy3fffSfr1q3zqr7kpW3btqYJ6fDhw9KgQQPTN+bEiRNex1j38+s3o8GnqOEHAABEeB8YK7wcOHBA1qxZY/q5FEY76JYoUUKqVq1q7icnJ5vh2vpaFu0crOEmr+YjAAAQWfyuwJw5c0YOHjzovn/o0CETQLQ/S40aNeTuu+82Q6iXLVsmOTk5ps+K0v3aVKQddDdv3iy33HKLGYmk9x9//HG599573eGkV69epj9L//79TR+aPXv2yMyZM2XGjBmB/N4BAECkBJht27aZ8GGx+p307dvXzNWydOlSc79ly5Zez/vwww+lffv2ppln0aJF5lgdOaSddTXAePZf0eHYq1atksGDB0urVq2kSpUqMmbMGIZQAwCAogUYDSHaMTc/Be1TOvros88+K/R9tHPvxx9/7O/pIQI5nNnSeMxKc3vfuFSJiw3q4DoAQAhgLSQAAGA7BBgAAGA7BBgAAGA7dBZAyHM4c0Sc2QXsz87zdmHKxESbZS4AAPZDgEHIu2nKh/KblPbp2NYT1vr8unT4BQD7ogkJAADYDhUYhCRt3rFsH91BJLZsvsdqs5FVedk2OqXAYdTaHNV6wpoAny0A4FIjwCAkefZNMYHEx7ld9FjmgQGA8EcTEgAAsB0CDAAAsB0CDAAAsB0CDAAAsB068cL2tNPu4cldivs0AACXEBUYAABgOwQYAABgOwQY2J/zrEhaQu6mtwEAYY8AAwAAbIcAAwAAbIcAAwAAbIdh1Ah9Tofv+ws9NlvKyLnc2y5XAE4OAFAcCDAIfdPqB+zYOBH5snTubUfWEZFSCRd5cgCA4kATEgAAsB0qMAhNMXEiTx/17VhtNrIqL8MPisRqnSVvjrMZEjezYYBOEgBQXAgwCE1RUSKxZf1/noaXgp7nzL6o0wIAhAaakAAAgO0QYAAAgO0QYAAAgO0QYAAAgO3QiRf2p51209KL+ywAAJcQFRgAAGA7BBgAAGA7BBgAAGA7BBgAAGA7BBgAAGA7BBgAAGA7BBgAABD+AWbDhg3StWtXSUxMlKioKFmyZInXfpfLJWPGjJEaNWpImTJlpEOHDnLgwAGvY3755Rfp3bu3xMfHS4UKFaR///5y5swZr2N27dolN910k5QuXVqSkpJkypQpRf0eAQBApAeYs2fPSosWLeSll17Kc78GjRdffFHmzJkjmzdvlrJly0pqaqqcO3fOfYyGl71798rq1atl2bJlJhQNHDjQvT8jI0M6duwotWrVku3bt8vUqVMlLS1NXn755aJ+nwAAIJJn4u3UqZPZ8qLVlxdeeEFGjx4td955p3nsH//4h1SrVs1Uau655x758ssvZcWKFbJ161Zp3bq1OWbWrFnSuXNnmTZtmqnsLFiwQJxOp8ydO1diY2OlSZMmsnPnTpk+fbpX0AEAAJEpoH1gDh06JMePHzfNRpaEhARp27atbNq0ydzXr9psZIUXpceXKFHCVGysY9q1a2fCi0WrOPv375dff/01z/fOzMw0lRvPDQAAhKeABhgNL0orLp70vrVPv1atWtVrf8mSJaVSpUpex+T1Gp7vcaFJkyaZsGRt2m8GuCjOsyJpCbmb3gYAhIywGYU0atQoSU9Pd2/ff/99cZ8SAACwQ4CpXr26+XrixAmvx/W+tU+/njx50mt/dna2GZnkeUxer+H5HhcqVaqUGdXkuQEAgPAU0ABTp04dEzDWrl3rfkz7omjfluTkZHNfv546dcqMLrKsW7dOzp8/b/rKWMfoyKSsrCz3MTpiqUGDBlKxYsVAnjIAAIiEUUg6X8vBgwe9Ou7qCCHtw1KzZk0ZOnSoTJgwQa688koTaJ555hkzsqhbt27m+EaNGsltt90mAwYMMEOtNaQMGTLEjFDS41SvXr1k7NixZn6YkSNHyp49e2TmzJkyY8aMQH7viHRZDhFnAf8LOB153/ZFTJxIVFTRzw0AENgAs23bNrnlllvc94cNG2a+9u3bV+bPny9PPvmkmStGhztrpeXGG280w6Z1QjqLDpPW0JKSkmJGH3Xv3t3MHWPRTrirVq2SwYMHS6tWraRKlSpmcjyGUCOQ4mY29P3gafX9e/Gnj4rElvX7nAAAvoly6eQtYUibrjQIaYdeO/SHcZxJl7hpNXNvDz8iceUSivuUwpLndQ4qAgwABPXz2+8KDGBrMXHS6Nxcc3P76A4SF1tIE5JVeRl+UCQ2ruDX9jweABBUBBhElqgo+U3+05ypTTwFBRhPGl5oEgKAkBE288AAAIDIQYABAAC2Q4ABAAC2Qx8YID/a5yUtnesDACGICgwAALAdAgwAALAdAgwAALAdAgwAALAdAkwQOZzZUvup5WbT2wAAIDAIMAAAwHYIMAAAwHYIMAAAwHaYyM4PLpdLfsvK8fl4z34vhfWBcThzpJC1jgEAwH8QYPyg4aXxmJVSFK0nrC1wfxk5J1/+Z5FkAABQMJqQgHwwigwAQhcVmCLaNrqDxMVGF/oBaFVeto1OkbjYAi6386zItNybZWIKfl0AACIdAaaINLwUGEj+cHzJQo7/fV9UVFRRTwsAgIhAExIAALAdKjCIWDryK1CjyAxntnskmY5Yo44GAMFDgAkibTI6PLlLMN8CF6H1hDUBG0V24UgyHbEWV+pizg4AUBCakAAAgO1QgUFE0RFe+8al+nSsX6PI9PgzGSIvBuQ0AQCFIMAgougIL39Gj/k+ikxEChlWDwAIHJqQAACA7RBgAACA7dCEBOSDUWQAELqowAAAANshwAAAANshwAAAANshwADFQVcfT0vI3fQ2AMAvBBgAAGA7BBgAAGA7BBgAAGA7zAMDBEOWQ8RZwP9eTkfetwsTE6frIVzcuQFAGCDAAEEQN7Oh7wdPq+/7sU8fFYktW6RzAoBwEvAmpNq1a5sF8y7cBg8ebPa3b9/+D/sGDRrk9RpHjhyRLl26SFxcnFStWlVGjBgh2dnZgT5VAABgUwGvwGzdulVycnLc9/fs2SP/9V//JT169HA/NmDAABk3bpz7vgYViz5Xw0v16tVl48aNcuzYMenTp4/ExMTIxIkTA326QODExEmjc3PNze2jOxS8erU2G1mVl+EHRWLjfDsWABCcAHPZZZd53Z88ebLUq1dPbr75Zq/AogElL6tWrZJ9+/bJmjVrpFq1atKyZUsZP368jBw5UtLS0iQ2NjbP52VmZprNkpGREbDvCfBJVJT8JqVzb2szT0EBxpOGF5qFACB0RiE5nU5588035YEHHjBNRZYFCxZIlSpVpGnTpjJq1ChxOH7vxLhp0yZp1qyZCS+W1NRUE0j27t2b73tNmjRJEhIS3FtSUlIQvzMAABC2nXiXLFkip06dkn79+rkf69Wrl9SqVUsSExNl165dprKyf/9+Wbx4sdl//Phxr/CirPu6Lz8ahIYNG+a+r4GHEAMAQHgKaoB57bXXpFOnTiasWAYOHOi+rZWWGjVqSEpKinzzzTemqamoSpUqZTbAFrTJKC29uM8CAGwraE1I3333nenH8te//rXA49q2bWu+Hjx40HzVvjEnTpzwOsa6n1+/GQAAEFmCFmDmzZtnhkDriKKC7Ny503zVSoxKTk6W3bt3y8mTJ93HrF69WuLj46Vx48bBOl0AABDpTUjnz583AaZv375SsuTvb6HNRAsXLpTOnTtL5cqVTR+Yxx9/XNq1ayfNmzc3x3Ts2NEElfvuu0+mTJli+r2MHj3azCNDExEAAAhagNGmI52MTkcfedIh0LrvhRdekLNnz5pOtt27dzcBxRIdHS3Lli2Thx56yFRjypYta4KQ57wxAAAgsgUlwGgVxeVy/eFxDSzr168v9Pk6SumDDz4IxqkBAIAwwGrUAADAdggwAADAdggwAADAdggwAADAdggwweQ8K5KWkLvpbQAAEBAEGAAAYDsEGAAAYDsEGAAAYDtBXY067LhcUkbO5d42fVoKuXxOR963CzsWAAAUiADjjyyHfFn6P8sjTPPrmSLT6vv5BAAAkB+akAAAgO1QgSkix2NfSVzZ+MKbhazKy/CDIrFxvr14jI/HAQAQoQgwRaUhI7as78fH+nk8AADIF01IAADAdqjABJNWXNLSg/oWAABEIiowQDFwOLOl9lPLzaa3AQD+IcAAAADbIcAAAADboQ8MEAQOZ04h+39vNiq0CcmZLdbAepfLJVEFHntWZGJi7u2njzLyDUDYIsAAQdB6who/jl1b4H5dvuLL0rm3f8vKkbhSF3t2AGB/NCEBAADboQIDBEiZmGjZNy7Vp2O12ciqvGwbnSJxsfn/r+g4kyHy4n/uZDlEnCUDs4BoXpMzRhXYQAUAIYMAAwRIVFRUgUEkP/qcAp8XG/37sTMbBm8BUfrMALARmpAAAIDtUIEBQl1MnDQ6N9fc3D66Q8HVGn8XEPU8HgBshAADFAMNIYcnd/Ht4Kgo+U1K/748ha/NVCwgCiCM0YQEAABshwADAABshyYkIJywAjqACEEFBgAA2A4BBgAA2A4BBgAA2A4BBgAA2A4BBgAA2A4BBgAA2A4BBgAA2E7AA0xaWppZlddza9jw9xV0z507J4MHD5bKlStLuXLlpHv37nLixAmv1zhy5Ih06dJF4uLipGrVqjJixAjJzs4O9KkCAACbCspEdk2aNJE1a9b8/iYlf3+bxx9/XJYvXy7vvvuuJCQkyJAhQ+Suu+6STz/91OzPyckx4aV69eqyceNGOXbsmPTp00diYmJk4sSJwThdAABgM0EJMBpYNIBcKD09XV577TVZuHCh3HrrreaxefPmSaNGjeSzzz6T6667TlatWiX79u0zAahatWrSsmVLGT9+vIwcOdJUd2JjY4NxygAAINL7wBw4cEASExOlbt260rt3b9MkpLZv3y5ZWVnSoUMH97HavFSzZk3ZtGmTua9fmzVrZsKLJTU1VTIyMmTv3r35vmdmZqY5xnMDAADhKeABpm3btjJ//nxZsWKFzJ49Ww4dOiQ33XSTnD59Wo4fP24qKBUqVPB6joYV3af0q2d4sfZb+/IzadIk0yRlbUlJSYH+1gAAQLg2IXXq1Ml9u3nz5ibQ1KpVS9555x0pU6aMBMuoUaNk2LBh7vtagSHEAAAQnoI+jFqrLVdddZUcPHjQ9ItxOp1y6tQpr2N0FJLVZ0a/XjgqybqfV78aS6lSpSQ+Pt5rAwAA4SnoAebMmTPyzTffSI0aNaRVq1ZmNNHatWvd+/fv32/6yCQnJ5v7+nX37t1y8uRJ9zGrV682gaRx48bBPl0A+XGeFUlLyN30NgCEUxPS8OHDpWvXrqbZ6OjRo/Lss89KdHS09OzZ0/RN6d+/v2nqqVSpkgkljzzyiAktOgJJdezY0QSV++67T6ZMmWL6vYwePdrMHaNVFgAAgIAHmH//+98mrPz8889y2WWXyY033miGSOttNWPGDClRooSZwE5HDukIo//5n/9xP1/DzrJly+Shhx4ywaZs2bLSt29fGTduHD8tIJicDt/3F3asp5g4kaioop8XAOQhyuVyuSQMaSderfjo3DOB6g/jOJMucdNq5t4efkTiyiUE5HWBAn/vnNnSeMxKc3vfuFSJiw3g3x3aFDQxMbg/gKePisSWDe57AIi4z2/WQgIAALYTlJl4AdiENu9ohcQX2mw0rX7u7eEHRWLjfDsWAIKAAANEMu2bUpTmHQ0vNAsBKEY0IQFh1l+m9lPLzaa3ASBcUYEB4ButuKSlc7UAhAQCDGAjDmdOIft/r7r4W4EpExMtUQx3BmATBBjARlpPWOPHsb/PeO2LgA/RBoAgog8MAACwHf7cAkKcNu1odcQX2mxkVV62jU4ptKKiTVL+VHUAIFQQYIAQp/1SitK0o8+hSQhAuKIJCQAA2A4VGCCMaMXl8OQuxX0aABB0VGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAAIDtEGAAFD/nWZG0hNxNbwNAIZiJF0BwOR3+HePL8ZaYOF0sqmjnBcDWCDAAgmta/eAd//RRkdiyfp8SAPujCQkAANgOFRgAgadNO1od8ZU2G1mVl+EHRWLjfDsWQMQiwADwicOZLY3HrDS3941LNStf50v7pRS1aUfDC81CAApBExIAALAdKjAADIczp9AKTF63C1MmJlqiChsppBWXtHR+EgB8RoABYLSesMbnK9F6wlqfjy20uQkAioAmJAAAYDv8WQREMG3e0QqJL7TZyKq8bBudUmBVRZuj/KnoAIC/CDBABNO+KUVp3tHn0CwEoDgRYAD4RAPL4clduFoAQgJ9YAAAgO0QYAAAgO0QYAAAgO0QYAAAgO0EPMBMmjRJrr32WilfvrxUrVpVunXrJvv37/c6pn379mb0g+c2aNAgr2OOHDkiXbp0kbi4OPM6I0aMkOxs32f/BAAA4Svgo5DWr18vgwcPNiFGA8fTTz8tHTt2lH379knZsr8v7jZgwAAZN26c+74GFUtOTo4JL9WrV5eNGzfKsWPHpE+fPhITEyMTJ04M9CkDAIBIDzArVqzwuj9//nxTQdm+fbu0a9fOK7BoQMnLqlWrTOBZs2aNVKtWTVq2bCnjx4+XkSNHSlpamsTGxv7hOZmZmWazZGRkBPT7AmBDzrMiExNzbz99lFWugTAS9D4w6em5C7RVqlTJ6/EFCxZIlSpVpGnTpjJq1ChxOBzufZs2bZJmzZqZ8GJJTU01oWTv3r35Nl0lJCS4t6SkpKB9TwBChNORG1Ly3Rx+HOuxuVzF+V0BKO6J7M6fPy9Dhw6VG264wQQVS69evaRWrVqSmJgou3btMpUV7SezePFis//48eNe4UVZ93VfXjQEDRs2zH1fww4hBghz0+oH51iqNUBkBxjtC7Nnzx755JNPvB4fOHCg+7ZWWmrUqCEpKSnyzTffSL169Yr0XqVKlTIbAAAIf0ELMEOGDJFly5bJhg0b5Iorrijw2LZt25qvBw8eNAFG+8Zs2bLF65gTJ06Yr/n1mwEQIWLiciskvtBmI6vyMvygSGycb8cCiLw+MC6Xy4SX9957T9atWyd16tQp9Dk7d+40X7USo5KTk2X37t1y8uRJ9zGrV6+W+Ph4ady4caBPGYCdREXldsb1afMILHrb12MBRF4FRpuNFi5cKO+//76ZC8bqs6Ida8uUKWOaiXR/586dpXLlyqYPzOOPP25GKDVv3twcq8OuNajcd999MmXKFPMao0ePNq9NMxEAn2kwScsdSAAgvAS8AjN79mwz8kgnq9OKirW9/fbbZr8Ogdbh0RpSGjZsKE888YR0795d/vnPf7pfIzo62jQ/6Vetxtx7771mHhjPeWMAAEDkKhmMJqSC6MggneyuMDpK6YMPPgjgmQEAgHDBWkgAAMB2CDAAAMB2CDAAip3DmS21n1puNr0NAMU6kR0AOJw5hV4Ez9DiT4ApExNtVrMvFqyzBBQrAgyAoGo9YY2fx6/1+dh941IlLjYI/4x5rqHkyzG+HO85EV9xhS4gjBBgAOBC/s7IyzpLwCVHgAEQcNq0o9URX2mzkVV52TY6pcCqijZJ+VvVARB+CDAAAk77pfjTtKPHHp7cxT5rLCnWWQKKFQEGQNjSyk7jMSt96y9jrbFUFNY6SwAuGQIMgLAd4RTU0U2sswQUKwIMANvypy9MSIxu8mXEkj9NU54Y3YQIQ4ABgEvJnxFLjG4C8kWAARC2I5wY3QSELwIMgLAd4RQSo5v8HeHE6CbAJwQYAAi2oo5wYnQTkC8CDACEEkY3AT4hwABAEReg9HX24JBZgBIIIwQYAAiHBSiBCMP/RQAQDnxdQZs5ZhAmCDAAEA4LULKCNiIMAQYA7LoAJRDBCDAAYFesoI0IRoABgEiZX4Yh2ggjBBgAQPEsQKlYhBJFRIABABTPApRKl1goyizFiHglIv4KAAAA26ECAwAhNMtvyMzwG6wFKC88HigiAgwAXEL+zAdTrDP8+tNB+GI6B/syAV9R0Lcm7BFgAADFJ1iVGPrWhD0CDACE2Cy/vrokM/wCIYoAAwAhNstvsFbQLirb9K+hb01EIcAAQBgIZiXGtv1rENYIMACA8BOszsGKDsIhgQADADYVrL41F/avsVXzlCWYw7T9acryZ2ZigpFfCDAAYFOXom9NsJunto3uIHGx0YF5MWe2+LGIQejNTBzokVMul0hWECtRxRy6CDAAgGIT2HDkkjIy19zaboJRAD/iLkUH4UA3ezkvwTkX43D1kA4wL730kkydOlWOHz8uLVq0kFmzZkmbNm2K+7QAIOxdquapwIqS36R07nuYrwGq7KiYWJHhRwLe7OVynpWoaVfm3mF24vAIMG+//bYMGzZM5syZI23btpUXXnhBUlNTZf/+/VK1atXiPj0ACGvBbJ66FPPi2KXZy+HMlCoSfD8/tFfKlC0fmBfLckjczIbmpsvlkuJpQArhADN9+nQZMGCA3H///ea+Bpnly5fL3Llz5amnnirWc3NERclv2b+JZMUU63kAgG0F4VMvqkSOSJTTt4Nd+p/Y/9xx+nU+rZ/7IICvrc1ec3zu0nJO4s3t0pLhV9eT317YUvBF9+Ocy8g52V46SuJcLvktK0fiSkmxCMkA43Q6Zfv27TJq1Cj3YyVKlJAOHTrIpk2b8nxOZmam2Szp6bnzBmRkZATsvBxnMiQ70yXX1bpcZNGtAXtdAEBgxNUK3yupmaKonZTjAnwurSVRPvvu3+LIyJDs84FNo9bntlZ3bBdgfvrpJ8nJyZFq1ap5Pa73v/rqqzyfM2nSJBk7duwfHk9KSgrCGX4ZhNcEAMA+EvQ/k2sG7fVPnz4tCQnmXewTYIpCqzXaZ8Zy/vx5+eWXX6Ry5coBnWNAk6GGou+//17i43NLeQgOrvWlwXXmOocTfp/tf5218qLhJTExscDjQjLAVKlSRaKjo+XEiRNej+v96tWr5/mcUqVKmc1ThQoVgnaO+gMjwFwaXGuuczjh95nrHE7ig/RZWFDlxVJCQlBsbKy0atVK1q5d61VR0fvJycnFem4AAKD4hWQFRmlzUN++faV169Zm7hcdRn327Fn3qCQAABC5QjbA/OUvf5Eff/xRxowZYyaya9mypaxYseIPHXsvNW2mevbZZ//QXAWutV3xO811Dif8PkfOdY5yFTZOCQAAIMSEZB8YAACAghBgAACA7RBgAACA7RBgAACA7RBg/PTSSy9J7dq1pXTp0maV7C1bdIEsFJUuAXHttddK+fLlzSrj3bp1MyuOezp37pwMHjzYzKpcrlw56d69+x8mOYR/Jk+ebGaoHjp0KNc5wH744Qe59957ze9rmTJlpFmzZrJt2zb3fh03oaMra9SoYfbrGm8HDhwI9GmENV1q5plnnpE6deqYa1ivXj0ZP36819o5XOei2bBhg3Tt2tXMgqv/RixZssRrvy/XVWfB7927t5ngTieU7d+/v5w5c0YCTkchwTeLFi1yxcbGuubOnevau3eva8CAAa4KFSq4Tpw4wSUsotTUVNe8efNce/bsce3cudPVuXNnV82aNV1nzpxxHzNo0CBXUlKSa+3ata5t27a5rrvuOtf111/PNS+iLVu2uGrXru1q3ry567HHHuM6B9Avv/ziqlWrlqtfv36uzZs3u7799lvXypUrXQcPHnQfM3nyZFdCQoJryZIlri+++MJ1xx13uOrUqeP67bffAnkqYe25555zVa5c2bVs2TLXoUOHXO+++66rXLlyrpkzZ7qP4ToXzQcffOD629/+5lq8eLGmQdd7773ntd+X63rbbbe5WrRo4frss89cH3/8sat+/fqunj17ugKNAOOHNm3auAYPHuy+n5OT40pMTHRNmjQp4D+YSHXy5EnzP8369evN/VOnTrliYmLMP1CWL7/80hyzadOmYjxTezp9+rTryiuvdK1evdp18803uwMM1zkwRo4c6brxxhvz3X/+/HlX9erVXVOnTnU/pte+VKlSrrfeeitAZxH+unTp4nrggQe8HrvrrrtcvXv3Nre5zoFxYYDx5bru27fPPG/r1q3uY/71r3+5oqKiXD/88IMrkGhC8pHT6ZTt27ebcpmlRIkS5v6mTZsCXxqLUOnp6eZrpUqVzFe95llZWV7XvWHDhlKzZk2uexFoU1yXLl28rifXOXCWLl1qZg/v0aOHaRK9+uqr5ZVXXnHvP3TokJmY0/P665ov2hzNvyO+u/76683SMl9//bW5/8UXX8gnn3winTp14joHkS+/v/pVm430/wOLHq+fl5s3b46MmXhDzU8//WTaXS+cCVjvf/XVV8V2XuFE17vSPhk33HCDNG3a1Dym/7Po2lgXLsyp1133wXeLFi2SHTt2yNatW/+wj+scGN9++63Mnj3bLIXy9NNPm2v96KOPmt9hXRrF+p3N698Rfp9999RTT5nVkPWPGV34V/9tfu6550y/C+v3mesceL5cV/2q4d1TyZIlzR+lgf4dJ8AgpKoDe/bsMX9JIbB0yfvHHntMVq9ebTqgI3ghXP/ynDhxormvFRj9nZ4zZ44JMAiMd955RxYsWCALFy6UJk2ayM6dO80fP9rxlOscOWhC8lGVKlVM0r9w9Iver169ejB+NhFlyJAhsmzZMvnwww/liiuucD+u11ab706dOuV1PNfdP9oUd/LkSbnmmmvMX0O6rV+/Xl588UVzW/+C4jpfPB2Z0bhxY6/HGjVqJEeOHDG3rX8r+Hfk4owYMcJUYe655x4zyuu+++6Txx9/3Ixq5DoHjy+/v/pV/63xlJ2dbUYmBfqzkgDjIy0Bt2rVyrS7ev61pfeTk5MD+kOJJNpPTMPLe++9J+vWrTPDIj3pNY+JifG67jrMWj8QuO6+S0lJkd27d5u/VK1NKwVacrduc50vnjZ/XjgNgPbTqFWrlrmtv9/6j7jn77M2hWjfAH6ffedwOEyfCk/6B6b+m8x1Dh5ffn/1q/7BqX80WfTfdv3ZaF+ZgApol+AIGEatva3nz59veloPHDjQDKM+fvx4cZ+abT300ENmSN5HH33kOnbsmHtzOBxew6h1aPW6devMMOrk5GSz4eJ4jkLiOgduiHrJkiXNMN8DBw64FixY4IqLi3O9+eabXsNQ9d+N999/37Vr1y7XnXfeyTBqP/Xt29d1+eWXu4dR65DfKlWquJ588kmucwBGKn7++edm04gwffp0c/u7777z+fdXh1FfffXVZiqBTz75xIx8ZBh1CJg1a5b5MNX5YHRYtY5zR9Hp/yB5bTo3jEX/x3j44YddFStWNB8Gf/rTn0zIQWADDNc5MP75z3+6mjZtav7Yadiwoevll1/22q9DUZ955hlXtWrVzDEpKSmu/fv3B+jdI0NGRob53dV/i0uXLu2qW7eumbskMzPTfQzXuWg+/PDDPP9N1tDo63X9+eefTWDRuXni4+Nd999/vwlGgRal/wlsTQcAACC46AMDAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADAABshwADwNZq164tL7zwQnGfBoBLjAADAABsh7WQAIS09u3bS9OmTc3tN954Q2JiYuShhx6ScePGyS233CLr16/3Op7l3YDIQAUGQMh7/fXXpWTJkrJlyxaZOXOmTJ8+XV599VVZvHixXHHFFSbMHDt2zGwAIkPJ4j4BAChMUlKSzJgxQ6KioqRBgwaye/duc3/AgAESHR0t5cuXl+rVq3MhgQhCBQZAyLvuuutMeLEkJyfLgQMHJCcnp1jPC0DxIcAAAADbIcAACHmbN2/2uv/ZZ5/JlVdeaZqPYmNjqcQAEYgAAyDkHTlyRIYNGyb79++Xt956S2bNmiWPPfaYex6YDRs2yA8//CA//fRTcZ8qgEuEYdQAQn4YdZMmTeT8+fOycOFCU3XRYdQTJkww/WK0GvPggw+acJOZmckwaiBCEGAAhHyAadmyJbPtAvBCExIAALAdAgwAALAdmpAAAIDtUIEBAAC2Q4ABAAC2Q4ABAAC2Q4ABAAC2Q4ABAAC2Q4ABAAC2Q4ABAAC2Q4ABAABiN/8P+uFzuhuU3ckAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "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.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": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 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=, cbar=, text=[])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAApkAAAG2CAYAAAA0vZ0sAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMP9JREFUeJzt3Q14FNW5wPF3E0hChASRjxAMgqJ8iHwIBYN6BUWRy4NSW0upFyIq3lJoKfELUEFEidWCUI0iKCpeEMQKVYtwKRC41AASoaIiCqJJIwEihUCARHbnPufYrFlI8CzMZHY2/5/PMZnJ7M6wk03evO/58FmWZQkAAABgoxg7nwwAAAAgyAQAAIAjyGQCAADAdgSZAAAAsB1BJgAAAGxHkAkAAADbEWQCAADAdgSZAAAAsB1BJgAAAGxHkAkAAIDaG2Q+//zz0qlTJ0lKStItPT1d3nvvvdM+ZvHixdKuXTtJSEiQyy67TJYtW1Zj1wsAAFCbeSbIPP/88+WJJ56QvLw82bx5s1x77bVy8803yyeffFLl8e+//74MGTJE7rzzTtmyZYsMGjRIt48//rjGrx0AAKC28VmWZYlHNWrUSJ566ikdSJ5s8ODBUlpaKu+++25w3xVXXCFdunSRWbNm1fCVAgAA1C51xIP8fr8uhasgUpXNq5KbmyuZmZkh+/r16ydLly497XOXlZXpViEQCMiBAwfkvPPOE5/PZ9O/AAAAOEnl0A4fPiypqakSE2NP4fb48eNSXl4ukSAuLk53B4xkngoyt23bpoNKdZPr168vS5YskQ4dOlR5bFFRkTRr1ixkn9pW+08nKytLJk+ebOt1AwAAdxQUFOgud2dLxR6tL6gvRfv8EglSUlJk9+7dER1oeirIbNu2rWzdulUOHTokb775pmRkZMjatWurDTTPxPjx40MyoOpcLVu2lEfWXCkJ9Z17uZ5bfb04rWf3zx19/o/2porTaqJzx7H9iY6fI+kz5996JW2c/UF47iex4rgauN9HU5w/R/KXzv9DGuw+6ujz+2rgzfevduc4fo7Gq/MdP0dgDlWvSHDiaLmsu3WuNGjQwJbnUxlMFWB+nddKkhq4O6Sl5HBALuj2lb4mgkwbU8Nt2rTRn3fr1k0++OADmTlzprzwwgtVRvh79+4N2ae21f7TiY+P1+1kKsB0MsiMqef8XyJ1z4lz9PljE0993bwYZNbEvYiNdz7IjKnnbJAZGxcdQWas89+2Uqeu8/+QOnUCng8yY+Ocf+/ViXH+hgfOIciMJHZ3davfwKebmwLije8xz4wur4rqL1m5/2Rlqqy+atWqkH0rV66stg8nAADAj/FbgYhoXuCZcrkqY/fv31+XrlVH3gULFkhOTo6sWLFCf33YsGHSokUL3adSGTNmjFxzzTUybdo0GTBggCxcuFBPfTR79myX/yUAAMCrAmLp5vY1eIFngsx9+/bpQHLPnj2SnJysJ2ZXAeb113/flzE/Pz9k9FivXr10IPrQQw/JhAkT5OKLL9Yjyzt27OjivwIAAKB28EyQ+dJLL5326yqrebJbb71VNwAAADsE9H/uCrh+BVEWZAIAALjNb1m6uX0NXuDpgT8AAACITGQyAQAADDHwxxxBJgAAQBhBpp/R5UYolwMAAMB2ZDIBAAAMUS43R5AJAABgiNHl5iiXAwAAwHZkMgEAAAypadDdngo9IN5AkAkAAGDIHwGjy/2sXQ4AABBd/Nb3ze1r8AL6ZAIAAMB2lMsBAAAM0SfTHEEmAACAoYD4xC8+16/BCyiXAwAAwHZkMgEAAAwFrO+bmwIeGfhDkAkAAGDIHwHlcj/lcgAAANRWZDIBAAAMkck0R5AJAABgKGD5dHNTwOXzm2J0OQAAAGxHJhMAAMAQ5XJzBJkAAACG/BKjm5v84g0EmQAAAIasCOiTadEnEwAAALUVmUwAAABD9Mk0R5AJAABgyG/F6OYmv0eWlWQKIwAAANiOTCYAAIChgPgk4HKOLiDeSGUSZAIAABiiT6Y5yuUAAACwHZlMAAAATw38scQLCDIBAADC6pPp7mTsAZfPb4pyOQAAAGxHJhMAAMBQIALWLg8wuhwAACC60CfTHJlMAACAMDKZzJNphj6ZAAAAsB2ZTAAAAEN+y6ebm/wun98UQSYAAIAhfwQM/PF7ZOAP5XIAAIAotm7dOhk4cKCkpqaKz+eTpUuXnvb4t956S66//npp0qSJJCUlSXp6uqxYsSLs8xJkAgAAGApYMRHRwlFaWiqdO3eW7Oxs46BUBZnLli2TvLw86dOnjw5St2zZEtZ5KZcDAABEcbm8f//+upmaMWNGyPbUqVPlL3/5i7zzzjvStWtX4+chyAQAAPCgkpKSkO34+Hjd7BYIBOTw4cPSqFGjsB5HuRwAAMA04Ko0wtytFvj3taSlpUlycnKwZWVlOXIf//jHP8qRI0fkF7/4RViPI5MJAADgqcnYY/THgoICPTCnghNZzAULFsjkyZN1ubxp06ZhPZYgEwAAwIOSkpJCgky7LVy4UO666y5ZvHix9O3bN+zHE2QCAAB4au3yGMfP8frrr8sdd9yhA80BAwac0XMQZAIAABgKiOoT6e6KO4Ewz6/6U+7cuTO4vXv3btm6daseyNOyZUsZP368FBYWyrx584Il8oyMDJk5c6b07NlTioqK9P569erpvp+mGPgDAAAQZibT7RaOzZs366mHKqYfyszM1J9PnDhRb+/Zs0fy8/ODx8+ePVtOnDgho0aNkubNmwfbmDFjwjovmUwAAIAo1rt3b7Gs6ufWfOWVV0K2c3JybDkvQSYAAICnJmOPES8gyAQAADAUUPNUWi73ybTcPb8pb4TCAAAA8BQymQAAAGFMhO52uTrgkRwhQSYAAIChgBWjm5sCLp/flDeuEgAAAJ5CJhMAAMCQX3y6ucnv8vlNEWQCAAAYolxujnI5AAAAbEcmEwAAwJA/AsrVfvEGgkwAAABDlMvNEWQCAAAY8lsxurnJzxRGAAAAqK3IZAIAABiyxCcBl/tkWkxhBAAAEF0ol5tjCiMAAADYjnI5AACAoYDl081NAZfPb4ogEwAAwJBfYnRzk98jhWhvXCUAAAA8hUwmAACAIcrl5ggyAQAADAUkRjc3BTxSiPbGVQIAAMBTyGQCAAAY8ls+3dzkZ3Q5AABAdKFPpjkymQAAAIYsK0YCVozr1+AF3rhKAAAAeAqZTAAAAEN+8enmJr/L5zdFkAkAAGAoYLm/rGPAEk+gXA4AAADbkckEAAAwFIiAgT8BBv7Yb926dTJw4EBJTU0Vn88nS5cuPe3xOTk5+riTW1FRkQNXBwAAol1AfBHRvMBT5fLS0lLp3LmzZGdnh/W4HTt2yJ49e4KtadOmjl0jAAAAPFYu79+/v27hUkFlw4YNHbkmAABQe7DiT5RmMs9Uly5dpHnz5nL99dfL3//+d7cvBwAAeLxPptvNCzyVyQyXCixnzZol3bt3l7KyMnnxxReld+/esnHjRrn88surfIw6TrUKJSUlNXjFAAAA0SGqg8y2bdvqVqFXr16ya9cuefrpp+W1116r8jFZWVkyefLkU/bfnrRbkho495fD8J/NFKflltV39PlbJFwqTjvw3TmOn2Nng8aOn6MgqZHj55Bv4x19+qPNxHGJNTBG73iK3/FzHDsY6/g56u2r6+jzW7HeGGjwY060bOL4OWLHlDt+jhPJCY4+/7EmceK0pE8POPr8MX5n3nd64I3b82SKN96P3si32qhHjx6yc+fOar8+fvx4OXToULAVFBTU6PUBAIDIZUXAyHLLI0FmVGcyq7J161ZdRq9OfHy8bgAAACdTWUzXM5kWQabtjhw5EpKF3L17tw4aGzVqJC1bttRZyMLCQpk3b57++owZM6R169Zy6aWXyvHjx3WfzNWrV8v//u//2n9xAAAA8GYmc/PmzdKnT5/gdmZmpv6YkZEhr7zyip4DMz8/P/j18vJyueeee3TgmZiYKJ06dZK//e1vIc8BAABgKhJGdwcYXW4/NTLcsqpfFV4FmpXdf//9ugEAANiBcrm5WjfwBwAAAM7zVLkcAADATZGwdniA0eUAAADRhXK5OcrlAAAAsB3lcgAAAENkMs0RZAIAABgiyDRHuRwAACCKrVu3TgYOHCipqani8/lk6dKlP/qYnJwcufzyy/UqiG3atDllmkgTBJkAAABhZjLdbuEoLS2Vzp07S3Z2ttHxakXFAQMG6MVr1MqKv//97+Wuu+6SFStWhHVeyuUAAACGrAiYQsgK8/j+/fvrZmrWrFl6We5p06bp7fbt28v69evl6aefln79+hk/D0EmAACAB/tklpSUhOxXpW3VzlZubq707ds3ZJ8KLlVGMxyUywEAADwoLS1NkpOTgy0rK8uW5y0qKpJmzZqF7FPbKqg9duyY8fOQyQQAAPBgJrOgoECSkpKC++3IYtqJIBMAAMCDQWZSUlJIkGmXlJQU2bt3b8g+ta3OVa9ePePnoVwOAACAoPT0dFm1atUPO0Rk5cqVen84yGQCAAB4MJNp6siRI7Jz586QKYrU1ESNGjWSli1byvjx46WwsFDmzZunv/7rX/9ann32Wbn//vvljjvukNWrV8sbb7whf/3rXyUcBJkAAACGLMunm5usMM+/efNmPedlhczMTP0xIyNDT7K+Z88eyc/PD35dTV+kAsqxY8fKzJkz5fzzz5cXX3wxrOmLFIJMAACAKNa7d2+xrOpn16xqNR/1mC1btpzVeQkyAQAADKmJ2N2ejD3g8vlNEWQCAABEcZ9MtzC6HAAAALYjkwkAABDFA3/cQpAJAABgiHK5OYJMAAAAQ2QyzdEnEwAAALYjkwkAABBGJtPt0d0WfTIBAACii5rS/DTzmtfYNXgB5XIAAADYjnI5AABAGKvtqP/cFGDFHwAAgOjC6HJzlMsBAABgO8rlAAAAhtTIch9rlxshyAQAADCkRpa7PrrcEk+gXA4AAADbkckEAAAwxMAfcwSZAAAAhggyzRFkAgAAGGLgjzn6ZAIAAMB2ZDIBAAAMMbrcHEEmAABAWEGmu8tKWkxhBAAAgNqKTCYAAIAhRpebI8gEAAAwpCrVblerLfEGRpcDAADAdmQyAQAADFEuN0eQCQAAYIp6uTGCTAAAAFOWz/UpjMTt8xuiTyYAAABsRyYTAADAECv+mCPIBAAAMMTAH3OUywEAAGA7MpkAAADhDLpxe+CN5Y2BPwSZAAAAhuiTaY5yOQAAAGxHJhMAAMAUk7EbI8gEAAAwxOhyc5TLAQAAYDsymQAAAOGWzPGjCDIBAAAMUS43R5AJAABgioE/xuiTCQAAANuRyQQAADCmVttxe8Udn3gBQSYAAIApyuXGKJcDAADAdmQyAQAATJHJNEYmEwAAwJTli4wWpuzsbGnVqpUkJCRIz549ZdOmTac9fsaMGdK2bVupV6+epKWlydixY+X48eNhnZMgEwAAIIotWrRIMjMzZdKkSfLhhx9K586dpV+/frJv374qj1+wYIGMGzdOH799+3Z56aWX9HNMmDAhrPMSZAIAABiyrMho4Zg+fbqMGDFChg8fLh06dJBZs2ZJYmKizJ07t8rj33//fbnyyivlV7/6lc5+3nDDDTJkyJAfzX6ejCATAAAg3D6ZbjcRKSkpCWllZWWnXG55ebnk5eVJ3759fwj+YmL0dm5ubpX/xF69eunHVASVX375pSxbtkz+8z//M6zvE4JMAAAAD0pLS5Pk5ORgy8rKOuWY4uJi8fv90qxZs5D9aruoqKjK51UZzEcffVSuuuoqqVu3rlx00UXSu3fvsMvljC4HAAAwdYYDb2z17/MXFBRIUlJScHd8fLwtT5+TkyNTp06V5557Tg8S2rlzp4wZM0amTJkiDz/8sPHzEGQCAAAY8lnfNzf5/n1+FWBWDjKr0rhxY4mNjZW9e/eG7FfbKSkpVT5GBZJDhw6Vu+66S29fdtllUlpaKnfffbc8+OCDutxugnI5AACAKbf7Ylo/9Mk0ERcXJ926dZNVq1YF9wUCAb2dnp5e5WOOHj16SiCpAlX9zw9j1BGZTAAAgCiWmZkpGRkZ0r17d+nRo4eeA1NlJtVoc2XYsGHSokWLYJ/OgQMH6hHpXbt2DZbLVXZT7a8INk0QZAIAAHiwT6apwYMHy/79+2XixIl6sE+XLl1k+fLlwcFA+fn5IZnLhx56SHw+n/5YWFgoTZo00QHm448/LuEgyAQAAIjyZSVHjx6tW3UDfSqrU6eOnohdtbNBn0wAAADYjkwmAABAlGcy3UCQCQAAYIogs2aCzE8//VR3FlVLFlV20003nc3TAgAAwOPOKMhUa1j+9Kc/lW3btunRRxVzJqnPFbV8EQAAQNTx4Ohyt5zRwB+1tFDr1q1l3759kpiYKJ988omsW7dOz7908gglAACAaFvxx+0WtZnM3NxcWb16tV6qSM2rpJpaRF1N4vm73/1OtmzZYv+VAgAAwDPOKJOpyuENGjTQn6tA85tvvtGfX3DBBbJjxw5xUnZ2trRq1UoSEhL0LPSbNm067fGLFy+Wdu3a6ePV2pvLli1z9PoAAEAUc3s5yUgYeORkkNmxY0f5xz/+oT9Xgd6TTz4pf//73+XRRx+VCy+8UJyyaNEivTSSmhz0ww8/lM6dO0u/fv102b4q77//vgwZMkTuvPNOnV0dNGiQbh9//LFj1wgAAIAzDDLVMkNqcXVFBZa7d++Wq6++WmcJZ86c6djrqtbRHDFihF5rs0OHDjJr1izdJ3Tu3LlVHq+u5cYbb5T77rtP2rdvL1OmTJHLL79cnn32WceuEQAARC815Mbt/pg+ieI+mSp7WKFNmzby2WefyYEDB+Tcc88NjjC3m5omKS8vT8aPHx/cp/qC9u3bV/cRrYrarzKfJ1/70qVLqz1PWVmZbhVKSkpsuX4AAIDa5IyCzDvuuENnCSv6ZSqNGjWS0tJS+e1vf1ttZvFsFBcX676gFYu5V1DbKsitiloEvqrj1f7qqMFLkydPPmV/vK+OxPtixcsuizvo6PPvSywQp/3zu0aOn2PdVxc5fg7/wTjHz1Gn1NlVY4+ff0KcFlvu/HoRyZ86/74+3Ob7yo+T6pQ6+z3lTxDHNV9R/c9mu5R0bur4OcqSEh0/R5O3P3f0+X2L64nTSqenOvr8J747LuLEy8QURsbO6LfQq6++KseOHTtlv9o3b9488TKVKT106FCwFRQ4HzgBAACPcHvAj+WdgT9hpQlU6VhNvK7a4cOH9YjtCirLqPpkNm3qzF+JahR7bGys7N27N2S/2k5JSanyMWp/OMcr8fHxugEAAKCGMpkNGzbUZXHV7/KSSy7RfTArmgoCVRl91KhR4oS4uDjp1q2brFq1KrhPDT5S2+np6VU+Ru2vfLyycuXKao8HAAA4LbczmFaUZjLXrFmjs5jXXnut/PnPf9YBZ+UgUM2TmZrqXB8LNYgnIyNDryzUo0cPmTFjhu4HqkabK8OGDZMWLVrofpUVKxNdc801Mm3aNBkwYIAsXLhQNm/eLLNnz3bsGgEAQPSKhBV3fNEYZKqATVFTFuXn58sLL7wgu3btkjfffFMHd6+99ppeblKt/uOEwYMHy/79+2XixIl68E6XLl1k+fLlwcE96prUiPMKvXr1kgULFugplyZMmCAXX3yxHlmu5vkEAACAc85o6KbKBg4dOlRuu+02Pcl5xZQ/aqDM1KlTHV1VZ/To0bpVpap102+99VbdAAAAzloklKstid7R5Y899pieCH3OnDlSt27d4P4rr7xSr8QDAAAQldzui2lFeZCp1if/j//4j1P2Jycny8GDzs7FCAAAgCgNMtUUQDt37jxl//r16x1duxwAAMBNri8paXln4M8ZBZlq/XA1cnvjxo16OqNvvvlG5s+fL/fee6+MHDnS/qsEAACIBBUr/rjdonXgz7hx4/Qcldddd50cPXpUl87VBOYqyFTLSgIAAESlSOgTaUn0Bpkqe/nggw/Kfffdp8vmR44ckQ4dOkj9+vXtv0IAAAB4zhkFmZUnYFfBJQAAQG0QCX0ifdGcyQQAAKiVKJc7O/AHAAAAOB0ymQAAAKYioFwubp/fEEEmAACAKcrlxiiXAwAAwHZkMgEAAEyRyTRGkAkAAGCIKYzMUS4HAACA7QgyAQAAYDvK5QAAAKbok2mMIBMAAMAQfTLNUS4HAACA7chkAgAAROGKO24jyAQAADBFn0xjlMsBAABgOzKZAAAAhhj4Y44gEwAAwBTlcmOUywEAAGA7MpkAAACGKJebI8gEAAAwRbncGOVyAAAA2I5MJgAAgCkymcbIZAIAAITZJ9PtFq7s7Gxp1aqVJCQkSM+ePWXTpk2nPf7gwYMyatQoad68ucTHx8sll1wiy5YtC+ucZDIBAACiOJO5aNEiyczMlFmzZukAc8aMGdKvXz/ZsWOHNG3a9JTjy8vL5frrr9dfe/PNN6VFixby9ddfS8OGDcM6L0EmAABAFJs+fbqMGDFChg8frrdVsPnXv/5V5s6dK+PGjTvleLX/wIED8v7770vdunX1PpUFDRflcgAAgHAzmW43ESkpKQlpZWVlVWYl8/LypG/fvj8EfzExejs3N7fKf+Lbb78t6enpulzerFkz6dixo0ydOlX8fn9Y3ycEmQAAAIYiqU9mWlqaJCcnB1tWVtYp11tcXKyDQxUsVqa2i4qKqvw3fvnll7pMrh6n+mE+/PDDMm3aNHnsscfC+j6hXA4AAOBBBQUFkpSUFNxWA3TsEAgEdH/M2bNnS2xsrHTr1k0KCwvlqaeekkmTJhk/D0EmAACABwf+JCUlhQSZVWncuLEOFPfu3RuyX22npKRU+Rg1olz1xVSPq9C+fXud+VTl97i4OKPLpFwOAABgyGtTGMXFxelM5KpVq0IylWpb9busypVXXik7d+7Ux1X4/PPPdfBpGmAqBJkAAABRLDMzU+bMmSOvvvqqbN++XUaOHCmlpaXB0ebDhg2T8ePHB49XX1ejy8eMGaODSzUSXQ38UQOBwkG5HAAAwIPlclODBw+W/fv3y8SJE3XJu0uXLrJ8+fLgYKD8/Hw94ryCGlC0YsUKGTt2rHTq1EnPk6kCzgceeEDCQZAJAAAQxUGmMnr0aN2qkpOTc8o+VUrfsGGDnA3K5QAAALAdmUwAAABDvn83N/nEGwgyAQAAorxc7gaCTAAAAEPhTiHkBLfPb4o+mQAAALAdmUwAAABTlMuNEWQCAACEwyPlardRLgcAAIDtyGQCAAAYYuCPOYJMAAAAU/TJNEa5HAAAALYjkwkAAGCIcrk5gkwAAABTlMuNUS4HAACA7chkAgAAGKJcbo4gEwAAwBTlcmMEmQAAAKYIMo3RJxMAAAC2I5MJAABgiD6Z5ggyAQAATFEuN0a5HAAAALYjkwkAAGDIZ1m6ucnn8vlNEWQCAACYolxujHI5AAAAbEcmEwAAwBCjy80RZAIAAJiiXG6McjkAAABsRyYTAADAEOVycwSZAAAApiiXGyPIBAAAMEQm0xx9MgEAAGA7MpkAAACmKJcbI8gEAAAIs2SOH0e5HAAAALYjkwkAAGDKsr5vbrK8kUolyAQAADDE6HJzlMsBAABgOzKZAAAAphhdbowgEwAAwJAv8H1zk8/l85uiXA4AAADbkckEAAAwRbk8+jKZjz/+uPTq1UsSExOlYcOGRo+5/fbbxefzhbQbb7zR8WsFAADRPbrc7eYFnslklpeXy6233irp6eny0ksvGT9OBZUvv/xycDs+Pt6hKwQAAFGPeTKjL8icPHmy/vjKK6+E9TgVVKakpDh0VQAAAPB0ufxM5eTkSNOmTaVt27YycuRI+fbbb097fFlZmZSUlIQ0AAAAxe0yuY9yeWRQpfJbbrlFWrduLbt27ZIJEyZI//79JTc3V2JjY6t8TFZWVjBrWtm9e7pL3OG6jl1rx3MKxWn5Zec5+vwrC9uK0/518BzHz3HeygTHz3G4pc/xc8SUO/v8zd8Xx1kxfsfPUeeo8+c43iTO8XMca+rs86fd8LWzJxCRE9nO/xxsULTP8XPUfyfZ8XMEhjn7cypGj25xVr17nL3fJ0rLRJY78MQM/PFGJnPcuHGnDMw5uX322Wdn/Py//OUv5aabbpLLLrtMBg0aJO+++6588MEHOrtZnfHjx8uhQ4eCraCg4IzPDwAAUFu5GmTec889sn379tO2Cy+80Lbzqedq3Lix7Ny587R9OJOSkkIaAACA4tVyeXZ2trRq1UoSEhKkZ8+esmnTJqPHLVy4UCf9VLLOUwN/mjRpoltN+ec//6n7ZDZv3rzGzgkAAKKIB0eXL1q0SDIzM2XWrFk6wJwxY4b069dPduzYocetVOerr76Se++9V66++uroHviTn58vW7du1R/9fr/+XLUjR44Ej2nXrp0sWbJEf67233fffbJhwwb9Iq1atUpuvvlmadOmjX5hAQAAaoPp06fLiBEjZPjw4dKhQwcdbKp5x+fOnVvtY1Ssddttt+lxKmdaVfbMFEYTJ06UV199NbjdtWtX/XHNmjXSu3dv/bmKyFU/SkUN7Pnoo4/0Yw4ePCipqalyww03yJQpU5grEwAAnJFIGN3t+/f5T54BR3X5O3k+cDXPeF5enh5zUiEmJkb69u2rB0JX59FHH9VZzjvvvFP+7//+L7qDTDU/5o/NkWlVSh/Xq1dPVqxYUQNXBgAAao0IGl2elpYWsnvSpEnyyCOPhOwrLi7WWclmzZqF7Ffb1Q2uXr9+vV74RlWMz4ZngkwAAAD8QM2AU3mAsh2rGh4+fFiGDh0qc+bM0YOlzwZBJgAAgAfL5UkGs+CoQFF1Idy7d2/IfrVd1YqIal5xNZZl4MCBwX2BQEB/rFOnju6aeNFFF0XXwB8AAADXBazIaIbi4uKkW7duegB05aBRbaenp59yvBpEvW3btuAAa9XUnON9+vTRn59coj8dMpkAAAAe7JNpSk1flJGRId27d5cePXroKYxKS0v1aHNl2LBh0qJFC73qoZpHs2PHjiGPb9iwof548v4fQ5AJAAAQxQYPHiz79+/XM/UUFRVJly5dZPny5cHBQGp6SDXi3G4EmQAAAIZ8lfpEunkN4Ro9erRuVTndctvKj83uUx2CTAAAgChe8cctDPwBAACA7chkAgAAeHAKo0hHkAkAABDFo8vdQrkcAAAAtiOTCQAAYMhnWbq5yeeRgT8EmQAAAKbUCovfr7LonoB4AuVyAAAA2I5MJgAAgCHK5eYIMgEAAEwxutwYQSYAAIApVvwxRp9MAAAA2I5MJgAAgCFW/DFHkAkAAGCKcrkxyuUAAACwHZlMAAAAQ77A981NPo9Mxk6QCQAAYIpyuTHK5QAAALAdmUwAAABTTMZujCATAADAEMtKmqNcDgAAANuRyQQAADDFwB9jBJkAAADh9MkMRMA1eABBJgAAgCH6ZJqjTyYAAABsRyYTAAAgrCmMXK5XW+IJBJkAAACmGPhjjHI5AAAAbEcmEwAAwJQaWe5z+eUKiCcQZAIAABhidLk5yuUAAACwHZlMAAAAUwz8MUaQCQAAYIog0xjlcgAAANiOTCYAAIApMpnGCDIBAABMMYWRMYJMAAAAQ0xhZI4+mQAAALAdmUwAAABT9Mk0RpAJAABgKmCpmrn71+ABlMsBAABgOzKZAAAApiiXGyPIBAAAMGZ9H2i6yhIvoFwOAAAA25HJBAAAMEW53BhBJgAAQFgjuxldboJyOQAAQJTLzs6WVq1aSUJCgvTs2VM2bdpU7bFz5syRq6++Ws4991zd+vbte9rjq0OQCQAAYMoKREYLw6JFiyQzM1MmTZokH374oXTu3Fn69esn+/btq/L4nJwcGTJkiKxZs0Zyc3MlLS1NbrjhBiksLAzntASZAAAAYffJdLuFYfr06TJixAgZPny4dOjQQWbNmiWJiYkyd+7cKo+fP3++/OY3v5EuXbpIu3bt5MUXX5RAICCrVq0K57QEmQAAAGH1yYyEZqi8vFzy8vJ0ybtCTEyM3lZZShNHjx6V7777Tho1aiThYOAPAACAB5WUlIRsx8fH61ZZcXGx+P1+adasWch+tf3ZZ58ZneeBBx6Q1NTUkEDVBH0yAQAATEVQuTwtLU2Sk5ODLSsry/b7+MQTT8jChQtlyZIletBQOMhkAgAAmNIzGLk8hZH1/YeCggJJSkoK7j45i6k0btxYYmNjZe/evSH71XZKSsppT/PHP/5RB5l/+9vfpFOnTmFfJplMAAAAD0pKSgppVQWZcXFx0q1bt5BBOxWDeNLT06t97ieffFKmTJkiy5cvl+7du5/R9ZHJBAAAiOIVfzIzMyUjI0MHiz169JAZM2ZIaWmpHm2uDBs2TFq0aBEst//hD3+QiRMnyoIFC/TcmkVFRXp//fr1dTNFkAkAAGAqoOaoDG+eSmeuwdzgwYNl//79OnBUAaOamkhlKCsGA+Xn5+sR5xWef/55PSr95z//ecjzqHk2H3nkEePzEmQCAABEudGjR+tW3eTrlX311Ve2nJMgEwAAIIrL5W4hyAQAADBFkGmM0eUAAACwHZlMAAAAU3pJR5fL1QHK5QAAAFHFsgK6uX0NXkAmEwAAIJw+mW5nEi1vZDLpkwkAAADbkckEAAAIK4tIJtMEQSYAAEA4q+34XO4TaXmjTyblcgAAANiOTCYAAIApyuXGCDIBAAAMWYGAWC6Xyy3K5QAAAKityGQCAACYolxujCATAADAlJqI3ccURiYYXQ4AAADbkckEAAAIq1zu9jyZlngBQSYAAIAhK2CJ5XK53CLIBAAAiDJ6+iC3M5kB8QLP9Mm86aabpGXLlpKQkCDNmzeXoUOHyjfffHPaxxw/flxGjRol5513ntSvX19+9rOfyd69e2vsmgEAAGorzwSZffr0kTfeeEN27Nghf/7zn2XXrl3y85///LSPGTt2rLzzzjuyePFiWbt2rQ5Kb7nllhq7ZgAAEIXl8ghoXuCZPpkqYKxwwQUXyLhx42TQoEHy3XffSd26dU85/tChQ/LSSy/JggUL5Nprr9X7Xn75ZWnfvr1s2LBBrrjiihq9fgAAEAUol0dfkFnZgQMHZP78+dKrV68qA0wlLy9PB6B9+/YN7mvXrp0uuefm5lYbZJaVlelWOVhVyku/Eycds06I08rKnP03+I/+8Lo5JXA01vFz+MsdP4X4y3yOn8Ny+N9x4jvnv2etGOdfJznhd/wU/uM10H/K4bffidIaeH9bzv6MUiyn3xjqHDXwWuHHnTha7sggmRPynYjLicQT6hq8wPKQ+++/30pMTFS31rriiius4uLiao+dP3++FRcXd8r+n/zkJ/p5qjNp0iT9/DReA74H+B7ge4DvAb4HvP89sGvXLltikGPHjlkpKSmu/3vk301di7qmSOZT/3MrwFUl7z/84Q+nPWb79u06A6kUFxfrLObXX38tkydPluTkZHn33XfF5zs126HK5MOHDw/JSio9evTQ/TurO+/JmcyDBw/q8nx+fr4+H9xTUlIiaWlpUlBQIElJSdwK7gV4b0Qcfk5FDlWJVNXLf/3rX9KwYUNbnlMNKC4vr4Fyl4G4uDg9GDqSuVouv+eee+T2228/7TEXXnhh8PPGjRvrdskll+i+lSrgUP0r09PTT3lcSkqK/kZQQWLlby41ulx9rTrx8fG6nUwFmAQ2kUHdB+5FZOBeRBbuR+TgXkSOmBj7xjiroC7SA7tI4mqQ2aRJE93ORCDwfR+nkzOVFbp166b7a65atUpPXaSokekqI1lVUAoAAIBaNvBn48aN8sEHH8hVV10l5557rp6+6OGHH5aLLrooGDAWFhbKddddJ/PmzdMlcZV5vPPOOyUzM1MaNWqk/6r87W9/q49nZDkAAICzPBFkJiYmyltvvSWTJk2S0tJSPRn7jTfeKA899FCwtK1GkqtM5dGjR4OPe/rpp3WaXGUyVcazX79+8txzz4V1bvX86rxVldBRs7gXkYN7EVm4H5GDexE5uBfuc3XgDwAAAKKTZ1b8AQAAgHcQZAIAAIAgEwAAAJGPTCYAAABsR5B5GtnZ2dKqVSs98WrPnj1l06ZN9t8BhMjKypKf/OQn0qBBA2natKkMGjRIzxpw8ooLo0aNkvPOO0/q16+vZw9Qk+zDWU888YReXev3v/8998Ilaqq2//qv/9Lf+/Xq1ZPLLrtMNm/eHPy6Gsc5ceJEPQOH+nrfvn3liy++cOtyo5bf79fT6LVu3Vq/zmo6vSlTpoSskc29cMa6detk4MCBkpqaqn8eLV26NOTrJq+7Wjnwtttu01MbqsVa1HSHR44cceiKazeCzGosWrRIz7Gppi/68MMPpXPnznoKpH379tXsHapl1q5dqwNItZLTypUr9dRUN9xwg566qsLYsWPlnXfekcWLF+vjv/nmG7nllltcve5op+apfeGFF6RTp04h+7kXNUctjXfllVfqRSbee+89+fTTT2XatGl67uAKTz75pPzpT3+SWbNm6fmFzznnHP1zS/1hBvuoZYmff/55efbZZ/XSx2pbvfbPPPMM98Jh6neB+n2skkBVMXkPqADzk08+0b9j1NLUKnC9++67nb702sntxdMjVY8ePaxRo0YFt/1+v5WammplZWW5el21zb59+1RqwFq7dq3ePnjwoFW3bl1r8eLFwWO2b9+uj8nNzXXxSqPX4cOHrYsvvthauXKldc0111hjxozR+7kXNeuBBx6wrrrqqmq/HggErJSUFOupp54K7lP3KD4+3nr99ddr6CprhwEDBlh33HFHyL5bbrnFuu222/Tn3IuaoX7uL1myJLht8rp/+umn+nEffPBB8Jj33nvP8vl8VmFhYQ1dee1BJrMKas3zvLw8nWavoCZ1V9u5ubk1+TdArXfo0CH9GqhVmxR1X1R2s/K9adeunbRs2ZJ74xCVWR4wYEDIa869qHlvv/22dO/eXW699VbdlaRr164yZ86c4Nd3794tRUVFIfdJrXymuvrwc8tevXr10ksWf/7553r7H//4h6xfv1769+/PvXCRyXtAfVQlcvVeqqCOV7/jVeYTtXDFn5pWXFys+9w0a9YsZL/a/uyzz1y7rtpGrU+v+v+pEmHHjh31PvUDJC4uTv+QOPneqK/BXgsXLtTdRVS5/GTci5r15Zdf6hKt6sYzYcIEfU9+97vf6fdDRkZG8Pu/qp9bvDfsNW7cOCkpKdF/4MbGxurfF48//rguwyrcC3eYvO7qo/ojrbI6deroRAbvE/sRZCKiM2gff/yxzhCg5hUUFMiYMWN0vyU1+A3u/9Glsi9Tp07V2yqTqd4fqu+ZCjJRc9544w2ZP3++LFiwQC699FLZunWr/oNYDUbhXgA/oFxehcaNG+u/Tk8esay2U1JSqnoIbDZ69GjdIXvNmjVy/vnnB/er1191Zzh48CD3xmGqa4Ia6Hb55Zfrv/RVUwOtVKd69bnKDnAvao4aLduhQ4eQfe3bt5f8/Hz9ecXPJn5uOe++++7T2cxf/vKXeoT/0KFD9SA4NTsG98I9Ju8B9fHkAbwnTpzQI875/W4/gswqqPJTt27ddJ+bylkEtZ2enu7AbUAF1ZdbBZhLliyR1atX6ylCKlP3RY2urXxv1BRH6hct98Ze1113nWzbtk1naSqayqSpkmDF59yLmqO6jZw8nZfqE3jBBRfoz9V7Rf2SrPzeUCVd1c+M94a9jh49qvvwVaYSE+r3BPfCPSbvAfVRJSnUH9EV1O8ade9U303YzO2RR5Fq4cKFekTaK6+8okej3X333VbDhg2toqIity8tqo0cOdJKTk62cnJyrD179gTb0aNHg8f8+te/tlq2bGmtXr3a2rx5s5Wenq4bnFd5dDn3omZt2rTJqlOnjvX4449bX3zxhTV//nwrMTHR+p//+Z/gMU888YT+OfWXv/zF+uijj6ybb77Zat26tXXs2LEavtrolpGRYbVo0cJ69913rd27d1tvvfWW1bhxY+v+++8PHsO9cG62iy1btuimQpjp06frz7/++mvj1/3GG2+0unbtam3cuNFav369nj1jyJAhDl1x7UaQeRrPPPOMDmbi4uL0lEYbNmyouTtTS6kfGlW1l19+OXiM+mHxm9/8xjr33HP1L9mf/vSnOhBFzQeZ3Iua9c4771gdO3bUfwC3a9fOmj17dsjX1RQuDz/8sNWsWTN9zHXXXWft2LGjhq8y+pWUlOj3gfr9kJCQYF144YXWgw8+aJWVlQWP4V44Y82aNVX+jlCBv+nr/u233+qgsn79+lZSUpI1fPhwHbzCfj71P7uzowAAAKjd6JMJAAAA2xFkAgAAwHYEmQAAALAdQSYAAABsR5AJAAAA2xFkAgAAwHYEmQAAALAdQSYAAABsR5AJACLSqlUrmTFjBq8FANiEIBMAAAC2Y1lJALVC7969pWPHjvrz1157TerWrSsjR46URx99VPr06SNr164NOZ4VdwHg7JDJBFBrvPrqq1KnTh3ZtGmTzJw5U6ZPny4vvviivPXWW3L++efrgHPPnj26AQDOTp2zfDwAeEZaWpo8/fTT4vP5pG3btrJt2za9PWLECImNjZUGDRpISkqK25cJAFGBTCaAWuOKK67QAWaF9PR0+eKLL8Tv97t6XQAQjQgyAQAAYDuCTAC1xsaNG0O2N2zYIBdffLEulcfFxZHRBAAbEWQCqDXy8/MlMzNTduzYIa+//ro888wzMmbMmOA8mevWrZPCwkIpLi52+1IBwPOYwghArZnC6NJLL5VAICALFizQ2Us1hdFjjz2m+2mqrOZ///d/6wC0rKyMKYwA4CwRZAKoNUFmly5dWNUHAGoI5XIAAADYjiATAAAAtqNcDgAAANuRyQQAAIDtCDIBAABAkAkAAIDIRyYTAAAAtiPIBAAAgO0IMgEAAGA7gkwAAADYjiATAAAAtiPIBAAAgNjt/wGkbAT7WgPlNQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "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": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGwCAYAAAC3qV8qAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAARX5JREFUeJzt3Ql4U1XawPG3K21KKVBk07KIArKJglZGRRwYKiAjijoKCiojo4IKCIO4YAEFRAQZx9Fxxc8BGf0+FkUGWWRR2WGQIoqgII5sDkuhpHTN97ynJCTYpkmbtL3J//c8l9zkniS3t6V5+573nBPhcDgcAgAAYCGRlX0CAAAA/iKAAQAAlkMAAwAALIcABgAAWA4BDAAAsBwCGAAAYDkEMAAAwHKiJUQVFhbK/v37JTExUSIiIir7dAAAgA90erqTJ09Kw4YNJTIyMvwCGA1eUlJSKvs0AABAGfz0009ywQUXhF8Ao5kX5wWoUaNGZZ8OAADwwYkTJ0wCwvk5HnYBjLPbSIMXAhgAAKyltPIPingBAIDlEMAAAADLIYABAACWE7I1MACA8FRQUCB5eXmVfRooQUxMjERFRUl5EcAAAEJm/pCDBw/K8ePHK/tUUIqaNWtK/fr1yzVPGwEMACAkOIOXunXris1mYxLTKhpk2u12OXz4sLnfoEGDMr8WAQwAICS6jZzBS3JycmWfDryIj483txrE6PerrN1JFPECACzPWfOimRdUfc7vU3lqlQhgAAAhg7Xvwuf7RAADAAAshwAGAABYDgEMAAAB1qVLFxk2bBjXNYgIYAAAqEQrV640NSEVPX9Nenq6tG/fXqyKAAYAAFgOAQwAAOVw6tQpGTBggFSvXt1MzPbiiy96HH/vvfekY8eOkpiYaGaf7devn2sit71798r1119v9mvVqmUyMffcc4+5v3jxYrnmmmvMrLU6t82NN94o33//vet1c3NzZejQoeY94+LipHHjxjJp0iTXcc3o/PGPf5TzzjtPatSoIb/97W/lq6++Msdmzpwp48aNM/f1PXXTx6yEAAYe7Hl2aftuW7PpPgDAu1GjRsmqVatkwYIFsmTJEtMltGXLFtdxnetkwoQJJliYP3++CVqcQUpKSor83//9n9nfuXOnHDhwQGbMmOEKjEaMGCGbNm2S5cuXS2RkpNx8881SWFhojv/lL3+Rjz76SD744APz3FmzZkmTJk1c73vbbbeZQOlf//qXbN68WS6//HLp2rWrHD16VP7whz/IY489Jq1btzbvqZs+ZiXMxAsP2XkFHvu2GC4QAJQkKytL3nrrLfnHP/5hggP17rvvygUXXOBqc99997n2L7zwQhN4XHHFFea5mrWpXbu2Oaaz0mq2xalv374e7/X222+bbMqOHTukTZs2sm/fPrn44otNlkYzKJqBcfriiy9kw4YNJoCpVq2aeWzq1KkmgPrf//1fGTx4sHnv6OhokxWyIjIwAACUkXbpaFdOamqq6zENSFq0aOG6r9mP3r17S6NGjUw30nXXXWce1wDEm127dsmdd95pgh7tAnJmV5zPu+eee2Tr1q3mvR555BGT/XHSbI8GSNr1pIGKc9uzZ49HN5SVkYEBACBItBsoLS3NbNrFoxkUDUD0vgY+3mjQo1mVN954Qxo2bGi6jjTz4nze5ZdfbgIS7SJatmyZ3H777dKtWzeTYdHgRWtjtDvrXO5ZHisjgAEAoIyaNWsmMTExsn79epNhUceOHZPvvvvOZFq+/fZbOXLkiEyePNnUuyitaXEXGxvrWpDSSZ+jdS0avFx77bWubqFz1ahRw9Su6HbrrbfKDTfcYGpcNLjR1bm1i8i9Lubc93V/T6shgAEAoIy0W2bQoEGmkFe7a7SO5cknnzQFt0qDGg0UXn75ZXnggQdk+/btpqDXnWZZtIZl4cKF0rNnT7Nas45I0td7/fXXTSZFszaPP/64x/OmTZtmjl122WXm/T788ENTz6IZFs3EdOrUSfr06SNTpkyR5s2by/79++WTTz4xhcA6KkoDG83gaDeU1uxo95azXsYKqIEBAKAcXnjhBZMl0S4fDRy0qLZDhw7mmHYZ6fBkDS5atWplMjFaTOvu/PPPN0OaNUCpV6+eGRqtAcmcOXNM/Yx2Gw0fPty8j7vExEQTnGgwokXBOrpp0aJF5rkaEOl+586d5d577zUBzB133CE//vijeQ9nkbBmbHQYt57n+++/b6mfgwiHw+GQEHTixAlJSkqSzMxMk2ILV/rtzc7P9rn9EXuW9JxfVEm/qM9ySbZV9+l58dHxrAILoNKcPn3aZBOaNm1q5kSBdb9fvn5+04UU4nQul6vev8rn9o7CaIk4k5frMTdNIiLzfXreujvXSUJsQllPEwAAv9CFFOJO5xeG1PsAAKDIwISRrO+eEkdhUbV7yXIlseWzZ9prwVjJ7SMic6V686K2AABUJAKYMLJ6VFqpNS1H7Cel54KioOTz0d0l2ZbopW2Wqy0AABWJACaM2GKjxBbr/VuenX/2uLb11j47Pyqg5wcAgK+ogQEAAJZDAGNBrBgNAAh3BDAAALix5+ZLk8c/MZvuo2oigAEAAJZDAAMAQBj4+uuvzfIBugaSLjXw0ksvSVgFMKtXrzbrPejS3noB5s+f73FcHytuc1/DwXnx3DddH8Ldtm3bzNoSOsWwruCp6z2gSHZeQbH7gaDDpjMGZpjN2xBqAIC12O12ufDCC83nrS76aHV+BzCnTp2SSy+9VF555ZVijx84cMBje/vtt02AolGfu/Hjx3u0e/jhhz3WQejevbtZoVMXstLgJz093azKCQBAqDl58qT0799fEhISzArT06dPly5dusiwYcPM8ZycHBk5cqRZ+FHbpKamysqVK13P1wUjdRXqTz/9VC655BKzSrYu1Kifr0664KN+nuqijlZadTpg88D06NHDbCU5N6pbsGCBWelSo75zV9EsKQKcNWuW5ObmmuBHlyFv3bq1We5blw4fPHhwsc/Rb65u7kEQAEDCezHbMmSp3Qt3y1rEGx8T5dcCtyNGjJAvv/xSPvroI7Na9NixY2XLli3Svn17c1xXqN6xY4dZoVp7QObNm2cClIyMDLn44ouLztVuNytdv/fee2ZF6rvuussEPfqZGoqCOpHdoUOH5JNPPpF33333V8c0hTVhwgRp1KiR9OvXzywVHh1ddDpr1641S4Br8OKUlpYmzz//vBw7dkxq1ar1q9ebNGmSWY48HFaNtudle+zH55U8oZw/K1EDQCjR4KXV2E/L9Rodn11epuftGJ9W6sSh7tkX/ZycPXu2dO3a1Tz2zjvvmEBF7du3z9zXW+djI0eOlMWLF5vHJ06caB7Ly8uT1157TZo1a+YKerS3I1QFNYDRb4hmWm655RaPxx955BG5/PLLpXbt2rJmzRoZM2aMSXNphkUdPHjQLLHtTiNS57HiAhh9DY1g3TMwWjtjFRpopM5ODeqK0QCAqueHH34wwceVV17peiwpKUlatGhh9jXLUlBQIM2bN/d4Xk5OjiQnJ7vu22w2V/CitCvq8OHDEqqCGsBoF5D26Wkhrjv3QKNdu3Ym0/KnP/3JZFHK2i+nzwuFPj0AQGBoN45mQvyl3UbOzMump7r6nEk5970DJSsrS6KiokxNqN66q1797Pp2MTExHse0C0uz+6EqaAHM559/Ljt37pR//vOfpbbVYqT8/HzZu3eviTi1Nka7n9w574dC5XRx3H/I/nXzcomPji+x7TF7ltz8SXezP7/3IqnlZYFGe26BXPv8CrMfF+UZSAJAKNMP8LIEH+5KWxMuELRGVIOPjRs3mrIKlZmZKd99950pp7jssstMBkazKTo6F0WC9l156623pEOHDmbEUmm0QFcLjurWrWvud+rUSZ588kmTUnNGlEuXLjXBTXHdR6HgdH6ha/+aSV+KOM7W//xKRK4ktiza7fbiOu9tjaLj/hSUAQAqhpZaDBw4UEaNGmVKK/Sz8JlnnjGfi/p7W7uOtDdjwIAB8uKLL5qA5pdffpHly5ebXoxevXr59D46OEYLgZ37P//8s/n81SzORRddJCE/jFpTWfoF66b27Nlj9rW4yL3+5MMPP5Q//vGPv3q+Fujq5DlfffWV6ffT6mgt4NVqaWdwokW92q00aNAgM/GOZnFmzJjh0fUEAECo0BpQ/eP9xhtvlG7dusnVV19thkM7SzC0WFcDmMcee8z8Md+nTx+PjI0v9u/fb4If3bTuVEcs6X5xn9WW4PDTihUrtK/jV9vAgQNdbf7+97874uPjHcePH//V8zdv3uxITU11JCUlOeLi4hyXXHKJY+LEiY7Tp097tPvqq68c11xzjaNatWqO888/3zF58mS/zjMzM9Ocl95awX9PnXC0mdnGbPuOHXOcyskrcdt37Khb26Ne27pvhYWFQTtn3QeAypKdne3YsWOHuS0v/X3ZePRCs+l+ZcjKyjKfk2+++aYj3L5fmT5+fvvdhaQT65RWFKRztZQ0X4uOPlq3bl2p76NpMa2jCUe22Civfa7Z+WePaTtvbQEAVd+///1v+fbbb81IJK1/cQ5/vummmyr71KosPvkAAHCjfxTunexbXUkgaZeODn7REgqtIdU/4uvUqcP3pgQEMAAAVDKtRdFh0vAdq1EDAADLIQNjQTpHzMlvilbv9jZfDAAAoYoMDAAAsBwyMBZUWQVmAABUFWRgAACA5RDAAADgLveUSHpS0ab7qJIIYAAAgOUQwAAAEAbeeOMNs5q1rjuom665tGHDBrEqAhgAAMLAypUr5c4775QVK1aYhZVTUlKke/fuZlVqKyKAAQCgkp08eVL69+8vCQkJ0qBBA5k+fbpZe3DYsGHmeE5OjowcOVLOP/980yY1NdUEJE4zZ86UmjVryqeffmpWsa5evbrccMMNZtVpp1mzZslDDz0k7du3l5YtW8qbb74phYWFsnz5crEiAhgAQGjShYe1CNfvzX72NXS/LK9RyqLH5xoxYoR8+eWX8tFHH8nSpUvNOkhbtmxxHR86dKjJmsyZM0e2bdsmt912mwlQdu3a5Wpjt9vNekrvvfeerF69Wvbt22eCnpJo+7y8PKldu7ZYEfPAAABCU55dZGLD8r3G1IvK9rwn9ovEJvicfXn33Xdl9uzZ0rVrV/PYO++8Iw0bFp27BiJ6X2+dj40cOVIWL15sHp84caJ5TIOR1157TZo1a+YKepyrWhdn9OjR5vW0FsaKCGAAAKhEP/zwgwk+rrzyStdjSUlJ0qJFC7OfkZEhBQUF0rx5c4/n5eTkSHJysuu+zWZzBS9Ku6IOHz5c7HtOnjzZZHO0GyouLk6siAAGABCaYmxFmRB/abeRM/MycrdIrK1s7x0gWVlZEhUVZVar1lt31atXP/uWMTEexyIiIsRRTFeWdjNpALNs2TJp166dWBUBDCqEPTdfWo391OzvGJ9mlkMAgKCKiPC5G6dEGryU9zVKceGFF5rgY+PGjdKoUSPzWGZmpnz33XfSuXNnueyyy0wGRrMpOgy6PKZMmSLPPfecKfbt2LGjWBmfIgAAVKLExEQZOHCgjBo1yhTU1q1bV5555hmJjIw0WRTtOtIRSgMGDJAXX3zRBDS//PKLGT2kGZRevXxbG+/555+XsWPHmlqbJk2ayMGDB11ZHPdMjlUwCgkAgEo2bdo06dSpk9x4442mqPbqq682w6Gd9SlarKsBzGOPPWZqY/r06eORsfHFq6++Krm5uXLrrbea+hjnpl1KVkQGBgCAKpCF0XlanE6dOiXjxo2TwYMHm/vaxaT3dSvOPffcYzZ3GuS418Ds3btXQgkBDAAAlezf//63fPvtt2Ykkta/OIc/33TTTZV9alUWAQwAAO60aDc9s8KviXbl7Ny5U2JjY6VDhw5mMrs6derwvSkBAQwCwp5bIPHR+V6O5xe774v4mChTyAYAoUoLc3WYNHxHAIOAuPb5FSKOWJ/adnzWv3U3GHYNADgXo5AAAIDlkIFBmcVFn41/Nz/dTeKj40tsq91GzszLpqe6ljqRnXZJdXx2Gd8dAECxCGBQZu51KVqnYovx7cdJgxdm4gUAlAddSAAAwHIIYAAAcGPPs0vbd9uaTfdRNRHAAAAAyyGAAQAgDLzxxhtmNetatWqZTddc2rBhQ8Dfp0uXLjJs2DAJNgIYVAgt2t07uZfZKOAFgIq3cuVKufPOO2XFihWydu1aSUlJke7du8vPP/9c4nN08ceqigAGAIBKdvLkSenfv78kJCSYFaKnT5/ukcnIycmRkSNHyvnnn2/apKammoDEaebMmVKzZk359NNPzSrW1atXlxtuuEEOHDjgaqOLRT700EPSvn17admypbz55ptSWFgoy5efnVy0SZMmMmHCBLPydY0aNcxikrp69dChQ11t9Jx0FKqu3eQMcvScli1bZhaUXLVqlcyYMcO00S1Yi0gSwAAAQpKuxKxFuP5u2fnZrtfQ/bK8hvsq0L4YMWKEfPnll/LRRx/J0qVLzTpIW7ZscR3XAEKzJnPmzJFt27bJbbfdZgKUXbt2udrY7XazntJ7770nq1evln379pmgpyTaPi8vT2rXru3xuL7GpZdeahaYfPrpp+W6667zCJY0QNE1mpyPbdy40bzOb37zGxO4dOrUSe6//34TPOmmmZ5gYB4YVAj9D506O9Xsr++3XmwxNq48gKDS4MP5e6esunzQpUzP8+f3nGZf3n33XZk9e7Z07drVPPbOO+9Iw4YNzb4GInpfb52PjRw5UhYvXmwenzhxonlMg4jXXntNmjVr5gp6nKtaF2f06NHm9bQWxt1vf/tbeeyxx1z3NRP06KOPyi+//CLR0dGyY8cOE9hoAPPAAw+Y2yuuuEJstqKvVxej1P369etLMBHAAABQiX744QcTfFx55ZWux5KSkqRFixZmPyMjQwoKCqR58+Yez8vJyZHk5GTXfQ0anMGL0q6ow4cPF/uekydPNtkcDT7i4uI8jnXs2NHjfps2bUyWRjMvGpzowpM33nijvPLKK+a4Pq5BTkXzO4DRtNQLL7xgVs3U1NC8efOkT58+ruPa/6WRpLu0tDQTKTodPXpUHn74Yfn4448lMjJS+vbta9JO2mfnpCmyIUOGmNTUeeedZ9r/+c9/LvtXCgAIK7q8iWZCypK5cWZeVt6+0usyKd7eO1CysrIkKirKfO7qrbvqbp+bMTExHse0/qS4riztItIARmtW2rVr96vjWs9y7ut07tzZBDvVqlUzwYo+TwOo7du3y5o1a7x2VVWZAObUqVOmb+y+++6TW265pdg22i+naS0n/YLdaaGSBj/az6dR57333msKhTR9pk6cOGEqozWtpekwjT71/bRASdsBAFAa/eAtb3e1BiLB7vK+8MILTfChf7A3atTIPJaZmSnfffedCRw046EZGM2m6DDo8pgyZYo899xzptj33EyLN1oHo8Ow9fNcn6/JBz03TWhoIHP11Ve72mqWRs+3ygUwPXr0MJs3+gWW1Pf1zTffmGyMfqOcF+/ll1+Wnj17mqhQ++O0Ulqrmt9++21zIVq3bi1bt26VadOmEcAAAEJKYmKiDBw4UEaNGmW6aurWrSvPPPOMCRI0CNOuI/3DX0cGvfjiiyag+eWXX8zoIc2E9OrVy6f3ef7552Xs2LEmWaCjjQ4ePOjK4rhncoqjWZfhw4ebz+RrrrnG9ZhmXrT+xT1ro6+9fv16M/pIX1e/Jv1aLDEKSdNM+g3Q/rsHH3xQjhw54jqmVdSaSXGP/DTTol+cfsHONhrZ6YVy74bauXOnHDt2rNj31AhQMzfuGwAAVqB/oOvoHa0t0c9EzWjocGhnfYr2amgAo8W1+tnap08fj4yNL1599VWTHNBh0Vof49w0eVCatm3bms9uHYLtDHY0gNFMy7n1LxrUaFdXq1atTAmIFh9boohXu4+0a6lp06by/fffyxNPPGEyNhqU6BekEZ8GNx4nER1tIjRnNKi3+nx39erVcx3TGQTPNWnSJBk3blygvxz4yH3YYWnHS2tb1KZAJKJoAiV/hyMCgBWzMNr74F6uoZ9pzrIJ7WLS+yV9zt1zzz1mc6dBjvvvT1/mYympjSYZtH7VnQYzxf1+1oyRfuYHW8ADmDvuuMMjYtP0llZFa1bGOTwsGMaMGWPG0TtpBiZYY89RvqGGvrZNbFl0e7qgmySIZ3EaAIQSnXNFJ4bTkUha/+Ic/nzTTTdV9qlVWUEfRq3FSTrhze7du00Ao7Ux5w7rys/PN5Gds25Gbw8dOuTRxnm/pNoarbs5t1gYAAB/adFuxsCMCr9w2pWjpRJaPtGhQwczmZ1+fqKSApj//Oc/pgZG+9mU9vEdP37cDAfTb5D67LPPzHTGOjWys82TTz5pRig5h4XpiCXt9yuu+wgBlHtKZGLRREnyxH6RWM/hdGUdoujvsMSj2VnSY17wMnYAUJVoYa5+LiKIAYyOR9dsitOePXvMCCGtYdFN++d0XhfNlGgNjM7dctFFF5kiXKVFSVono9MM6xBpDVJ0tkDtenLOMNivXz/zOoMGDTIzBeo4c50nRteGgPWHKPoyLDE7L/hD8AAA1uX3KKRNmzaZSFE3pXUnuq9Ds7RIVyeg+/3vf2+KeDQAcabB3Lt3tFBJF5LSLiUdPq1Dsl5//XWPGQiXLFligiN9vlZd6+szBwwAwBuK/sPn++R3BkaHS3l7Y50cpzSaqXFOWlcSLf7VwAcAgNI4yw10gcL4+MDNgovg0O9TcbMH+4O1kAAAlqc9ADpPiXOQiK4LpN3cqIIrhNvt5vuk369zl0bwBwEMACAkOEeplrSAIaoODV7Ku1o1AUyo0+6+vKJUnU9y7cXvl0aLcr39teN+Drof5LVFAIQfzbjoiFedLFUHiKBq0m6j8mRenAhgQp0GC85h0f6aepHvbUsZcm2LjpeMPWemkw7gKq0AcC79cAzEBySqtqCshQQAABBMZGDCycjdIrGldN1ot5Ez81Jae/e2AABUIAKYcKLBiJdunnK3BwCgghDAIDBKK/j1tzhYlzRwYjVqAMA5CGAQGP50JfnQNl5HNDVJcRvBVKMcJwcACDUU8QIAAMshA4Oy07lcdPi0L/wpDtbFHDMPiSy6ie8OAKBYBDDwpEW76Zm+XRXt5ilLka8vxcFMdAcA8IIuJAAAYDkEMAAAwHIIYAAAgOUQwFiRzpGSnlS0uc+XAgBAmKCIF1WvOBgAgFKQgQEAAJZDAAMAACyHAAYAAFgONTBVhfuChVqYGx0ZmIURfVk4EQAAiyGAqSrMgoVF4me09H0FZn8WUQQAIETQhQQAACyHDEwVlP3QFrEl1QvYwogurC8EAAgRBDBVkQYavi6S6MvCiAAAhBi6kAAAgOUQwAAAAMshgAEAAJZDAAMAACyHIl4rYmFEAECYIwMDAAAshwAGAABYDgEMAACwHAIYAABgOQQwAADAcghgAABA6Acwq1evlt69e0vDhg0lIiJC5s+f7zqWl5cno0ePlrZt20pCQoJpM2DAANm/f7/HazRp0sQ8132bPHmyR5tt27bJtddeK3FxcZKSkiJTpkwpz9eJEGbPs0vbd9uaTfcBAKHP7wDm1KlTcumll8orr7zyq2N2u122bNkiTz/9tLmdO3eu7Ny5U37/+9//qu348ePlwIEDru3hhx92HTtx4oR0795dGjduLJs3b5YXXnhB0tPT5fXXXy/L1wgAAMJ9IrsePXqYrThJSUmydOlSj8f++te/ypVXXin79u2TRo0auR5PTEyU+vXrF/s6s2bNktzcXHn77bclNjZWWrduLVu3bpVp06bJ4MGD/T1lAAAQYoJeA5OZmWm6iGrWrOnxuHYZJScny2WXXWYyLPn5+a5ja9eulc6dO5vgxSktLc1kc44dO1bs++Tk5JjMjfsGAABCU1CXEjh9+rSpibnzzjulRo0arscfeeQRufzyy6V27dqyZs0aGTNmjOlG0gyLOnjwoDRt2tTjterVq+c6VqtWrV+916RJk2TcuHHB/HJQSbLzT3utbcnOzy52vzTx0fEmuAYAWE/QAhgt6L399tvF4XDIq6++6nFsxIgRrv127dqZTMuf/vQnE4RUq1atTO+nQZD762oGRot/YX09Puntc9suH3Txue36fuvFFmMr41kBAEIugHEGLz/++KN89tlnHtmX4qSmppoupL1790qLFi1MbcyhQ4c82jjvl1Q3o4FPWYMfAAAQ5gGMM3jZtWuXrFixwtS5lEYLdCMjI6Vu3brmfqdOneTJJ580rxUTE2Me0+JgDW6K6z5C6ImLipP1e38y+9mPfivxCSUHwdpt5My8rLx9peka8qUtACCMApisrCzZvXu36/6ePXtMAKL1LA0aNJBbb73VDKFeuHChFBQUmJoVpce1q0gLdNevXy/XX3+9GYmk94cPHy533XWXKzjp16+fqWcZNGiQqaHZvn27zJgxQ6ZPnx7Irx1VmNam2ByOojvR8T539WjwQrcQAIQ+vwOYTZs2meDDyVl3MnDgQDNXy0cffWTut2/f3uN5mo3p0qWL6eaZM2eOaasjh7RYVwMY9/oVHY69ZMkSGTJkiHTo0EHq1KkjY8eOZQg1AAAoWwCjQYgW5pbE2zGlo4/WrVtX6vtoce/nn3/u7+khDGXnFXjs24p6HQEAIYy1kAAAgOUQwAAAAMsJ6kR2QEXQwt2T3xQtBuptBBIAIHQQwKDKs+cWiOTmezmeX+x+cbLzC3yu1wIAVF0EMKjyrp2yQrIlzqe2HZ9d7r1BRK4ktizaPZ1fKAlnl9sCAFgINTAAAMByyMCgSoqPiXLtb36qm0hsQolttdvImXnZ9FRXscWW/GN9xJ4lPRcE+GQBABWOAAZVkvsq0SYg8RKUuNO23gKY7PyzgREAwLroQgIAAJZDAAMAACyHAAYAAFgOAQwAALAcinhheVq0u3dyr8o+DQBABSIDAwAALIcABgAAWA4BDKwv95RIelLRpvsAgJBHAAMAACyHAAYAAFgOAQwAALAchlGj6su1+3681LZuNTIORzlPDABQWQhgUPVNvShgbeN1kcgmKUV38jTYqVHOkwMAVAa6kAAAgOWQgUHVFGMTeWK/b22128iZeRm5WyTWVmLT7MxDIotuCtBJAgAqCwEMqibt6olN8P95Grx4e54GRgAAy6MLCQAAWA4BDAAAsBwCGAAAYDkEMAAAwHIo4oX1adFuemZlnwUAoAKRgQEAAJZDAAMAACyHAAYAAFgOAQwAALAcAhgAAGA5BDAAAMByCGAAAEDoBzCrV6+W3r17S8OGDSUiIkLmz5/vcdzhcMjYsWOlQYMGEh8fL926dZNdu3Z5tDl69Kj0799fatSoITVr1pRBgwZJVlaWR5tt27bJtddeK3FxcZKSkiJTpkwp69cIAADCPYA5deqUXHrppfLKK68Ue1wDjb/85S/y2muvyfr16yUhIUHS0tLk9OnTrjYavHz99deydOlSWbhwoQmKBg8e7Dp+4sQJ6d69uzRu3Fg2b94sL7zwgqSnp8vrr79e1q8TAACE80y8PXr0MFtxNPvy0ksvyVNPPSU33XSTeex//ud/pF69eiZTc8cdd8g333wjixcvlo0bN0rHjh1Nm5dffll69uwpU6dONZmdWbNmSW5urrz99tsSGxsrrVu3lq1bt8q0adM8Ah0AABCeAloDs2fPHjl48KDpNnJKSkqS1NRUWbt2rbmvt9pt5AxelLaPjIw0GRtnm86dO5vgxUmzODt37pRjx44V+945OTkmc+O+AQCA0BTQAEaDF6UZF3d633lMb+vWretxPDo6WmrXru3RprjXcH+Pc02aNMkES85N62aAcsk9JZKeVLTpPgCgygiZUUhjxoyRzMxM1/bTTz9V9ikBAAArBDD169c3t4cOHfJ4XO87j+nt4cOHPY7n5+ebkUnubYp7Dff3OFe1atXMqCb3DQAAhKaABjBNmzY1Acby5ctdj2ktita2dOrUydzX2+PHj5vRRU6fffaZFBYWmloZZxsdmZSXl+dqoyOWWrRoIbVq1QrkKQMAgHAYhaTztezevdujcFdHCGkNS6NGjWTYsGHy7LPPysUXX2wCmqefftqMLOrTp49pf8kll8gNN9wg999/vxlqrUHK0KFDzQglbaf69esn48aNM/PDjB49WrZv3y4zZsyQ6dOnB/JrR7jLs3uvbcm1F7/vixibSERE2c8NABDYAGbTpk1y/fXXu+6PGDHC3A4cOFBmzpwpf/7zn81cMTrcWTMt11xzjRk2rRPSOekwaQ1aunbtakYf9e3b18wd46RFuEuWLJEhQ4ZIhw4dpE6dOmZyPIZQI5Di/3a5jv33rfHUi/x78Sf2i8QmlOm8AACli3Do5C0hSLuuNBDSgl4r1MMcOX5Quiz4ndlfedNSSa5ZfK0PAned1+/9SWzB+vEngAGAoH5++52BASxNu3bOyH70W7HZEktuq91GzszLyN0isWefW2p7AEBQEcAgvLjXpWgXj6/dPBq80CUEAFVGyMwDAwAAwgcBDAAAsBwCGAAAYDnUwAAl0ZqX9EyuDwBUQWRggsieZ5e277Y1m+4DAIDAIIABAACWQwADAAAshwAGAABYDkW8ftBVF7Lzs31uf8Se5bnvZSLX7PzT/pwKAABhjQDGDxq8pM5O9bm9ozBaIs7kuHrMTZOIyHy/v0EAAODX6EICAACWQwamjFbevlLio+O9ttFuo57zu5r9f93yqSTbqpfYNvvUCYmf0dLsO6LiynpaAACEBQIYP2tgzt6JLdq8PiHGc99Le0dhrNjOvL7dfcFBAADwKwQwfjidX+ja7zBhWekBTESuJBYlVeTa51d6bR8vp+UbEi8AAPiEGhigBPbcfGny+Cdm030AQNVBBqaMPh99vdeaFnXEflJ6LnC27yLJtsSSG+eeEplatBsfE1XW0wIAICwQwJSRLTZKbLHeL192/tnj2tZ7+7PHIqiBqTLD5hMvefzMfmexxXoJQAEAFYoAJoh0lNLJbya79lG12HMLJD665K4h924j3ffW1sjNd81VqAXflGIDQPAQwCBsXfv8Cu+F2H4UYZ9biJ2dVyC2aoE8WwCAOwKYINIuo72TewXzLVAOEZG54jg7sKwYuZ77paZUcl1D4D2G3AMAAo4ABmElLvrswLvqzZ8tdSmIs20n+7QURKqkmNuVBacloVxnCgDwhmHUCCsUSANAaCADg7CixdTr+633qa0/S0Goo5mHpccnvQNyngAA7whgEHYZGFuMc6yQd9kxBa59W0x8qc/LjmYqZQCoKHQhAQAAyyGAAUrgPiMysyMDQNVCFxJQAu0yyhiYwfUBgCqIDAwAALAcAhgAAGA5BDBAZdDVx9OTijbdBwD4hQAGAABYDgEMAACwHAIYAABgOQyjBoIhz+69tiXXXvx+aXQ24DMrXgNAOCOAAYIg/m+XizgcvjWeepHvL/zEfpFY1rkGgIB3ITVp0sSsN3PuNmTIEHO8S5cuvzr2wAMPeLzGvn37pFevXmKz2aRu3boyatQoyc/P57sFAACCk4HZuHGjFBScXQRv+/bt8rvf/U5uu+0212P333+/jB8/3nVfAxUnfa4GL/Xr15c1a9bIgQMHZMCAARITEyMTJ04M9OkCgeO22GP2o9+KzZZYclvtNnJmXkbuFom1+dYWABCcAOa8887zuD958mRp1qyZXHfddR4BiwYoxVmyZIns2LFDli1bJvXq1ZP27dvLhAkTZPTo0ZKeni6xsbHFPi8nJ8dsTidOnAjY1wT4xL02Rbt5fO3q0eCFbiEAqDqjkHJzc+Uf//iH3HfffaaryGnWrFlSp04dadOmjYwZM0bs9rNFjGvXrpW2bdua4MUpLS3NBCRff/11ie81adIkSUpKcm0pKSlB/MoAAEDIFvHOnz9fjh8/Lvfcc4/rsX79+knjxo2lYcOGsm3bNpNZ2blzp8ydO9ccP3jwoEfwopz39VhJNBAaMWKE674GPAQxAACEpqAGMG+99Zb06NHDBCtOgwcPdu1rpqVBgwbStWtX+f77701XU1lVq1bNbIAlaJdRemZlnwUAWFbQupB+/PFHU8fyxz/+0Wu71NRUc7t7925zq7Uxhw4d8mjjvF9S3QwAAAgvQQtg3nnnHTMEWkcUebN161Zzq5kY1alTJ8nIyJDDhw+72ixdulRq1KghrVq1CtbpAgCAcO9CKiwsNAHMwIEDJTr67FtoN9Hs2bOlZ8+ekpycbGpghg8fLp07d5Z27dqZNt27dzeByt133y1TpkwxdS9PPfWUmUeGLiIAABC0AEa7jnQyOh195E6HQOuxl156SU6dOmWKbPv27WsCFKeoqChZuHChPPjggyYbk5CQYAIh93ljAABAeAtKAKNZFEcx06hrwLJq1apSn6+jlBYtWhSMUwMAACGA1agBAIDlEMAAAADLIYABAACWQwADAAAshwAmmHJPiaQnFW26DwAAAoIABgAAWA4BDAAAsJygLuYYctznttEuoehS4r9ce/H7pbUFAABeEcD4I+9skBE/o6VnQFOaqRf59VYAAKBkdCEBAADLIQNTRtkPbRFbUr3Su4WcmZeRu0Vibb69eIyP7QAACFMEMGWlQUZsgu/tY/1sDwAASkQXEgAAsBwCGAAAYDl0IQWTdhmlZwb1LWBN9jy7pM5ONfvr+60XG3VPAOAXMjBAJcjOKyh2HwDgGwIYAABgOXQhAUGQnZ8t9ryoEo/b87I99uO9tJX8bJGICLMb73BI0R4AhDcCGCAIeszr6vW4ozBaIs7kP3vMTZOIyHzvL9gkxdysy8+WhGrVS26nS1xMbFi0/8R+hu4DCFl0IQEWcpp6GQAwyMAAARIXFScnvx3vY+tcSWz5rNnL+u5xHbJWYsv4yBMS3XxqgM4SAEIDAQwQILbYaNkxrrdPbY/YT0rPBUUBzOeju0uyLbHktscPSc9/TT27oKh2EwViBfRz6VDuM7U2AFDVEcAAARIREWGCGF9k559tp8/x9rzs2LMFvvF/u9z3VdD9XQGdmhkAFkIAA1SC+Oh4OfnNZNc+AMA/BDBAJdCMy97JvXxr7DZLb/aj34rNS3eT3yugu7cHAAshgAGqOve6FF2ewtdVzVkBHUAIYxg1AACwHDIwQChhAVEAYYIMDAAAsBwCGAAAYDkEMAAAwHIIYAAAgOUQwAAAAMshgAEAAJZDAAMAACyHAAYAAFhOwAOY9PR0syqv+9ayZUvX8dOnT8uQIUMkOTlZqlevLn379pVDhw55vMa+ffukV69eYrPZpG7dujJq1CjJz88P9KkCAACLCspMvK1bt5Zly5adfZPos28zfPhw+eSTT+TDDz+UpKQkGTp0qNxyyy3y5ZdfmuMFBQUmeKlfv76sWbNGDhw4IAMGDJCYmBiZOHFiME4XAABYTFACGA1YNAA5V2Zmprz11lsye/Zs+e1vf2see+edd+SSSy6RdevWyVVXXSVLliyRHTt2mACoXr160r59e5kwYYKMHj3aZHdiY2OLfc+cnByzOZ04cSIYXxoAAAjVGphdu3ZJw4YN5cILL5T+/fubLiG1efNmycvLk27durnaavdSo0aNZO3atea+3rZt29YEL05paWkmIPn6669LfM9JkyaZjI5zS0lJCcaXBgAAQjGASU1NlZkzZ8rixYvl1VdflT179si1114rJ0+elIMHD5oMSs2aNT2eo8GKHlN66x68OI87j5VkzJgxJsPj3H766adAf2kAACBUu5B69Ojh2m/Xrp0JaBo3biwffPCBxMfHS7BUq1bNbAAAIPQFfRi1ZluaN28uu3fvNnUxubm5cvz4cY82OgrJWTOjt+eOSnLeL66uBkAFyT0lkp5UtOk+AIRyAJOVlSXff/+9NGjQQDp06GBGEy1fvtx1fOfOnaZGplOnTua+3mZkZMjhw4ddbZYuXSo1atSQVq1aBft0AQBAOHYhjRw5Unr37m26jfbv3y/PPPOMREVFyZ133mmKawcNGiQjRoyQ2rVrm6Dk4YcfNkGLjkBS3bt3N4HK3XffLVOmTDF1L0899ZSZO4YuIiCIcu2+Hy+trbsYm0hERNnPCwAqIoD5z3/+Y4KVI0eOyHnnnSfXXHONGSKt+2r69OkSGRlpJrDTYc86wuhvf/ub6/ka7CxcuFAefPBBE9gkJCTIwIEDZfz48YE+VQDupl4UnLZP7BeJTeBaA6jaAcycOXO8Ho+Li5NXXnnFbCXR7M2iRYsCfWoAACBEBGUiOwAWod07miHxhXYbOTMvI3eLxNp8awsAQUAAA4QzrU0pS/eOBi90CwGoRKxGDYQQe26+NHn8E7PpPgCEKjIwAHyjGZf0TK4WgCqBAAawEHtugcRHl5xZcc+6+JuBiY+JkgiGOwOwCAIYwEKufX6FiKP4FdnP1fHZsxNG+mLH+DSxxfIrAYA1UAMDhJKIXEm85HGz6T4AhCr+3AKquLjos39nfDHmaomPLnlR1GP2LLn5k6L9ZY9dJbVs1UvtkjJZHRFxOByBOmUACDoCGKCKc69L6TGvq9e2jsJoiTgT7/T5uKdERJZeB5PYsuj2dEE3SZCYcp4tAFQMupAAAIDlkIEBqjjtMlrfb71PbbPzs6XLB13M/qo7lnntblJHs7NKzeoAQFVEAANYoAvJplP++0mDl9Kel51XUI4zA4DKQwADhBANWDIGZlT2aQBA0FEDAwAALIcABgAAWA4BDAAAsBwCGAAAYDkEMAAAwHIIYAAAgOUQwACofLmnRNKTijbdB4BSEMAAAADLYSI7AMGVa/evjS/tnXSmYbfFLgGEDwIYAME19aLgtX9iv0hsgt+nBMD66EICAACWQwYGQOBp145mR3yl3UbOzMvI3SKxNt/aAghbBDAAfGLPzZdWYz81+zvGp4kt1suvD61LKWvXjgYvdAsBKAUBDADDnlsg8dH5XgOY4vZLEx8TJRGlFdpqwJKeyXcCgM8IYAAY1z6/QsQR69PV6Pjscp+vWqnZGgAoA4p4AfgmIlcSL3ncbLoPAJWJP4uAMBYXffZvmC/GXC3x0fEltj1mz5KbPynaX/bYVVLLVt1rd5TJ6IiIw+EI5CkDgEEAA4Qx99qUHvO6em3rKIyWiDPxTp+Pe0pEpPc6mMSWRbenC7pJgsQE4GwB4Cy6kAAAgOWQgQHCmHYZre+33qe22fnZ0uWDLmZ/1R3LvHY3Hc3OKjWjAwDlQQADhHkXkk0nnfOTBi/enpedV1DOMwMA7whgAPhEA5aMgRlcLQBVAjUwAADAcgIewEyaNEmuuOIKSUxMlLp160qfPn1k586dHm26dOliUtfu2wMPPODRZt++fdKrVy+x2WzmdUaNGiX5+b7P/gkAAEJXwLuQVq1aJUOGDDFBjAYcTzzxhHTv3l127NghCQln10a5//77Zfz48a77Gqg4FRQUmOClfv36smbNGjlw4IAMGDBAYmJiZOLEiYE+ZQAAEO4BzOLFiz3uz5w502RQNm/eLJ07d/YIWDRAKc6SJUtMwLNs2TKpV6+etG/fXiZMmCCjR4+W9PR0iY31bbpzAGEu95TIxIZF+7o6NotEAiEj6DUwmZlFC7TVrl3b4/FZs2ZJnTp1pE2bNjJmzBix2+2uY2vXrpW2bdua4MUpLS1NTpw4IV9//XWx75OTk2OOu28AQlyuvShIKXGz+9HWbWP2YCC8RyEVFhbKsGHD5OqrrzaBilO/fv2kcePG0rBhQ9m2bZvJrGidzNy5c83xgwcPegQvynlfj5VUezNu3LhgfjkAqpqpFwWnLdkaILwDGK2F2b59u3zxxRcejw8ePNi1r5mWBg0aSNeuXeX777+XZs2alem9NIszYsQI133NwKSkpJTj7AEAQNgFMEOHDpWFCxfK6tWr5YILLvDaNjU11dzu3r3bBDBaG7NhwwaPNocOHTK3JdXNVKtWzWwAQpxOoKcZEl9ot5Ez8zJyt0iszbe2AMKvBkZXntXgZd68efLZZ59J06ZNS33O1q1bza1mYlSnTp0kIyNDDh8+7GqzdOlSqVGjhrRq1SrQpwzASnQBSi3G9WlzC1h039e2AMIvA6PdRrNnz5YFCxaYuWCcNStJSUkSHx9vuon0eM+ePSU5OdnUwAwfPtyMUGrXrp1pq8OuNVC5++67ZcqUKeY1nnrqKfPaZFkA+EwDk/SigQQAQkvAMzCvvvqqGXmkk9VpRsW5/fOf/zTHdQi0Do/WIKVly5by2GOPSd++feXjjz92vUZUVJTpftJbzcbcddddZh4Y93ljAABA+IoORheSN1pYq5PdlUZHKS1atCiAZwYAAEIFizkCCKrs/Gyx50WV2qbLB13M/srbV5rVrn2h7XQpEgDhhwAGQFD1mNfVr/bOQMYX6/utN6tkAwg/rEYNoNI5CqOL3QeAkvCbAkDAxUXFyclvi4ruPx99vdhivXchHbNnyc2fdDf783svklq26l67m5xZndJq7gCELgIYAAFn6lIcRYuuXjv5Sx+fNdn82+2bLaW8eK4ktizaPZ1fKAnBWNvVfQ2lQEyS5067vKjbAcqNAAYAzuXvjLysswRUOAIYAAEXHxMlO8an+dzenpsvHZ9dbvY3PdVVbLEl/2o6Ys+SngsCcpoALIwABkBQupC8BSHe6PO8PTc733s9TYWssaRYZwmoVAQwACqdBix7J/cK+OtqZqfV2E/NvmaEvAZVzjWWysK5zhKACkMAA8CyjtqzSh3dlHjJ42Z//4nLvY5uclcrLkEiI0uZZYJ1loBKRQADwLJuOTP0uiQ6p0zEmTikz8c9JSIy36fXXXnbGkm2JUqlyD0lMrFh0b52aZHZAYpFAAMAFam0Idrux30Zzu3E8GyEGQIYAJai3TuaIfGFPxPkaXdUaRmdgPBnyDXDs4ESEcAAsBStTfG1e0fbZQzMCPo5Aah4BDAAEGz+DNFmeDbgEwIYAAg2f4ZoM7oJ8AkBDAAUs2CkPc/7hHlmRuD5RYtKLuqzXJJ9HKIdHx1ftFYUgHIhgAGAczhXu/Z1iHaPuWk+D9Fed+c6SWBoNFBuBDAAUIFYQRsIDAIYABCRuKg4OfnteN+vhUP/iT1zJ1fES69QRGSuVG/+bHCvMytoI8wQwADAmfWYdozrXcYVtHv4sIJ2kAMYIMwQwABAGVbQ9mcBSlbQBgKPAAYArMrfFbQZoo0QQgADAPDOnzWZ/MUaTigjAhgAQGALhP3BitsoozOzGAAAAFgHGRgAqED23AKJj873cXRTV58Li+NjogI7w2+w1m86tz1QRgQwAFCBOr/wqTgKnfPHeJ9fpuOz//I6v4y7r9NvlIRqMRIwrN+EKo4ABgAqUGkT2hUtUZD/q/3SHMu5WiIifVuPyV9BXb+ptAJhf7M7ThQHhzwCGAAIsrho38sN3QMWX4MXX9dvKqv1/daLTQOCYPCnK8mfthQHhzwCGAAIMv3w1yAg0I5mZwU1cHFyOEy/FlClEMAAQEXM8huEDIajMKZC1m8K+AKU/hQI+4Pi4LBCAAMAFmXqUhxljSxizwQ0xXMUStWZQbgsmHwv5BHAAIBF6dDpHePTgvLaRQtQ+jb0uzwCOfxbu7qyna/14sWlNRZnTsyUEftxDvEjvpOIal4CMAqPKwQBDACEyQKU/nBfgNKfod+ldU2da+MTPSShWmC+huzsLOnSJMWntvGFhbLhx/+Y/S6NL5DsSN8LrddPay42X+uCKDwOGgIYAEClDP1WV07ODVhwFBGZI9Wb+/a+GrC0bdpIKpt+ea6skf2ISH52UN4nPiouoEPhNds1b8FH0qVrmtROTpbKUKUDmFdeeUVeeOEFOXjwoFx66aXy8ssvy5VXXlnZpwUAIa8ihn4HMzia22uJ1LYFbl6c7Dy79Jjfzezbh20rKkQOALv9iFz/ya1Fdxb0kmBZ+eN/JN5b1sjPLrXthwqk76tZkpiYKI8++qgMHz5cateuLRWpygYw//znP2XEiBHy2muvSWpqqrz00kuSlpYmO3fulLp161b26QFASAvW0G9lz7PL9R9eH9Tg6HcvritHgXNxJ5IriS2Ldq+fH7xAw5u4wkI5faarK85t3xfaTRZI2aKZoizp3r27TJs2TWbMmFHhgUyEo4oO8Neg5YorrpC//vWv5n5hYaGkpKTIww8/LI8//nipzz9x4oQkJSVJZmam1KhRIyDndOT4Qemy4Hdmf+VNSyW5Zv2AvC4AhBNTbBuErpLsvALpMGHZmTfRZRUiAhzAjPV5eHvWzglmv3qLpyUiMs+n5+V/N1LES3eaflqflqLPszg5UXrdcWSuRDef6uM5+5fpyt6bLd+nfy+bN2+WCy64QKZOnWp6TaKiosodyPj6+V0lMzC5ubnmoowZM8b1WGRkpHTr1k3Wrl1b7HNycnLM5qRfuPNCBMrJEyelILvAtR8TGaSZKQEAfot2OGTdY9cF5crZc/Kly4ul//F81nHz74mMx3x+xsrHfie2ABU0nz3nwH/Mx0WelILY5133tVdkypQpMnLkSBPIODMyDzzwgEk6aDDiD+fndqn5FUcV9PPPP+tZO9asWePx+KhRoxxXXnllsc955plnzHPYuAb8DPAzwM8APwP8DEiFXIPNmzf/6vM4IyPDcd5555X7tX/66SevsUKVzMCUhWZrtGbGSbucjh49KsnJyQGtvNbIULuyfvrpp4B1TYFrXZn4meY6hxJ+nivGl19+KT179vR47PDhw2bgzd/+9jfTlTRq1KgyZWA083Ly5Elp2LCh13ZVMoCpU6eO+eIPHTrk8bjer1+/+LqTatWqmc1dzZo1g3aOGrwQwFQMrjXXOZTw88x1DgUJCQklBi6aTChvMa8vQY/vJcwVKDY2Vjp06CDLly/3yKjo/U6dOlXquQEAgCITJ06Upk2byt///ncTuOzdu1cmTJhQISORqmQGRumFGDhwoHTs2NHM/aLDqE+dOiX33ntvZZ8aAABhLSZGR3mJLFmyJCAZl5AKYP7whz/IL7/8ImPHjjUT2bVv314WL14s9erVq9Tz0m6qZ5555lfdVeBaWxU/01znUMLPc8W47LLL5PbbbzfJhQYNGkhlqLLzwAAAAFiqBgYAAMAbAhgAAGA5BDAAAMByCGAAAIDlEMD4SReratKkicTFxZkFJzds2BCc70yYmDRpklm0U5dk1/U0+vTpY1Ycd3f69GkZMmSImVW5evXq0rdv319Ncgj/TJ482cxQPWzYMK5zgP38889y1113mZ/X+Ph4adu2rWzatMl1XMdN6OhKHbmhx3WNt127dgX6NEJaQUGBPP3002b+Eb2GzZo1M3OPuI9J4TqXzerVq6V3795mFlz9HTF//nyP475cV50Fv3///mbSRp1QdtCgQZKVlSUBV75Vi8LLnDlzHLGxsY63337b8fXXXzvuv/9+R82aNR2HDh2q7FOzrLS0NMc777zj2L59u2Pr1q2Onj17Oho1auTIyspytXnggQccKSkpjuXLlzs2bdrkuOqqqxy/+c1vKvW8rWzDhg2OJk2aONq1a+d49NFHXY9zncvv6NGjjsaNGzvuuecex/r16x0//PCD49NPP3Xs3r3b1Wby5MmOpKQkx/z58x1fffWV4/e//72jadOmjuzs7ACcQXh47rnnHMnJyY6FCxc69uzZ4/jwww8d1atXd8yYMcPVhutcNosWLXI8+eSTjrlz55r1iObNm+dx3JfresMNNzguvfRSx7p16xyff/6546KLLnLceeedjkAjgPGDLiQ5ZMgQ1/2CggJHw4YNHZMmTQr4NyZcHT582PynWbVqlbl//PhxR0xMjPkF5fTNN9+YNmvXrq3EM7WmkydPOi6++GLH0qVLHdddd50rgOE6B8bo0aMd11xzTYnHCwsLHfXr13e88MILrsf02lerVs3x/vvvB+gsQl+vXr0c9913n8djt9xyi6N///5mn+scGOcGML5c1x07dpjnbdy40dXmX//6lyMiIsIs1BxIdCH5KDc3VzZv3mzSZU6RkZHm/tq1awOfGgtTmZmZ5tY5o6Ne87y8PI/r3rJlS2nUqBHXvQy0K65Xr14e15PrHDgfffSRmT38tttuM12iOtnXG2+84Tq+Z88eMzGn+/XXNV+0O5rfI777zW9+Y5aW+e6778z9r776Sr744gvp0aMH1zmIfPn51VvtNtL/B07aXj8v169fHx4z8VY1//3vf02/67kzAev9b7/9ttLOK5Toeldak3H11VdLmzZtzGP6n0XXxjp3YU697noMvpszZ45s2bJFNm7c+KtjXOfA+OGHH+TVV181U6s/8cQT5lo/8sgj5mdYl0Zx/swW93uEn2ffPf7442bVaf1jRhcP1N/Nzz33nKm7cP48c50Dz5frqrcavLuLjo42f5QG+mecAAZVKjuwfft285cUAuunn36SRx99VJYuXWoK0BG8IFz/8tQF7pRmYPRn+rXXXjMBDALjgw8+kFmzZsns2bOldevWsnXrVvPHjxaecp3DB11IPqpTp46J9M8d/aL369evH4zvTVgZOnSoLFy4UFasWCEXXHCB63G9ttp9d/z4cY/2XHf/aFecLnl/+eWXm7+GdFu1apX85S9/Mfv6FxTXufx0ZEarVq08Hrvkkktk3759Zt/5u4LfI+UzatQok4W54447zCivu+++2ywmqKMauc7B48vPr97q7xp3+fn5ZmRSoD8rCWB8pCngDh06mH5X97+29H6nTp0C+k0JJ1onpsHLvHnz5LPPPjPDIt3pNddVT92vuw6z1g8ErrvvunbtKhkZGeYvVeemmQJNuTv3uc7lp92f504DoHUajRs3Nvv6862/xN1/nrUrRGsD+Hn2nd1uNzUV7vQPTP2dzHUOHl9+fvVW/+DUP5qc9He7fm+0ViagAloSHAbDqLXaeubMmabSevDgwWYY9cGDByv71CzrwQcfNEPyVq5c6Thw4IBrs9vtHsN7dWj1Z599ZoZRd+rUyWwoH/dRSFznwA1Rj46ONsN8d+3a5Zg1a5bDZrM5/vGPf3gMQ9XfGwsWLHBs27bNcdNNNzGM2k8DBw50nH/++a5h1Drkt06dOo4///nPXOcAjFT897//bTYNEaZNm2b2f/zxR59/fnUY9WWXXWamEvjiiy/MyEeGUVcBL7/8svkw1flgdFi1jnNH2el/kOI2nRvGSf9jPPTQQ45atWqZD4Obb77ZBDkIbADDdQ6Mjz/+2NGmTRvzx07Lli0dr7/+usdxHYr69NNPO+rVq2fadO3a1bFz584AvXt4OHHihPnZ1d/FcXFxjgsvvNDMXZKTk+Nqw3UumxUrVhT7O1mDRl+v65EjR0zAonPz1KhRw3HvvfeawCjQIvSfwOZ0AAAAgosaGAAAYDkEMAAAwHIIYAAAgOUQwAAAAMshgAEAAJZDAAMAACyHAAYAAFgOAQwAALAcAhgAAGA5BDAALK1Jkyby0ksvVfZpAKhgBDAAAMByWAsJQJXWpUsXadOmjdl/7733JCYmRh588EEZP368XH/99bJq1SqP9izvBoQHMjAAqrx3331XoqOjZcOGDTJjxgyZNm2avPnmmzJ37ly54IILTDBz4MABswEID9GVfQIAUJqUlBSZPn26RERESIsWLSQjI8Pcv//++yUqKkoSExOlfv36XEggjJCBAVDlXXXVVSZ4cerUqZPs2rVLCgoKKvW8AFQeAhgAAGA5BDAAqrz169d73F+3bp1cfPHFpvsoNjaWTAwQhghgAFR5+/btkxEjRsjOnTvl/fffl5dfflkeffRR1zwwq1evlp9//ln++9//VvapAqggDKMGUOWHUbdu3VoKCwtl9uzZJuuiw6ifffZZUxej2Zg//elPJrjJyclhGDUQJghgAFT5AKZ9+/bMtgvAA11IAADAcghgAACA5dCFBAAALIcMDAAAsBwCGAAAYDkEMAAAwHIIYAAAgOUQwAAAAMshgAEAAJZDAAMAACyHAAYAAIjV/D/niEWBpJ4ZZAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "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(dataset=\"gen2rwt\", pt=ptvals, eta=etavals, weight=corr(ptvals, etavals))\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://coffea-hep.readthedocs.io/en/latest/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": { "tags": [] }, "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 numpy as np\n", "\n", "ext = extractor()\n", "ext.add_weight_sets(\n", " [\n", " \"* * data/Fall17_17Nov2017_V32_MC_L2Relative_AK4PFPuppi.jec.txt\",\n", " \"* * data/Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi.junc.txt\",\n", " ]\n", ")\n", "ext.finalize()\n", "\n", "jec_stack_names = [\"Fall17_17Nov2017_V32_MC_L2Relative_AK4PFPuppi\", \"Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi\"]\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": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "starting columns: {'nConstituents', 'rawFactor', 'electronIdx2', 'nElectrons', 'btagDeepFlavC', 'Rho', 'pt_gen', 'chHEF', 'bRegCorr', 'btagDeepFlavB', 'partonFlavour', 'genJetIdxG', 'genJetIdx', 'muonIdx2', 'btagDeepC', 'electronIdxG', 'electronIdx1G', 'pt_raw', 'jercCHF', 'muonIdx1', 'electronIdx1', 'mass', 'cleanmask', 'electronIdx2G', 'neEmEF', 'qgl', 'area', 'muonIdx2G', 'pt', 'muonIdxG', 'chEmEF', 'bRegRes', 'btagDeepB', 'hadronFlavour', 'jetId', 'btagCSVV2', 'btagCMVA', 'muonSubtrFactor', 'mass_raw', 'phi', 'nMuons', 'muEF', 'jercCHPUF', 'neHEF', 'puId', 'eta', 'muonIdx1G'}\n", "new columns: {'jet_energy_uncertainty_jes', 'pt_jec', 'pt_orig', 'JES_jes', 'jet_energy_correction', 'mass_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(Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi=evaluator[\"Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi\"])\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": { "tags": [] }, "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)\n", "print(\"untransformed mass ratios\", jets.mass / jets.mass_raw)\n", "\n", "print(\"transformed pt ratios\", corrected_jets.pt / corrected_jets.pt_raw)\n", "print(\"transformed mass ratios\", corrected_jets.mass / corrected_jets.mass_raw)\n", "\n", "print(\"JES UP pt ratio\", corrected_jets.JES_jes.up.pt / corrected_jets.pt_raw)\n", "print(\"JES DOWN pt ratio\", corrected_jets.JES_jes.down.pt / corrected_jets.pt_raw)" ] }, { "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://coffea-hep.readthedocs.io/en/latest/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": { "tags": [] }, "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))\n", "print(\"systematic +:\", btag_sf.eval(\"up\", events.Jet.hadronFlavour, abs(events.Jet.eta), events.Jet.pt))\n", "print(\"systematic -:\", btag_sf.eval(\"down\", events.Jet.hadronFlavour, abs(events.Jet.eta), events.Jet.pt))" ] }, { "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": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
CorrectionSet (schema v2)\n",
       "my custom corrections\n",
       "📂\n",
       "└── 📈 gen2_to_gen1 (v0)\n",
       "    Reweights gen2 to agree with gen1\n",
       "    Node counts: MultiBinning: 1\n",
       "    ╭──────────── ▶ input ─────────────╮ ╭──────────── ▶ input ────────────╮\n",
       "    │ pt (real)                        │ │ eta (real)                      │\n",
       "    │ pt                               │ │ eta                             │\n",
       "    │ Range: [0.0, 100.0), overflow ok │ │ Range: [-3.0, 3.0), overflow ok │\n",
       "    ╰──────────────────────────────────╯ ╰─────────────────────────────────╯\n",
       "    ╭─── ◀ output ───╮\n",
       "    │ out (real)     │\n",
       "    │ No description │\n",
       "    ╰────────────────╯\n",
       "
\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\n", "import 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.model_dump_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": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([1.34972203, 0.98072826, 1.34972203, ..., 1.20142093, 0.51586373,\n", " 1.32139074], shape=(10000,))" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ceval = cset.to_evaluator()\n", "\n", "ceval[\"gen2_to_gen1\"].evaluate(ptvals, etavals)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
[[1, 0.514, 0.787, 0.91, 0.91],\n",
       " [0.398, 0.765, 0.981, 0.941, 0.91, 0.941, 1.11, 0.91],\n",
       " [1, 0.505, 0.738, 0.646, 0.854],\n",
       " [0.621, 0.613, 0.941],\n",
       " [0.426, 0.271, 0.787, 0.941, 0.941],\n",
       " [0.698, 0.765, 0.774, 0.88, 0.906, 0.941, 0.91, 1.2],\n",
       " [1, 0.777, 0.613, 0.491],\n",
       " [0.613, 0.491, 0.906, 0.941],\n",
       " [0.906],\n",
       " [0.994, 0.609, 0.765, 0.774, 0.491, 0.787, 0.91, 1.11, 0.91],\n",
       " ...,\n",
       " [0.861, 0.88],\n",
       " [0.906, 0.91, 0.91, 0.91],\n",
       " [0.514, 0.523, 0.491, 1.14, 1.2, 0.941, 1.11, 1.11, 0.91],\n",
       " [0.791, 0.787, 0.88],\n",
       " [1.14, 1.2],\n",
       " [0.981, 0.854, 1.2],\n",
       " [1.07],\n",
       " [0.459, 0.646, 0.854, 0.906, 0.854, 1.2],\n",
       " [0.941, 1.11]]\n",
       "--------------------------------------------------------------\n",
       "backend: cpu\n",
       "nbytes: 1.8 kB\n",
       "type: 40 * var * float64
" ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def myJetSF(jets):\n", " # Older correctionlib did not handle jagged arrays\\\n", " # the solution for flattening and unflattening is left commented below\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", " return ceval[\"gen2_to_gen1\"].evaluate(jets.pt, jets.eta)\n", "\n", "\n", "myJetSF(events.Jet)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In dask mode, the only difference is that you need to call compute on stuff" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "fname = \"coffea/tests/samples/nano_dy.root\"\n", "events = NanoEventsFactory.from_root(\n", " {fname: \"Events\"},\n", " schemaclass=NanoAODSchema,\n", " metadata={\"dataset\": \"DYJets\"},\n", " mode=\"dask\",\n", ").events()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "starting columns: {'nConstituents', 'rawFactor', 'electronIdx2', 'nElectrons', 'btagDeepFlavC', 'Rho', 'pt_gen', 'chHEF', 'bRegCorr', 'btagDeepFlavB', 'partonFlavour', 'genJetIdxG', 'genJetIdx', 'muonIdx2', 'btagDeepC', 'electronIdxG', 'electronIdx1G', 'pt_raw', 'jercCHF', 'muonIdx1', 'electronIdx1', 'mass', 'cleanmask', 'electronIdx2G', 'neEmEF', 'qgl', 'area', 'muonIdx2G', 'pt', 'muonIdxG', 'chEmEF', 'bRegRes', 'btagDeepB', 'hadronFlavour', 'jetId', 'btagCSVV2', 'btagCMVA', 'muonSubtrFactor', 'mass_raw', 'phi', 'nMuons', 'muEF', 'jercCHPUF', 'neHEF', 'puId', 'eta', 'muonIdx1G'}\n", "new columns: {'jet_energy_uncertainty_jes', 'pt_jec', 'pt_orig', 'JES_jes', 'jet_energy_correction', 'mass_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(Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi=evaluator[\"Fall17_17Nov2017_V32_MC_Uncertainty_AK4PFPuppi\"])\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": "code", "execution_count": 18, "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": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dask.awkward" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "myJetSF(events.Jet)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
[[1, 0.514, 0.787, 0.91, 0.91],\n",
       " [0.398, 0.765, 0.981, 0.941, 0.91, 0.941, 1.11, 0.91],\n",
       " [1, 0.505, 0.738, 0.646, 0.854],\n",
       " [0.621, 0.613, 0.941],\n",
       " [0.426, 0.271, 0.787, 0.941, 0.941],\n",
       " [0.698, 0.765, 0.774, 0.88, 0.906, 0.941, 0.91, 1.2],\n",
       " [1, 0.777, 0.613, 0.491],\n",
       " [0.613, 0.491, 0.906, 0.941],\n",
       " [0.906],\n",
       " [0.994, 0.609, 0.765, 0.774, 0.491, 0.787, 0.91, 1.11, 0.91],\n",
       " ...,\n",
       " [0.861, 0.88],\n",
       " [0.906, 0.91, 0.91, 0.91],\n",
       " [0.514, 0.523, 0.491, 1.14, 1.2, 0.941, 1.11, 1.11, 0.91],\n",
       " [0.791, 0.787, 0.88],\n",
       " [1.14, 1.2],\n",
       " [0.981, 0.854, 1.2],\n",
       " [1.07],\n",
       " [0.459, 0.646, 0.854, 0.906, 0.854, 1.2],\n",
       " [0.941, 1.11]]\n",
       "--------------------------------------------------------------\n",
       "backend: cpu\n",
       "nbytes: 1.8 kB\n",
       "type: 40 * var * float64
" ], "text/plain": [ "" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "myJetSF(events.Jet).compute()" ] } ], "metadata": { "kernelspec": { "display_name": "awkward-coffea", "language": "python", "name": "awkward-coffea" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.9" } }, "nbformat": 4, "nbformat_minor": 4 }