L2L3_profiles_demo_HATS21.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# JEC profile plots\n",
    "\n",
    "Coffea 0.7.1 https://github.com/cms-jet/CoffeaJERC/blob/master/genL2L3.ipynb \n",
    "\n",
    "and [nanoevents.ipynb](https://github.com/CoffeaTeam/coffea/blob/master/binder/nanoevents.ipynb) to illustrate profile plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/cms.mmorris/coffeajerc/CoffeaJERC/coffea-jerc/lib/python3.6/site-packages/coffea/nanoevents/schemas/nanoaod.py:195: RuntimeWarning: Missing cross-reference index for FatJet_genJetAK8Idx => GenJetAK8\n",
      "  RuntimeWarning,\n"
     ]
    }
   ],
   "source": [
    "import awkward as ak\n",
    "import numpy as np\n",
    "import time\n",
    "from coffea.nanoevents import NanoEventsFactory, NanoAODSchema\n",
    "fname = \"root://cmseos.fnal.gov//store/user/hats/2021/JEC/RunIISummer19UL17NanoAOD_QCD_Pt-15to7000_NANOAODSIM_JMECustomTuples.root\"\n",
    "\n",
    "events = NanoEventsFactory.from_root(fname, schemaclass=NanoAODSchema).events()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from coffea import processor, hist\n",
    "class FancyJECL2L3Processor(processor.ProcessorABC):\n",
    "    def __init__(self):\n",
    "        dataset_axis = hist.Cat(\"dataset\", \"Primary dataset\")\n",
    "        eta_axis = hist.Bin(\"eta\", r\"$\\eta$\", 20, -5, 5)\n",
    "        pt_axis = hist.Bin(\"pt\", r\"$p_{T}$ [GeV]\", \n",
    "                           np.array([5,10,15,20,25,30,35,40,45,50,60,70,80,90,\n",
    "                                     100,120,140,160,180,\n",
    "                                     200,250,300,350,400,450,500,\n",
    "                                     600,700,800,900,1000,\n",
    "                                     1500,2000,3000,4000,5000]))\n",
    "        dr_axis = hist.Bin(\"dr\", r\"$\\delta (\\eta)$\", 20, 0., 1)\n",
    "        m_axis = hist.Bin(\"m\", r\"$p_{T}$ [GeV]\", 200, 0, 500)\n",
    "        r_axis = hist.Bin(\"r\", \"RECO / GEN response\", 200, 0, 5)\n",
    "        \n",
    "        self._accumulator = processor.dict_accumulator({\n",
    "            'pt':hist.Hist(\"Counts\", dataset_axis, pt_axis),\n",
    "            'eta':hist.Hist(\"Counts\", dataset_axis, eta_axis),\n",
    "            'dr':hist.Hist(\"Counts\", dataset_axis, dr_axis),\n",
    "            'r_pt_ptveta':hist.Hist(\"Counts\", dataset_axis, pt_axis, eta_axis, r_axis),\n",
    "            'r_m_ptveta':hist.Hist(\"Counts\", dataset_axis, pt_axis, eta_axis, r_axis),\n",
    "            'r_m_ptvm':hist.Hist(\"Counts\", dataset_axis, pt_axis, m_axis, r_axis),\n",
    "            'cutflow': processor.defaultdict_accumulator(int),\n",
    "        })\n",
    "    \n",
    "    @property\n",
    "    def accumulator(self):\n",
    "        return self._accumulator\n",
    "    \n",
    "    def process(self, events):\n",
    "        output = self.accumulator.identity()\n",
    "        output['cutflow']['all events'] += len(events)\n",
    "        #print(len(events))\n",
    "        \n",
    "        selectedEvents = events[\n",
    "            (ak.num(events.Jet) > 2)\n",
    "        ]\n",
    "\n",
    "        #leading jet\n",
    "        #print(selectedEvents.Jet[:,0])\n",
    "        #leading two jets\n",
    "        jet = selectedEvents.Jet[:,0:2]\n",
    "        jet = ak.flatten(jet)\n",
    "        #only with genmatch\n",
    "        jet = jet[~ak.is_none(jet.matched_gen)]\n",
    "        #only with good deltaR match\n",
    "        jet = jet[jet.delta_r(jet.matched_gen)<0.2]\n",
    "        \n",
    "        ptresponse = jet.pt/jet.matched_gen.pt\n",
    "        \n",
    "        output['dr'].fill(dataset=selectedEvents.metadata[\"dataset\"],\n",
    "                            dr=jet.delta_r(jet.matched_gen))\n",
    "        output['pt'].fill(dataset=selectedEvents.metadata[\"dataset\"],\n",
    "                            pt=jet.pt)\n",
    "        output['eta'].fill(dataset=selectedEvents.metadata[\"dataset\"], \n",
    "                                 eta=jet.eta)\n",
    "        output['r_pt_ptveta'].fill( dataset=selectedEvents.metadata[\"dataset\"], pt=jet.pt, eta=jet.eta, r=ptresponse)\n",
    "        \n",
    "        return output\n",
    "\n",
    "    def postprocess(self, accumulator):\n",
    "        return accumulator\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6656c9e7c560457daff02ef99fbe90df",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value='Preprocessing'), FloatProgress(value=0.0, max=1.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2643ee4ae21b4e27b9f0584715bc8980",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value='Processing'), FloatProgress(value=0.0, max=6.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/cms.mmorris/coffeajerc/CoffeaJERC/coffea-jerc/lib/python3.6/site-packages/coffea/nanoevents/schemas/nanoaod.py:195: RuntimeWarning: Missing cross-reference index for FatJet_genJetAK8Idx => GenJetAK8\n",
      "  RuntimeWarning,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "{'pt': <Hist (dataset,pt) instance at 0x7f6e4d730cc0>, 'eta': <Hist (dataset,eta) instance at 0x7f6e4d730d30>, 'dr': <Hist (dataset,dr) instance at 0x7f6e4d730c88>, 'r_pt_ptveta': <Hist (dataset,pt,eta,r) instance at 0x7f6e4cbb4ba8>, 'r_m_ptveta': <Hist (dataset,pt,eta,r) instance at 0x7f6e4cbb4828>, 'r_m_ptvm': <Hist (dataset,pt,m,r) instance at 0x7f6e4cbb4748>, 'cutflow': defaultdict(<class 'int'>, {'all events': 626100})}\n",
      "Events/s: 22049.368940829772\n"
     ]
    }
   ],
   "source": [
    "samples = {\n",
    "    \"QCD\": [fname]\n",
    "}\n",
    "tstart = time.time() \n",
    "output = processor.run_uproot_job(\n",
    "    samples,\n",
    "    \"Events\",\n",
    "    FancyJECL2L3Processor(),\n",
    "    processor.iterative_executor,\n",
    "    {\"schema\": NanoAODSchema},\n",
    ")\n",
    "elapsed = time.time() - tstart\n",
    "print(output)\n",
    "print(\"Events/s:\", output['cutflow']['all events']/elapsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " output: \n",
      "\n",
      "{'pt': <Hist (dataset,pt) instance at 0x7f6e4d730cc0>, 'eta': <Hist (dataset,eta) instance at 0x7f6e4d730d30>, 'dr': <Hist (dataset,dr) instance at 0x7f6e4d730c88>, 'r_pt_ptveta': <Hist (dataset,pt,eta,r) instance at 0x7f6e4cbb4ba8>, 'r_m_ptveta': <Hist (dataset,pt,eta,r) instance at 0x7f6e4cbb4828>, 'r_m_ptvm': <Hist (dataset,pt,m,r) instance at 0x7f6e4cbb4748>, 'cutflow': defaultdict(<class 'int'>, {'all events': 626100})}\n",
      "\n",
      "histogram and dictionary names: \n",
      "\n",
      "pt\n",
      "eta\n",
      "dr\n",
      "r_pt_ptveta\n",
      "r_m_ptveta\n",
      "r_m_ptvm\n",
      "cutflow\n"
     ]
    }
   ],
   "source": [
    "print(\"\\n output: \\n\")\n",
    "print(output)\n",
    "print(\"\\nhistogram and dictionary names: \\n\")\n",
    "for key in output.keys():\n",
    "    print(key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mplhep as hep\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "EtaBins = output['r_pt_ptveta'].axis('eta')\n",
    "EtaBinNums = len(output['eta'].values()[('QCD',)])\n",
    "\n",
    "PtBins = output['r_pt_ptveta'].axis('pt')\n",
    "PtBinNums = len(output['pt'].values()[('QCD',)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot $\\eta$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEMCAYAAADeYiHoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfKklEQVR4nO3de3RV9d3n8feHoERRlCqrD3KRaLEYCwaIMN5Q643eiOOoIHUUR0t1pLWL2j507IOXts9qq1WrpQO0D0tsdZB6qbEPHaZWsHapGJCLTZSCeAtSS8FlpKAY/M4f5wSOyU44kOycXD6vtbI8+7d/v803p13nk71/Z/+2IgIzM7PGehS6ADMz65gcEGZmlsgBYWZmiRwQZmaWyAFhZmaJHBBmZpYo1YCQNF7SWknrJc1ood9/kxSSynPavpMdt1bS+WnWaWZmTfVM68CSioBZwLlALVAlqTIiahr1OxS4HliW01YKTAJOAI4CnpB0XETsSqteMzP7uDTPIMYA6yNiQ0TsBBYAFQn9vgf8CHg/p60CWBARH0TEq8D67PHMzKydpHYGAQwA3szZrgXG5naQNAoYFBH/KelbjcY+12jsgMb/gKSpwFSA3r17jx42bFgblW5m1j2sWLHiHxHRL2lfmgHRIkk9gDuAKft7jIiYC8wFKC8vj+XLl7dNcWZm3YSk15vbl2ZAbAQG5WwPzLY1OBT4DLBUEsC/AJWSJuQx1szMUpbmHEQVMFRSiaQDyUw6VzbsjIh3I+LIiBgSEUPIXFKaEBHLs/0mSeolqQQYCjyfYq1mZtZIamcQEVEvaRqwGCgC5kVEtaRbgeURUdnC2GpJC4EaoB64zt9gMjNrX+oqy317DsKsa/jwww+pra3l/fff33tny1txcTEDBw7kgAMO+Fi7pBURUZ40pmCT1GZmSWprazn00EMZMmQI2flJa6WIYMuWLdTW1lJSUpL3OC+1YWYdyvvvv88RRxzhcGhDkjjiiCP2+azMAWFmHY7Doe3tz3vqgDAzs0SegzBrwQPL3uCxVa2/BaeibACTxw5ug4q6p6KiIoYPH059fT3HH3888+fP5+CDD27S75RTTuGZZ54pQIVNnXnmmdx+++2UlyfO/wJw1113MXXq1MTfZX8tXbqUAw88kFNOOaXVx3JAmLXgsVUbqdlUR2n/Pon7azbVATS7H2DZq1tZ9urWZoMmn2NA9w6Zgw46iFWrVgHw5S9/mdmzZzN9+vTd++vr6+nZs2cq4dBw7DTcddddXHbZZW0eEIcccogDwqw9lPbvw4NfPXm/x+/tLGRvwQB7QqS7BkSu008/nTVr1rB06VL+7d/+jb59+/Lyyy/z17/+lUMOOYRt27axdOlSbrrpJg4//HBefPFFLrnkEoYPH85Pf/pTduzYwW9/+1uOPfZYHn/8cb7//e+zc+dOjjjiCO6//34++clPcvPNN/PKK6+wYcMGBg8ezMaNG7n77rspKysD4LTTTmPWrFmceOKJu+vasWMHV155JatXr2bYsGHs2LFj975rr72WqqoqduzYwUUXXcQtt9zC3XffzVtvvcVZZ53FkUceyZIlSxL7AcyYMYPKykp69uzJeeedx+23387mzZu55ppreOONN4BM2AwYMIDZs2dTVFTEr3/9a+655x5OP/30/X+zI6JL/IwePTrM2tols5+JS2Y/U+gyOkwd7aGmpqZJW+/evSMi4sMPP4wJEybEz3/+81iyZEkcfPDBsWHDhib9lixZEocddli89dZb8f7778dRRx0VM2fOjIiIu+66K66//vqIiNi6dWt89NFHERHxi1/8IqZPnx4RETfddFOMGjUqtm/fHhER99577+4xa9eujaTPm5/85Cdx5ZVXRkTE6tWro6ioKKqqqiIiYsuWLRERUV9fH2eccUasXr06IiKOPvro2Lx58+5jJPX7xz/+Eccdd9zuOt95552IiLj00kvj6aefjoiI119/PYYNG7a79ttuuy3v95bMjcuJn6uepDazDm/Hjh2UlZVRXl7O4MGDueqqqwAYM2ZMs9/rP+mkk+jfvz+9evXi2GOP5bzzzgNg+PDhvPbaa0Dmnovzzz+f4cOHc9ttt1FdXb17/IQJEzjooIMAuPjii/nd737Hhx9+yLx585gyZUqTf+9Pf/oTl112GQAjRoxgxIgRu/ctXLiQUaNGMXLkSKqrq6mpqWkyvrl+hx12GMXFxVx11VU88sgjuy9HPfHEE0ybNo2ysjImTJhAXV0d27Zt24d3de98icnMOrzcOYhcvXv3bnZMr169dr/u0aPH7u0ePXpQX18PwNe+9jWmT5/OhAkTWLp0KTfffHPisQ8++GDOPfdcHnvsMRYuXMiKFSvyrv3VV1/l9ttvp6qqir59+zJlypTE+xGa69ezZ0+ef/55/vjHP/LQQw/xs5/9jCeffJKPPvqI5557juLi4rxr2Vc+gzCzbuvdd99lwIDMo2bmz5/fYt+rr76ar3/965x00kn07du3yf5x48bxwAMPAPCXv/yFNWvWAFBXV0fv3r057LDDePvtt/n973+/e8yhhx7Ke++912K/bdu28e677/L5z3+eO++8k9WrVwNw3nnncc899+w+VkOA5h6ztRwQZtZt3XzzzVx88cWMHj2aI488ssW+o0ePpk+fPlx55ZWJ+6+99lq2bdvG8ccfz8yZMxk9ejQAJ554IiNHjmTYsGFMnjyZU089dfeYqVOnMn78eM4666xm+7333nt88YtfZMSIEZx22mnccccdANx9990sX76cESNGUFpayuzZswH40pe+xKOPPkpZWRlPP/10q94fL9Zn1oKJc54FaNW3mLpSHe3hpZde4vjjjy90GU289dZbnHnmmbz88sv06NE5/7ZOem9bWqyvc/6WZmbt6L777mPs2LH84Ac/6LThsD88SW1mtheXX345l19+eaHLaHfdJwrNzGyfOCDMzCyRA8LMzBKlGhCSxktaK2m9pBkJ+6+R9KKkVZL+LKk02z5E0o5s+ypJs9Os08y6p4lznt39DTFrKrWAkFQEzAI+B5QClzYEQI4HImJ4RJQBPwbuyNn3SkSUZX+uSatOM7O2VFtbS0VFBUOHDuWYY45h2rRpfPDBBwA8//zzjBs3jk9/+tOMHDmSq6++mu3bt3PvvffSr18/Ro4cydChQzn//PM7xLLlaZ5BjAHWR8SGiNgJLAAqcjtERF3OZm+ga9yUYWbdUkRw4YUXcsEFF7Bu3TrWrVvHjh07+Pa3v83bb7/NxRdfzI9+9CPWrl3LypUrGT9+/O67nidOnMjKlStZt24dM2bM4MILL+Sll14q6O+T5tdcBwBv5mzXAmMbd5J0HTAdOBD4bM6uEkkrgTrguxHRulsCzcxS9uSTT1JcXLz7buuioiLuvPNOjj76aIqKirjiiis4+eQ9NztedNFFicc566yzmDp1KnPnzuXOO+9sl9qTFPw+iIiYBcySNBn4LnAFsAkYHBFbJI0GfivphEZnHEiaCkwFGDzY6+SbWcYtj1dT81Zdk/aG52o02P5BZtG+4Tcv/lh70jM6So/qw01fOqHFf7e6unr3EhsN+vTpw5AhQ1i1ahXXX399XvUDjBo1ijlz5uTdPw1pXmLaCAzK2R6YbWvOAuACgIj4ICK2ZF+vAF4Bjms8ICLmRkR5RJT369evreo2Myu4jrAMUppnEFXAUEklZIJhEjA5t4OkoRGxLrv5BWBdtr0fsDUidkk6BhgKbEixVjPrQvb2l36Dtl7jqrS0lIceeuhjbXV1dfztb3/jnHPOYcWKFVRUVDQz+uNWrlxZ8DWpUjuDiIh6YBqwGHgJWBgR1ZJulTQh222apGpJq8jMQ1yRbR8HrMm2PwRcExFb06rVzKwtnH322Wzfvp377rsPgF27dvHNb36TadOmccMNNzB//nyWLVu2u/8jjzzC22+/3eQ4Tz31FHPnzuUrX/lKu9WeJNU5iIhYBCxq1DYz53XiBbmIeBh4OM3azMzamiQeffRRrrvuOr73ve+xefNmJk6cyI033gjAggULuOGGG/j73/9Ojx49GDduHOPHjwfgwQcf5M9//jPbt2+npKSEhx9+uOBnEAWfpDYz60oGDRpEZWUlAM888wyXXnopL7zwAqNGjeLkk09OfEbDlClTEh9jWmgOCDPrttJ+vsYpp5zC66+/nuq/kSavxWRmZokcEGbW4XSEr3h2NfvznjogzKxDKS4uZsuWLQ6JNhQRbNmyheLi4n0a5zkIM+tQBg4cSG1tLZs3by50KV1KcXExAwcO3KcxDggz61AOOOAASkpKCl2G4UtMZmbWDAeEmZklckCYmVkiB4SZmSVyQJiZWSIHhJmZJXJAmJlZIgeEmZklckCYmVkiB4SZmSVyQJiZWSIHhJmZJUo1ICSNl7RW0npJMxL2XyPpRUmrJP1ZUmnOvu9kx62VdH6adZqZWVOpBYSkImAW8DmgFLg0NwCyHoiI4RFRBvwYuCM7thSYBJwAjAd+nj2emZm1kzTPIMYA6yNiQ0TsBBYAFbkdIqIuZ7M30PCEkApgQUR8EBGvAuuzxzMzs3aS5vMgBgBv5mzXAmMbd5J0HTAdOBD4bM7Y5xqNHZAwdiowFWDw4MFtUrSZmWUUfJI6ImZFxLHAvwLf3cexcyOiPCLK+/Xrl06BZmbdVJoBsREYlLM9MNvWnAXABfs51szM2liaAVEFDJVUIulAMpPOlbkdJA3N2fwCsC77uhKYJKmXpBJgKPB8irWamVkjqc1BRES9pGnAYqAImBcR1ZJuBZZHRCUwTdI5wIfAO8AV2bHVkhYCNUA9cF1E7EqrVjMzayrNSWoiYhGwqFHbzJzX17cw9gfAD9KrzszMWlLwSWozM+uYHBBmZpbIAWFmZokcEGZmlsgBYWZmiRwQZmaWyAFhZmaJHBBmZpbIAWFmZokcEGZmlsgBYWZmiRwQZmaWyAFhZmaJHBBmZpbIAWFmZokcEGZmlsgBYWZmiRwQZmaWKNWAkDRe0lpJ6yXNSNg/XVKNpDWS/ijp6Jx9uyStyv5UplmnmZk1ldozqSUVAbOAc4FaoEpSZUTU5HRbCZRHxHZJ1wI/BiZm9+2IiLK06jMzs5aleQYxBlgfERsiYiewAKjI7RARSyJie3bzOWBgivWYmdk+SDMgBgBv5mzXZtuacxXw+5ztYknLJT0n6YKkAZKmZvss37x5c6sLNjOzPVK7xLQvJF0GlANn5DQfHREbJR0DPCnpxYh4JXdcRMwF5gKUl5dHuxVsZtYNpHkGsREYlLM9MNv2MZLOAW4EJkTEBw3tEbEx+98NwFJgZIq1mplZI2kGRBUwVFKJpAOBScDHvo0kaSQwh0w4/D2nva+kXtnXRwKnArmT22ZmlrLULjFFRL2kacBioAiYFxHVkm4FlkdEJXAbcAjwG0kAb0TEBOB4YI6kj8iE2A8bffvJzMxSluocREQsAhY1apuZ8/qcZsY9AwxPszYzM2uZ76Q2M7NEDggzM0vkgDAzs0QOCDMzS9QhbpQzS8MDy97gsVVNbr3ZJzWb6ijt36eNKmqdmk11TJzz7H6PrygbwOSxg9uwIuvqfAZhXdZjqzZSs6muxT41m+pa7FPavw8VZS2tENM+KsoGtBhUe/s9ajbVtTosrfvZ5zMISX2BQRGxJoV6zNpUaf8+PPjVkwtdRqtNHju4VX/9t+bMw7qvvM4gJC2V1EfSJ4AXgF9IuiPd0szMrJDyvcR0WETUARcC90XEWCDxJjczM+sa8g2InpL6A5cAv0uxHjMz6yDyDYhbyKyptD4iqrJLcK9LrywzMyu0fCepN0XEiIaNiNjgOQgzs64t3zOIe/JsMzOzLqLFMwhJJwOnAP0kTc/Z1YfMEt5mZtZF7e0S04FkntfQEzg0p70OuCitoszMrPBaDIiIeAp4StK9EfF6O9VkZmYdQL6T1L0kzQWG5I6JiM+mUZSZmRVevgHxG2A28EtgV3rlmJlZR5Hvt5jqI+J/R8TzEbGi4WdvgySNl7RW0npJMxL2T5dUI2mNpD9KOjpn3xWS1mV/rtiH38nMzNpAvgHxuKT/Kam/pE80/LQ0QFIRMAv4HFAKXCqptFG3lUB59h6Lh4AfZ8d+ArgJGAuMAW7KLhJoZmbtJN9LTA1/wX8rpy2AY1oYM4bMndcbACQtACqAmt0HiFiS0/854LLs6/OBP0TE1uzYPwDjgf+TZ71mZtZKeQVERJTsx7EHAG/mbNeSOSNozlXA71sY22RRfklTgakAgwf7QShmZm0pr4CQdHlSe0Tc1xZFSLoMKAfO2JdxETEXmAtQXl4ebVGLmZll5HuJ6aSc18XA2WSeC9FSQGwEBuVsD8y2fYykc4AbgTMi4oOcsWc2Grs0z1rNzKwN5HuJ6Wu525IOBxbsZVgVMFRSCZkP/EnA5EbHGQnMAcZHxN9zdi0G/j1nYvo84Dv51GpmZm1jnx85mvVPoMV5iYiolzSNzId9ETAvIqol3Qosj4hK4DYyS3n8RhLAGxExISK2SvoemZABuLVhwtrMzNpHvnMQj5P51hJkPuyPBxbubVxELAIWNWqbmfO62afSRcQ8YF4+9ZmZWdvL9wzi9pzX9cDrEVGbQj1mZtZB5HWjXHbRvpfJrOjaF9iZZlFmZlZ4eQWEpEuA54GLyTyXepkkL/dtZtaF5XuJ6UbgpIZvGknqBzxBZnkMMzPrgvJdi6lHo6+hbtmHsWZm1gnlewbxfyUtZs9aSBNp9O0kMzPrWvb2TOpPAZ+MiG9JuhA4LbvrWeD+tIszM7PC2dsZxF1k72COiEeARwAkDc/u+1KKtZmZWQHtbR7hkxHxYuPGbNuQVCoyM7MOYW8BcXgL+w5qwzrMzKyD2VtALJf0lcaNkq4G9vrIUTMz67z2NgfxDeBRSV9mTyCUAwcC/zXFuszMrMBaDIiIeBs4RdJZwGeyzf8ZEU+mXpmZmRVUvs+DWAIs2WtHMzPrMnw3tJmZJXJAmJlZIgeEmZklckCYmVmiVANC0nhJayWtlzQjYf84SS9Iqm/8fAlJuyStyv5UplmnmZk1le9qrvtMUhEwCzgXqAWqJFVGRE1OtzeAKcANCYfYERFladVnZmYtSy0ggDHA+ojYACBpAVAB7A6IiHgtu++jFOswM7P9kOYlpgHAmznbtdm2fBVLWi7pOUkXJHWQNDXbZ/nmzZtbUaqZmTXWkSepj46IcmAycJekYxt3iIi5EVEeEeX9+vVr/wrNzLqwNANiIzAoZ3tgti0vEbEx+98NwFJgZFsWZ2ZmLUtzDqIKGCqphEwwTCJzNrBXkvoC2yPiA0lHAqcCP06tUuuQHlj2Bo+tav5vippNdQCU9u/T7P7m9nVHNZvqmDjn2VYdo6JsAJPHDm6jiqyjSy0gIqJe0jRgMVAEzIuIakm3AssjolLSScCjQF/gS5JuiYgTgOOBOdnJ6x7ADxt9+8m6gcdWbWzxQ35vH/6l/ftQUbYv015dV1u8Dw2B7IDoPtI8gyAiFgGLGrXNzHldRebSU+NxzwDD06zNOofS/n148KsnF7qMTm/y2MGt/mBv7dmHdT4deZLazMwKyAFhZmaJHBBmZpbIAWFmZokcEGZmlsgBYWZmiRwQZmaWyAFhZmaJHBBmZpbIAWFmZokcEGZmlsgBYWZmiRwQZmaWyAFhZmaJHBBmZpbIAWFmZokcEGZmlsgBYWZmiVINCEnjJa2VtF7SjIT94yS9IKle0kWN9l0haV3254o06zQzs6ZSCwhJRcAs4HNAKXCppNJG3d4ApgAPNBr7CeAmYCwwBrhJUt+0ajUzs6bSPIMYA6yPiA0RsRNYAFTkdoiI1yJiDfBRo7HnA3+IiK0R8Q7wB2B8irWamVkjaQbEAODNnO3abFubjZU0VdJyScs3b96834WamVlTPQtdQGtExFxgLkB5eXkUuBzL8cCyN3hs1cZWHaNmUx2l/fu0UUXWFmo21TFxzrOtOkZF2QAmjx3cRhVZmtI8g9gIDMrZHphtS3usdQCPrdpIzaa6ZvfXbKprcT9Aaf8+VJTle9JpaasoG9BiYOfzv2nNprpW/+Fg7SfNM4gqYKikEjIf7pOAyXmOXQz8e87E9HnAd9q+REtTaf8+PPjVkwtdhrWRyWMHt/ov/9aefVj7Su0MIiLqgWlkPuxfAhZGRLWkWyVNAJB0kqRa4GJgjqTq7NitwPfIhEwVcGu2zczM2kmqcxARsQhY1KhtZs7rKjKXj5LGzgPmpVmfmZk1z3dSm5lZIgeEmZklckCYmVkiB4SZmSVyQJiZWSIHhJmZJXJAmJlZIgeEmZklckCYmVkiB4SZmSVyQJiZWSIHhJmZJXJAmJlZok79RDkz63z8VLrOwwFhTfhxoZaWvT0hsOGJdHt7ch3ggGgHDghrouFxoa35gPfjQi2Jn0rXuTggLJEfF2pmnqQ2M7NEqQaEpPGS1kpaL2lGwv5ekh7M7l8maUi2fYikHZJWZX9mp1mnmZk1ldolJklFwCzgXKAWqJJUGRE1Od2uAt6JiE9JmgT8CJiY3fdKRJSlVZ+ZmbUszTOIMcD6iNgQETuBBUBFoz4VwPzs64eAsyUpxZrMzCxPaQbEAODNnO3abFtin4ioB94FjsjuK5G0UtJTkk5PsU4zM0vQUb/FtAkYHBFbJI0GfivphIioy+0kaSowFWDwYH8n2sysLaUZEBuBQTnbA7NtSX1qJfUEDgO2REQAHwBExApJrwDHActzB0fEXGAuQHl5eaTxS3Q2vsnNzNpKmgFRBQyVVEImCCYBkxv1qQSuAJ4FLgKejIiQ1A/YGhG7JB0DDAU2pFhrl+Gb3Kw7aO1yHV6qIz+pBURE1EuaBiwGioB5EVEt6VZgeURUAv8B/ErSemArmRABGAfcKulD4CPgmojYmlatXY1vcrOurLV/vHipjvylOgcREYuARY3aZua8fh+4OGHcw8DDadZmZp1Ta5fr8FId+fOd1GZmlsgBYWZmiRwQZmaWyAFhZmaJHBBmZpbIAWFmZokcEGZmlsgBYWZmiRwQZmaWqKOu5tpttXaxPS+0Z7Z3rV3LCbrHek4+g+hgGhbb219eaM+sZRVlA1r9R1TNprpWr5rcGfgMogPyYntm6WntWk7QfdZz8hmEmZklckCYmVkiX2JqQ36am1n30R0mun0G0YZaO8EMnmQ26wy6y0S3zyDamCeYzbq+7jLR7YDI8uUhM2tvHf0ylS8xZfnykJm1p85wmSrVMwhJ44GfAkXALyPih4329wLuA0YDW4CJEfFadt93gKuAXcDXI2JxmrWCLw+ZWfvpDJepUgsISUXALOBcoBaoklQZETU53a4C3omIT0maBPwImCipFJgEnAAcBTwh6biI2NXcv7dh8z9b9Wb58pCZdUZtcZmqOWleYhoDrI+IDRGxE1gAVDTqUwHMz75+CDhbkrLtCyLig4h4FVifPV6z/rmzvlXF+vKQmXU2rb1MtezVrS3uT/MS0wDgzZztWmBsc30iol7Su8AR2fbnGo1t8uktaSowNbu5beE1p6xtTcELgS+35gAZRwL/aP1hugS/F3v4vdjD78UeHeG9OLq5HZ36W0wRMReYW+g6cklaHhHlha6jI/B7sYffiz38XuzR0d+LNC8xbQQG5WwPzLYl9pHUEziMzGR1PmPNzCxFaQZEFTBUUomkA8lMOlc26lMJXJF9fRHwZEREtn2SpF6SSoChwPMp1mpmZo2kdokpO6cwDVhM5muu8yKiWtKtwPKIqAT+A/iVpPXAVjIhQrbfQqAGqAeua+kbTB1Mh7rkVWB+L/bwe7GH34s9OvR7ocwf7GZmZh/nO6nNzCyRA8LMzBI5IFIk6ZuSQtKRha6lUCTdJullSWskPSrp8ELX1J4kjZe0VtJ6STMKXU+hSBokaYmkGknVkq4vdE2FJqlI0kpJvyt0Lc1xQKRE0iDgPOCNQtdSYH8APhMRI4C/At8pcD3tJme5mc8BpcCl2WVkuqN64JsRUQr8F+C6bvxeNLgeeKnQRbTEAZGeO4FvA936WwAR8f8iomEdlOfI3NPSXeSz3Ey3EBGbIuKF7Ov3yHwwdtu1bSQNBL4A/LLQtbTEAZECSRXAxohYXehaOpj/Afy+0EW0o6TlZrrth2IDSUOAkcCyApdSSHeR+QPyowLX0aJOvdRGIUl6AviXhF03Av+LzOWlbqGl9yIiHsv2uZHMZYb727M261gkHQI8DHwjIlr3AJZOStIXgb9HxApJZxa4nBY5IPZTRJyT1C5pOFACrM4sTMtA4AVJYyLib+1YYrtp7r1oIGkK8EXg7OheN954yZgckg4gEw73R8Qjha6ngE4FJkj6PFAM9JH064i4rMB1NeEb5VIm6TWgPCIKvWJjQWQfGnUHcEZEbC50Pe0pu77YX4GzyQRDFTA5IqoLWlgBZJfxnw9sjYhvFLicDiN7BnFDRHyxwKUk8hyEpe1nwKHAHyStkjS70AW1l+zkfMNyMy8BC7tjOGSdCvx34LPZ/x+syv4FbR2YzyDMzCyRzyDMzCyRA8LMzBI5IMzMLJEDwszMEjkgzMwskQPCzMwSOSDMzCyRA8IsRZL+VdKvJC2W9LqkGwpdk1m+HBBm6RpOZsXOzwPnAJcXthyz/HmxPrN0DQcujIhdknYBWwtdkFm+fAZhlpLs6qX9IuKVbNMI4MUClmS2TxwQZun5NB9/pGQZ4IdIWafhgDBLz3BgTc52GQ4I60S8mqtZO5G0HhgeETsKXYtZPnwGYdYOJB0O7HQ4WGfiMwgzM0vkMwgzM0vkgDAzs0QOCDMzS+SAMDOzRA4IMzNL5IAwM7NEDggzM0v0/wFdt1x7rO+lIgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = hist.plotgrid(output['eta'], overlay=\"dataset\", stack=False, density=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot $p_{T}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    " # your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## Plot JER and response correction for $\\eta$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEGCAYAAADIRPqpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAj4ElEQVR4nO3dfZxU5X338c83+ICJIqjUm/AQsG5tUBPFrYhprQaraJOs9VZR0ohWJWkgsWlMqrV3MTHe0aYtiW1ipUqB3FGwxATSYAkVzZOAohARjGHFp6WrEDBuEnwI9nf/ca6Bwzi7DLAzZ9j9vl+vee2Z37nOnN8chv3tdc4111FEYGZmVqS3FZ2AmZmZi5GZmRXOxcjMzArnYmRmZoVzMTIzs8LtV3QCjeKII46I4cOHF52Gmdk+5dFHH/15RAzc29dxMUqGDx/OihUrik7DzGyfIum57ngdn6YzM7PCuRiZmVnhXIzMzKxwvmZkZrYXfvOb39DW1sZrr71WdCo11bdvX4YMGcL+++9fk9d3MTIz2wttbW0ccsghDB8+HElFp1MTEcHmzZtpa2tjxIgRNdmHT9OZme2F1157jcMPP7zHFiIASRx++OE17f25GJmZ7aWeXIhKav0eXYzMzKxwLkZmDeSu5c8z/valOz3uWv580WnZPqCtrY2Wlhaampo46qijmDJlCq+//joADz/8MKeddhrHHHMMJ554IldeeSVbt25l5syZDBw4kBNPPJGmpibOPvtsHnrooULydzEyayDzV21gbXvHTo/5qzYUnZY1uIjg/PPP57zzzmPdunWsW7eOV199lc9+9rO89NJLXHjhhdxyyy089dRTrFy5knHjxvHLX/4SgPHjx7Ny5UrWrVvHtddey/nnn8+TTz5Z9/fg0XRmDWbkoH7M/egYAMbfvrTgbGxfsGTJEvr27cvll18OQJ8+fZg2bRrvete76NOnDxMnTmTMmDHb219wwQUVX+eMM85g0qRJTJ8+nWnTptUl9xIXIzOzbvK576xh7X93dOtrjnxnP6Z+8Ngu26xZs4aTTjppp1i/fv0YPnw4q1at4uqrr656f6NGjeL222/fo1z3hk/TmZnZdhFRyH7dMzJrcGvbOzj+hkVAdgoPoOWEwUwYPazItKyCXfVgamXkyJHMmzdvp1hHRwcvvvgiZ555Jo8++igtLS1VvdbKlSt597vfXYs0u+SekVkDazlhMCMH9dv+ADyowd5i7NixbN26ldmzZwPw5ptv8ulPf5opU6ZwzTXXMGvWLJYvX769/b333stLL730ltf5/ve/z/Tp07nqqqvqlnuJe0ZmDWzC6GFv6QF5UIOVk8S3vvUtJk+ezI033simTZsYP348119/PQBz5szhmmuuYePGjbztbW/jtNNOY9y4cQDMnTuXH/3oR2zdupURI0bwzW9+s5CeUc2KkaQZwAeAjRFxXNm6TwN/DwyMiJ8r+2rvV4Bzga3AZRHxWGo7EfibtOkXImJWip8EzAQOAhYCV0dESDoMmAsMB54FLoqIl2v1Ps3MGsHQoUNZsGABAA899BCXXHIJjz32GKNGjWLMmDH88Ic/fMs2l112GZdddlmdM62slqfpZgLjyoOShgJnAflv8p0DNKXHJOC21PYwYCowGjgZmCppQNrmNuCq3HalfV0L3B8RTcD96bmZWa9x6qmn8txzzzFq1KiiU6lazYpRRPwA2FJh1TTgs0B+yEYLMDsyy4D+kgYBZwOLI2JL6t0sBsaldf0iYllkQz9mA+flXmtWWp6Vi5uZWYOq6wAGSS3Ahoj4SdmqwcALuedtKdZVvK1CHODIiGhPyy8CR3ZP9mZmlRU1HLqeav0e61aMJL0d+Gvgb+u1z9Rr6vQISpokaYWkFZs2bapXWmbWg/Tt25fNmzf36IJUup9R3759a7aPeo6m+21gBPCTNBX5EOAxSScDG4ChubZDUmwDcHpZ/MEUH1KhPcBLkgZFRHs6nbexs4QiYjowHaC5ubnnfpLMrGaGDBlCW1sbPf0P2tKdXmulbsUoIlYDv1V6LulZoDmNplsATJE0h2ywwiupmCwC/m9u0MJZwHURsUVSh6RTgOXApcA/pTYLgInAzenn/Dq8PTPrpfbff/+a3f20N6nZaTpJdwNLgWMktUm6oovmC4H1QCvwr8DHASJiC3Aj8Eh6fD7FSG3uSNs8DdyX4jcDfyRpHXBmem5mZg2sZj2jiLhkF+uH55YDmNxJuxnAjArxFcBxFeKbgbG7ma6ZmRXI0wGZmVnhXIzMzKxwLkZmZlY4FyMzMyuci5GZmRXOxcjMzArnYmRmZoVzMTIzs8K5GJmZWeFcjMzMrHAuRmZmVjgXIzMzK5yLkZmZFc7FyMzMCudiZGZmhXMxMjOzwrkYmZlZ4VyMzMyscC5GZmZWuJoVI0kzJG2U9EQu9iVJP5X0uKRvSeqfW3edpFZJT0k6Oxcfl2Ktkq7NxUdIWp7icyUdkOIHpuetaf3wWr1HMzPrHrXsGc0ExpXFFgPHRcR7gJ8B1wFIGglcDBybtvmapD6S+gBfBc4BRgKXpLYAtwDTIuJo4GXgihS/Ang5xaeldmYN6a7lzzP+9qUcf8Mijr9hEWvbO4pOyawQNStGEfEDYEtZ7HsRsS09XQYMScstwJyIeD0ingFagZPTozUi1kfEG8AcoEWSgPcD89L2s4Dzcq81Ky3PA8am9mYNZ/6qDaxt72DkoH7bHy0nDC46LbO626/Aff8ZMDctDyYrTiVtKQbwQll8NHA48ItcYcu3H1zaJiK2SXoltf95d78Bs+4wclA/5n50TNFpmBWqkAEMkq4HtgHfKGL/uTwmSVohacWmTZuKTMXMrFerezGSdBnwAeDDEREpvAEYmms2JMU6i28G+kvaryy+02ul9Yem9m8REdMjojkimgcOHLiX78zMzPZUXYuRpHHAZ4EPRcTW3KoFwMVpJNwIoAl4GHgEaEoj5w4gG+SwIBWxB4AL0vYTgfm515qYli8AluSKnpmZNaCaXTOSdDdwOnCEpDZgKtnouQOBxWlMwbKI+FhErJF0D7CW7PTd5Ih4M73OFGAR0AeYERFr0i7+Cpgj6QvASuDOFL8T+LqkVrIBFBfX6j2amVn3qFkxiohLKoTvrBArtb8JuKlCfCGwsEJ8Pdlou/L4a8CFu5WsmZkVyjMwmJlZ4VyMzMyscC5GZmZWuCK/9Gpme2htewfH37AIYPusDRNGDys4K7M952Jkto8pny6oNJ+di5Hty1yMzPYxE0YP26nwjL99aYHZmHUPXzMyM7PCuRiZmVnhXIzMzKxwLkZmZlY4FyMzMyuci5GZmRXOxcjMzArnYmRmZoVzMTIzs8K5GJmZWeFcjMzMrHAuRmZmVjgXIzMzK5yLkZmZFa5mxUjSDEkbJT2Rix0mabGkdenngBSXpFsltUp6XNKo3DYTU/t1kibm4idJWp22uVWSutqHmZk1rlr2jGYC48pi1wL3R0QTcH96DnAO0JQek4DbICsswFRgNHAyMDVXXG4DrsptN24X+zAzswZVs2IUET8AtpSFW4BZaXkWcF4uPjsyy4D+kgYBZwOLI2JLRLwMLAbGpXX9ImJZRAQwu+y1Ku3DzMwaVL2vGR0ZEe1p+UXgyLQ8GHgh164txbqKt1WId7WPt5A0SdIKSSs2bdq0B2/HzMy6Q2EDGFKPJorcR0RMj4jmiGgeOHBgLVMxM7Mu1LsYvZROsZF+bkzxDcDQXLshKdZVfEiFeFf7MDOzBlXvYrQAKI2ImwjMz8UvTaPqTgFeSafaFgFnSRqQBi6cBSxK6zoknZJG0V1a9lqV9mFmZg1qv1q9sKS7gdOBIyS1kY2Kuxm4R9IVwHPARan5QuBcoBXYClwOEBFbJN0IPJLafT4iSoMiPk42Yu8g4L70oIt9mJlZg6pZMYqISzpZNbZC2wAmd/I6M4AZFeIrgOMqxDdX2oeZmTUuz8BgZmaFczEyM7PCuRiZmVnhXIzMzKxwLkZmZlY4FyMzMyuci5GZmRXOxcjMzArnYmRmZoVzMTIzs8Lt9nRAacLSoRHxeA3yMevR7lr+PPNXbdj+fG17ByMH9SswI7PGUFUxkvQg8KHU/lFgo6QfR8Rf1jA3sx5n/qoNrG3v2P585KB+tJwwuIstzHqHantGh0ZEh6QryW4PPlWSe0Zme2DkoH7M/eiYotMwayjVXjPaL92o7iLgP2qYj5mZ9ULVFqPPkd3orjUiHpF0FLCudmmZmVlvUu1puvaIeE/pSUSsl/SPNcrJzMx6mWp7Rv9UZczMzGy3ddkzkjQGOBUYKCk/cq4f0KeWiZmZWe+xq9N0BwAHp3aH5OIdwAW1SsrMzHqXLotRRHwf+L6kmRHxXHftVNKngCuBAFYDlwODgDnA4WTfZfpIRLwh6UBgNnASsBkYHxHPpte5DrgCeBP4ZEQsSvFxwFfIem93RMTN3ZW7mZl1v2qvGR0oabqk70laUnrsyQ4lDQY+CTRHxHFkBeNi4BZgWkQcDbxMVmRIP19O8WmpHZJGpu2OBcYBX5PUR1If4KvAOcBI4JLU1szMGlS1o+n+HfgX4A6yXkh37PcgSb8B3g60A+8HJqT1s4AbgNuAlrQMMA/4Z0lK8TkR8TrwjKRW4OTUrjUi1gNImpParu2GvM3MrAaqLUbbIuK27thhRGyQ9PfA88CrwPfITsv9IiK2pWZtQGmOlMHAC2nbbZJeITuVNxhYlnvp/DYvlMVHV8pF0iRgEsCwYcP27o2Zmdkeq/Y03XckfVzSIEmHlR57ssM00WoLMAJ4J/AOstNsdRcR0yOiOSKaBw4cWEQKZmZG9T2jiennZ3KxAI7ag32eCTwTEZsAJN0LvA/oL2m/1DsaApSmNt4ADAXaJO0HHEo2kKEUL8lv01ncrEda297B8TcsAnZMvjphtHv7tu+oqhhFxIhu3OfzwCmS3k52mm4ssAJ4gGy4+Byy4jc/tV+Qni9N65dEREhaANyVZoJ4J9AEPAwIaJI0gqwIXcyOa1FmPU75rN+lWcFdjGxfUu0tJC6tFI+I2bu7w4hYLmke8BiwDVgJTAe+C8yR9IUUuzNtcifw9TRAYQtZcSEi1ki6h2xgwjZgckS8mfKdQjaXXh9gRkSs2d08zfYVE0YP26nwjL99aYHZmO2Zak/T/V5uuS9Zb+Yxsu//7LaImApMLQuvZ8douHzb14ALO3mdm4CbKsQXAgv3JDczM6u/ak/TfSL/XFJ/stNpZmZme63a0XTlfk02Gs7MzGyvVXvN6Dtko+cguw7zbuCeWiVlZma9S7XXjP4+t7wNeC4i2mqQj5mZ9UJVnaZLE6b+lGzm7gHAG7VMyszMepeqipGki8i+w3MhcBGwXJJvIWFmZt2i2tN01wO/FxEbASQNBP6LbOJSMzOzvVLtaLq3lQpRsnk3tjUzM+tStT2j/5S0CLg7PR+Pv1RqZmbdpMtiJOlo4MiI+Iyk84HfT6uWAt+odXJmZtY77Kpn9GXgOoCIuBe4F0DS8WndB2uYm5mZ9RK7uu5zZESsLg+m2PCaZGRmZr3OropR/y7WHdSNeZiZWS+2q2K0QtJV5UFJV5LdKtzMzGyv7eqa0V8A35L0YXYUn2bgAOBPapiXmZn1Il0Wo4h4CThV0hnAcSn83YhYUvPMzMys16j2fkYPkN0W3MzMrNt5FgUzMyuci5GZmRWukGIkqb+keZJ+KulJSWMkHSZpsaR16eeA1FaSbpXUKulxSaNyrzMxtV8naWIufpKk1WmbWyWpiPdpZmbVKapn9BXgPyPid4H3Ak8C1wL3R0QTcH96DnAO0JQek4DbACQdBkwFRgMnA1NLBSy1uSq33bg6vCczM9tDdS9Gkg4FTgPuBIiINyLiF0ALMCs1mwWcl5ZbgNmRWQb0lzQIOBtYHBFbIuJlYDEwLq3rFxHLIiKA2bnXMjOzBlTtrN3daQSwCfg3Se8l+/7S1WRTD7WnNi8CR6blwcALue3bUqyreFuF+FtImkTW22LYsGF7/o7MOnHX8ueZv2rD9udr2zsYOahfgRmZNaYiTtPtB4wCbouIE4Ffs+OUHACpRxO1TiQipkdEc0Q0Dxw4sNa7s15o/qoNrG3v2P585KB+tJxQ8W8js16tiJ5RG9AWEcvT83lkxeglSYMioj2daivdzG8DMDS3/ZAU2wCcXhZ/MMWHVGhvVoiRg/ox96Njik7DrKHVvWcUES8CL0g6JoXGAmuBBUBpRNxEYH5aXgBcmkbVnQK8kk7nLQLOkjQgDVw4C1iU1nVIOiWNors091pmZtaAiugZAXwC+IakA4D1wOVkhfEeSVcAzwEXpbYLgXOBVmBraktEbJF0I/BIavf5iNiSlj8OzCSbWfy+9DAzswZVSDGKiFVkE66WG1uhbQCTO3mdGcCMCvEV7JhLz8zMGpxnYDAzs8K5GJmZWeFcjMzMrHAuRmZmVjgXIzMzK1xRQ7vNrIbWtndw/A2LgB2zPkwY7SmvrHG5GJn1MOXTDZWmI3IxskbmYmTWw0wYPWynwjP+9qUFZmNWHV8zMjOzwrkYmZlZ4VyMzMyscC5GZmZWOBcjMzMrnIuRmZkVzsXIzMwK52JkZmaFczEyM7PCuRiZmVnhXIzMzKxwhRUjSX0krZT0H+n5CEnLJbVKmivpgBQ/MD1vTeuH517juhR/StLZufi4FGuVdG3d35yZme2WIidKvRp4EuiXnt8CTIuIOZL+BbgCuC39fDkijpZ0cWo3XtJI4GLgWOCdwH9J+p30Wl8F/ghoAx6RtCAi1tbrjVnvdNfy55m/asNOsbXtHYwc1K+TLcyspJCekaQhwB8Dd6TnAt4PzEtNZgHnpeWW9Jy0fmxq3wLMiYjXI+IZoBU4OT1aI2J9RLwBzEltzWpq/qoN22/XsLa9Y3shKr+lg5m9VVE9oy8DnwUOSc8PB34REdvS8zag9D94MPACQERsk/RKaj8YWJZ7zfw2L5TFR1dKQtIkYBLAsGG+14vtvZGD+jH3o2OKTsNsn1P3npGkDwAbI+LReu+7XERMj4jmiGgeOHBg0emYmfVaRfSM3gd8SNK5QF+ya0ZfAfpL2i/1joYApZPvG4ChQJuk/YBDgc25eEl+m87iZmbWgOreM4qI6yJiSEQMJxuAsCQiPgw8AFyQmk0E5qflBek5af2SiIgUvziNthsBNAEPA48ATWl03gFpHwvq8NbMzGwPNdJtx/8KmCPpC8BK4M4UvxP4uqRWYAtZcSEi1ki6B1gLbAMmR8SbAJKmAIuAPsCMiFhT13diZma7pdBiFBEPAg+m5fVkI+HK27wGXNjJ9jcBN1WILwQWdmOqZmZWQ43UMzKzGlnb3sH425duf95ywmAmjPYIUmscLkZmPVz595xK34VyMbJG4mJk1sNNGD1sp8KT7yGZNQpPlGpmZoVzMTIzs8K5GJmZWeFcjMzMrHAuRmZmVjgXIzMzK5yLkZmZFc7fMzLbQ6U7u5a+RAr4rq5me8jFyGwPlQpRvgD5rq5me8bFyGwv+M6uZt3D14zMzKxw7hmZ9UKexdsajYuRWS/jWbytEbkYmfUynsXbGpGvGZmZWeFcjMzMrHB1L0aShkp6QNJaSWskXZ3ih0laLGld+jkgxSXpVkmtkh6XNCr3WhNT+3WSJubiJ0lanba5VZLq/T7NzKx6RfSMtgGfjoiRwCnAZEkjgWuB+yOiCbg/PQc4B2hKj0nAbZAVL2AqMBo4GZhaKmCpzVW57cbV4X2ZmdkeqvsAhohoB9rT8i8lPQkMBlqA01OzWcCDwF+l+OyICGCZpP6SBqW2iyNiC4CkxcA4SQ8C/SJiWYrPBs4D7qvD27MerDT9T0n57AtmtucKHU0naThwIrAcODIVKoAXgSPT8mDghdxmbSnWVbytQrzS/ieR9bYYNszDWq1rleah6ynT/5R/7wj83SOrr8KKkaSDgW8CfxERHfnLOhERkqLWOUTEdGA6QHNzc833Z/u+njj9T76g5ost+LtHVj+FFCNJ+5MVom9ExL0p/JKkQRHRnk7DbUzxDcDQ3OZDUmwDO07rleIPpviQCu3NrILy7x2Bv3tk9VfEaDoBdwJPRsQ/5lYtAEoj4iYC83PxS9OoulOAV9LpvEXAWZIGpIELZwGL0roOSaekfV2aey0zM2tARfSM3gd8BFgtaVWK/TVwM3CPpCuA54CL0rqFwLlAK7AVuBwgIrZIuhF4JLX7fGkwA/BxYCZwENnABQ9eMDNrYEWMpvsR0Nn3fsZWaB/A5E5eawYwo0J8BXDcXqRpZmZ15LnpzKwiz+xt9eRiZNaJ3nxb8dIIO4+us3pxMTLrRG++rbhn9rZ6czEy60JP/F6RWSNyMTKzqvgaktWSi5GZ7ZKvIVmtuRiZJZ4ItXO+hmS15mJklpQPWOhJE6HWgk/bWXdyMTLL8YCF6vi0nXU3FyMz220+bWfdzcXIei1fI+pePm1ne6OI246bNYSefLO8ems5YfBOhXxte8dOhd5sV9wzsl7N14i6R6XTdr57rO0OFyPrFcpPyYFPy9VSpR7m8me2sPyZLTv9O7g4WYmLkfUKleaZ82m52ql099jyPwhcnCxP2e2CrLm5OVasWFF0GtZNOhuc4FNyjaN8VvRfvrYNgNEjDtvexsWp8Ul6NCKa9/Z13DOyHqHSX92w4xebe0GNp7z35J5T7+aeUeKe0b5lV8UH/ItrX1fNv3El/nevr+7qGbkYJS5Gjc3FxyoNQilXbcHK8+dm77gY7YKkccBXgD7AHRFxc1ftXYyKs6e/ZPxLxMpV81nK25PitSd68mfVxagLkvoAPwP+CGgDHgEuiYi1nW3jYlSd3f3PXg2ffrGi1OLzXNLVwIye5J6PneoBDF04GWiNiPUAkuYALUCnxWj9pl97fq0q1OIvydEjDnOhsUJUGoLe3WpZ8HqSntozugAYFxFXpucfAUZHxJSydpOASenpccATdU10zxwB/LzoJKrgPLvPvpAjOM/utq/keUxEHLK3L9JTe0ZViYjpwHQASSu6o6tZa86ze+0Lee4LOYLz7G77Up7d8To9daLUDcDQ3PMhKWZmZg2opxajR4AmSSMkHQBcDCwoOCczM+tEjzxNFxHbJE0BFpEN7Z4REWt2sdn02mfWLZxn99oX8twXcgTn2d16VZ49cgCDmZntW3rqaTozM9uHuBiZmVnhekUxkjRO0lOSWiVdW2H9gZLmpvXLJQ3PrbsuxZ+SdHaBOf6lpLWSHpd0v6R35da9KWlVetR0oEYVeV4maVMunytz6yZKWpceEwvOc1oux59J+kVuXV2Op6QZkjZKqvj9NmVuTe/hcUmjcuvqeSx3leeHU36rJT0k6b25dc+m+KruGgK8F3meLumV3L/t3+bWdfl5qXOen8nl+ET6PB6W1tXleEoaKumB9DtnjaSrK7Tp3s9nRPToB9kAhqeBo4ADgJ8AI8vafBz4l7R8MTA3LY9M7Q8ERqTX6VNQjmcAb0/Lf17KMT3/VQMdy8uAf66w7WHA+vRzQFoeUFSeZe0/QTbIpd7H8zRgFPBEJ+vPBe4DBJwCLK/3sawyz1NL+wfOKeWZnj8LHNEgx/N04D/29vNS6zzL2n4QWFLv4wkMAkal5UPIplcr/7/erZ/P3tAz2j41UES8AZSmBsprAWal5XnAWElK8TkR8XpEPAO0ptere44R8UBEbE1Pl5F9d6reqjmWnTkbWBwRWyLiZWAxMK5B8rwEuLtGuXQqIn4AbOmiSQswOzLLgP6SBlHfY7nLPCPioZQHFPfZrOZ4dmZvPte7bTfzLOqz2R4Rj6XlXwJPAuU3BOvWz2dvKEaDgRdyz9t460Hd3iYitgGvAIdXuW29csy7guwvkpK+klZIWibpvBrkV1Jtnv87ddvnSSp9+bhex3K39pVOd44AluTC9Tqeu9LZ+6jnsdxd5Z/NAL4n6VFl028VbYykn0i6T9KxKdaQx1PS28l+iX8zF6778VR22eJEYHnZqm79fPbI7xn1ZJL+FGgG/jAXfldEbJB0FLBE0uqIeLqYDPkOcHdEvC7po2Q9zvcXlEs1LgbmRcSbuVgjHc99hqQzyIrR7+fCv5+O5W8BiyX9NPUMivAY2b/trySdC3wbaCool2p8EPhxROR7UXU9npIOJiuGfxERHbXaD/SOnlE1UwNtbyNpP+BQYHOV29YrRySdCVwPfCgiXi/FI2JD+rkeeJDsr5ha2GWeEbE5l9sdwEnVblvPPHMupuw0SB2P56509j4abrorSe8h+/duiYjNpXjuWG4EvkVtTnNXJSI6IuJXaXkhsL+kI2jA45l09dms+fGUtD9ZIfpGRNxboUn3fj5rfSGs6AdZ72892amY0sXJY8vaTGbnAQz3pOVj2XkAw3pqM4ChmhxPJLvI2lQWHwAcmJaPANZRo4uvVeY5KLf8J8Cy2HFR85mU74C0fFhReaZ2v0t2QVhFHM+0j+F0fsH9j9n5AvHD9T6WVeY5jOx66qll8XcAh+SWHyKbTb+oPP9X6d+a7Jf48+nYVvV5qVeeaf2hZNeV3lHE8UzHZTbw5S7adOvns2YHu5EeZKM+fkb2y/z6FPs8WQ8DoC/w7+k/1MPAUbltr0/bPQWcU2CO/wW8BKxKjwUpfiqwOv0HWg1cUfCx/CKwJuXzAPC7uW3/LB3jVuDyIvNMz28Abi7brm7Hk+yv3nbgN2Tn1a8APgZ8LK0X8NX0HlYDzQUdy13leQfwcu6zuSLFj0rH8SfpM3F9wXlOyX02l5ErnpU+L0XlmdpcRjZ4Kr9d3Y4n2anWAB7P/bueW8vPp6cDMjOzwvWGa0ZmZtbgXIzMzKxwLkZmZlY4FyMzMyuci5GZmRXOxch6nNys209I+o6k/ik+XNKruRmRV0m6NK07WNLtkp5OU608KGl0WjdE0vw0A/HTkr6i7Hb2ne3/PklvmZ9N2czrP02zLv9E0j+mLxaWz8a8StKtKT5T0gZJB6bnR0h6truPmVnRXIysJ3o1Ik6IiOPIvjg4Obfu6bSu9Jid4nektk0RcRJwOXBEmjD3XuDbEdEE/A5wMHBTpR1LOgg4PCLayuIfA84CTomI44HfAzYCB+WanZHL65O5+Jtk39uoSpra3/+3bZ/iD6z1dEvZxSSNkn4bGA38TUT8D0BEPBMR3yWbV++1iPi3FH8T+BTwZ2kiy3Knk00hVO564M8j4hfpdd6IiJujuvm+vgx8Kk1V1dl7GK7sfjyzgSeAocrui/NImrT2c6ndOyR9N/XMnpA0PsWflfR3qXf2sKSjc6+7RDvuozUsxWcqu5fNQ5LWS7ogxQdJ+kGuZ/oHKX6WpKWSHpP072nOM7PtXIysx5LUBxgL5G+Q99tlp+n+gGzap1Wx82SpJccCj+YDqYA8Dxxdof05wH+W5dEPODiy25B05YFcXp/KxZ8HfgR8ZBfbNwFfi4hjgWPS85OBE4CTJJ1GNgv0f0fEe1PPMZ/rK6nX9s9kBRDgn4BZEfEe4BvArbn2g8i+qf8B4OYUmwAsiogTgPcCq9L8b38DnBkRo4AVwF/u4r1YL+NZu60nOkjSKrIe0ZNk91MpeTr9otxO0oe6cd/vA67pqoGyOwbfAvQHJkTEQ2nVGRHx8042+yIwH/huFy/9XGT3lYHslOBZwMr0/GCy4vRD4B8k3UJ2o7kf5ra/O/dzWloeA5yflr8O/F2u/bdTT3KtpCNT7BFgRroW9u2IWCXpD8luVPnj7KwnB5D1WM22c8/IeqJXU8F5F9n8WZO7bs4a4L2pJ1VuLTtmHge293RKk4Pm40cBL0R2g7btUk/qV5JGpOelnsMTZL+Ydyki1pHND3ZRF81+nU8H+GLuGtTREXFnRPyM7C6jq4EvKHfrbbK5yCotd+b13LJSnj8gu5PpBmBmGiAisputlXIZGRFXVPH61ou4GFmPFdmdcT8JfLqr6y2R3atoBfC5NGChdK3kj4H7gbfnRt31Af4BmBk77rxb8pZTdDlfBG7LjewT2QS9u+MmdtHryllEdl3r4LS/wZJ+S9I7ga0R8f+AL5EVppLxuZ+lnstDZDPZA3yYrGfVKWU3K3wpIv6VbFDIKLJJSd+Xuw71Dkm/U+X7sF7Cp+msR4uIlZIeJ7t98w9J14xyTWZExK3AlWRFplXSq8DPgc9EREj6E+Brkv4P2R9wC4G/rrC7ccAnOknlNrJp/5dLeh34FfBjdpxGg+yaUem61eMRcWnZe1kj6TF2LiCdve/vSXo3sDTV118Bf0p2netLkv6HbNboP89tNiAdq9fJjhfp/fybpM8Am8hGGXbldOAzkn6T9nlpRGySdBlwd2mIOtk1pJ/t6n1Y7+FZu826Qfol++OIaC46lz2RvrvU3MU1K7Oacs/IrBtEdnfbfbIQmTUC94zMzKxwHsBgZmaFczEyM7PCuRiZmVnhXIzMzKxwLkZmZla4/w+l5q1mmh6LYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEICAYAAACuxNj9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuaElEQVR4nO3deXxU9bnH8c+TPWSHBEiAJIDsECBEpaJ1AfcoWtuqvVastt7b2lZre62tWm1rbW3vba1tr71uldYu3rpUQdSKW7VWKyQQIGxhSwIJBMgesj/3jznBmJJkkszMmUme9+s1L+ac2b55afLM+a2iqhhjjDF9CXM7gDHGmOBnxcIYY0y/rFgYY4zplxULY4wx/bJiYYwxpl8Rbgfwh9TUVM3OznY7hjHGhJT169cfVtW0Ez02LItFdnY269atczuGMcaEFBHZ19tj1gxljDGmX1YsjDHG9MuKhTHGmH5ZsTDGGNMvKxbGGGP6ZcXCGGNMv6xYGGOM6dewnGdhjFv+vK6Mww2txMdEkBAdQXx0BPExnn8TYj48jo4IdzuqMQNixcIYH6moPcZ/Pl3k1XOjwsOOF5GuApIQHUFibCT/ceZUZoxP8HNaYwbGioUxPlKwrwaAp25czJS0eBpa2mlobqe+pY2G5nbPcUs79c2eW0O38/XN7VTWNfPWjiriosO597J57v4wxvRgxcIYHyksrSY6IoyFmSlERYSRlhA94Pe45tH3jxcdY4KJ3zq4ReRxETkkIpu7nRstIq+KyE7n3xTnvIjIgyJSIiJFIpLb7TUrnOfvFJEV/sprzFAVltUwb0ISURGD/7XKzUxmW2UdjS3tPkxmzND5czTUE8AFPc7dDrymqtOA15xjgAuBac7tRuAh8BQX4G7gVOAU4O6uAmNMMGlt72TT/loWZiYP6X0WZqXQqbCxvMYnuYzxFb8VC1X9G3C0x+nlwErn/krgsm7nf6se7wHJIpIOnA+8qqpHVbUaeJV/LUDGuK64oo7W9k5yM4f2XSZ3kuf1haU1PkhljO8Eep7FOFWtcO5XAuOc+xOAsm7PK3fO9Xb+X4jIjSKyTkTWVVVV+Ta1Mf0o2FcNwMIhFoukUZFMTYs7/n7GBAvXJuWpqgLqw/d7WFXzVDUvLe2Ee3cY4zeFZTWkJ8UwPilmyO+Vm5lCQWk1nl8RY4JDoIvFQad5CeffQ875/cCkbs+b6Jzr7bwxQaWwtHrITVBdcrNSqG5qY8/hRp+8nzG+EOhi8QLQNaJpBfB8t/PXOqOiFgO1TnPVK8B5IpLidGyf55wzJmgcqm+mvPrYkDu3u3QVnQLrtzBBxJ9DZ/8I/AOYISLlInID8CPgXBHZCSxzjgHWALuBEuAR4EsAqnoU+D7wgXP7nnPOmKDR1Rntq2IxbWw8CdERFJRav4UJHn6blKeqV/fy0NITPFeBm3p5n8eBx30YzRifKiitJjJcmJOR5JP3CwsTFmQmWye3CSq26qwxQ1RYWsPsjCRiIn23OGBuZgo7DtbTYJPzTJCwYmHMELR3dFJUXkOuj5qguuR2Tc4rq/Hp+xozWFYsjBmCbZX1NLd1Dnl+RU8LJiUDWFOUCRpWLIwZgkKnE3qh88fdV5JiI5k2Nt46uU3QsGJhzBAUlNaQlhDNxJRYn793bmYKhWU1NjnPBAUrFsYMQWFpNQsnJSMiPn/v3Kxkapra2G2T80wQsGJhzCAdbWxl75EmcrP8sxDy8cl51m9hgoAVC2MGyV/9FV2mpsWTGGOT80xwsGJhzCAVltYQHibMm+ibyXg9eSbnpdjOeSYoWLEwZpAKSquZlZ7AqCj/7U68KDOFHYfqqWtu89tnGOMNKxbGDEJHp7KxrIaFk/y7cWNuVjJqk/NMELBiYcwg7DxUT2NrB7lZyX79nAWTkhHBmqKM66xYGDMIXX+8/X1lkRATyfSxCdbJbVxnxcKYQSgsrWZ0XBRZY0b5/bNys5IpLK2ms9Mm5xn3WLEwZhAK/DgZr6eFmSnUNbez+3CD3z/LmN5YsTBmgGqb2thV1eizzY768+HkvJqAfJ4xJ2LFwpgB2lBeA+DzlWZ7MyU1jqTYSOu3MK6yYmHMABXsq0YE5vtp5nZPYWHCwsxkKxbGVVYsjBmgwrIaZoxLID7af5PxevLsnNdA7TGbnGfcYcXCmAHo7FTPSrMBaoLq0tVvscEm5xmXWLEwZgB2H26gvrk9YJ3bXeZPSiJMbAVa4x4rFsYMQEFpDYDP99zuT0JMJNPH2eQ84x4rFsYMQGFpNYkxEUxJjQ/4Z+dmpbChrMYm5xlXWLEwZgAKS2tYkJlCWJj/J+P1lJuZQn1zOyVVNjnPBJ4VC2O81NDSzvaD9QFvgurS9bnWb2HcYMXCGC9tLKtBNXCT8XqanBpHyiibnGfcYcXCGC91baO6YGKyK58vIizMTDneyW5MIFmxMMZLBaU1nDQ2nqRRka5lyM1MpuRQA7VNNjnPBJYVC2O8oOpMxgvQEh+96ZqcV1hmTVEmsKxYGOOFfUeaqG5qIzfLnf6KLvMnJXsm51lTlAkwV4qFiHxNRLaIyGYR+aOIxIjIZBF5X0RKROQpEYlynhvtHJc4j2e7kdmMbF2dyoGeud1TXHQEM8cn2ogoE3ABLxYiMgH4KpCnqnOBcOAq4H7gZ6p6ElAN3OC85Aag2jn/M+d5xgRUYWkN8dERTBub4HYUcrOS2VBWQ4dNzjMB5FYzVAQQKyIRwCigAjgHeNp5fCVwmXN/uXOM8/hSCcT2ZMZ0U1BazfxJSYS7MBmvp9zMFBpa2tl5qN7tKGYECXixUNX9wH8BpXiKRC2wHqhR1XbnaeXABOf+BKDMeW278/wxPd9XRG4UkXUisq6qqsq/P4QZUZpa29lWWc/CSe72V3SxnfOMG9xohkrBc7UwGcgA4oALhvq+qvqwquapal5aWtpQ386Y4zaV19LRqa73V3TJGjOK0XFRNjnPBJQbzVDLgD2qWqWqbcCzwBIg2WmWApgI7Hfu7wcmATiPJwFHAhvZjGRdI4/cmrndk4iQazvnmQDrt1iIyCgRuUtEHnGOp4lI/hA+sxRY7LyvAEuBYuAN4JPOc1YAzzv3X3COcR5/XVWtZ88ETGFpNdnOt/lgsTAzhd1VjdQ0tbodxYwQ3lxZ/AZoAT7mHO8H7h3sB6rq+3g6qguATU6Gh4FvAreKSAmePonHnJc8Boxxzt8K3D7YzzZmoFSVgtKa4/0EweL45Dybb2ECxJtNhKeq6pUicjWAqjYNdTSSqt4N3N3j9G7glBM8txn41FA+z5jBKq8+xuGGlqDpr+jSNTKroLSas2eOdTuOGQG8ubJoFZFYQAFEZCqeKw1jhr1CZ8/rYOmv6DIqKoKZ423nPBM43hSLu4GXgUki8nvgNeA2v6YyJkgU7KsmJjKMmePdn4zXU25mChtKbXKeCYx+i4Wqvgp8ArgO+COemddv+jeWMcGhsKyGnInJRIQH3zJqi7JSaGztYHulTc4z/ufNaKglQLOqvggkA98WkSx/BzPGbc1tHRQfqA26zu0uxyfnWVOUCQBvvi49BDSJyHw8o5F2Ab/1aypjgsCWA7W0dQTPZLyeJo2OJTXeJueZwPCmWLQ78xqWA79S1V8BwdeAa4yPFR6fjJfsao7edO2cZ8NnTSB4UyzqReRbwDXAiyISBri3VZgxAVJQWs3ElFjGJsS4HaVXuZkp7DncyNFGm5xn/MubYnElnqGyN6hqJZ6lOH7i11TGBIHC0pqgGzLbU65z1VNoTVHGz7wZDVWpqj9V1bed41JVtT4LM6xV1B6jorb5+B/jYJUzMZkIZ3KeMf7kzWioT4jIThGpFZE6EakXkbpAhDPGLYVBtnhgb2KjwpmVnmjLlRu/86YZ6sfApaqapKqJqpqgqon+DmaMmwpLq4mKCGN2evD/r56bmczG8hraOzrdjmKGMW+KxUFV3er3JMYEkYLSGuZNSCIqIvgm4/WUm5VCU2sH2w/a5DzjP94sJLhORJ4C/kK3NaFU9Vl/hTLGTa3tnWzaX8u1i0Nj7umHk/NqmJOR5HIaM1x587UpEWgCzgMucW5D2c/CmKBWXFFHa3snuVnB3V/RZWJKLGkJ0RTss05u4z/9Xlmo6ucCEcSYYNE1DDVYJ+P1ZDvnmUDwZjTURBF5TkQOObdnRGRiIMIZ44bC0hrGJ8aQnhTrdhSv5WamsO9IE4cbbPcA4x/e7pT3ApDh3FY554wZlgpKq8nNSnY7xoB0NZnZ0h/GX7wpFmmq+htVbXduTwBpfs5ljCsO1TdTXn2MhZNCo7+iy7wJSTY5z/iVN8XiiIhcIyLhzu0a4Ii/gxnjhmBfPLA3MZHhzMlItE5u4zfeFIvrgU8Dlc7tk4B1epthqbC0hshwYe6E0BuCujAzhaLyWpucZ/zCm7Wh9qnqpaqa5twuU9XSQIQzJtAKS6uZnZ5ITGS421EGLDcrhWNtHWyznfOMH3gzGmqKiKwSkSpnNNTzIjIlEOGMCaS2jk6KymuDfj2o3nQterjemqKMH3jTDPUH4P+AdDyjof6MZy9uY4aVv5cc5lhbB0tOSnU7yqBMSI4le8woXi0+6HYUMwx5UyxGqervuo2GehII3t1gjBmk1UUVJMRE8PHpoVksRISLc9J5d9dhjth8C+Nj3hSLl0TkdhHJFpEsEbkNWCMio0VktL8DGhMILe0dvLKlkvNmjyc6IvT6K7rk52TQqfDS5kq3o5hhxpuFBD/t/PvvPc5fBShg/Rcm5P1tx2Hqm9u5ZH6621GGZOb4BKamxbG66ADXhMhCiCY0eLM21ORABDHGTas2HiBlVGTI9ld0ERHyczJ48PWdHKprZmyitRgb3/BmNNSnRCTBuX+niDwrIgv9H82YwDjW2sHarQe5YG46keHBv39Ffy6Zn44qrNlU4XYUM4x485txl6rWi8jpwDLgMeDX/o1lTOC8vu0QTa0dXJIT2k1QXU4am8DM8QmsLrJiYXzHm2LR4fx7MfCwqr4IRA3lQ0UkWUSeFpFtIrJVRD7mdJi/6uz3/aqIpDjPFRF5UERKRKRIRHKH8tnG9LS66ACp8dGcOmWM21F8Jj8nnXX7qjlQc8ztKGaY8KZY7BeR/wWuxDMKKtrL1/Xl58DLqjoTmA9sBW4HXlPVacBrzjHAhcA053Yj8NAQP9uY4xpa2nl92yEunjee8DBxO47PXJyTAVhTlPEdb/7ofxp4BThfVWuA0cB/DvYDRSQJ+Die5ixUtdV53+XASudpK4HLnPvLgd+qx3tAsogMj/YC47q1xQdpae/kkvkZbkfxqcmpcczJSGSVNUUZH/Fmbagm4BBwunOqHdg5hM+cDFQBvxGRQhF5VETigHGq2vV/diUwzrk/ASjr9vpy55wxQ7a66ADpSTHH97EeTvJzMthYVkPZ0Sa3o5hhwJvRUHcD3wS+5ZyKBJ4cwmdGALnAQ6q6EGjkwyYnAFRV8czh8JqI3Cgi60RkXVVV1RDimZGitqmNt3ZUkZ+TTtgwaoLqku902FtHt/EFb5qhLgcuxfNHHVU9ACQM4TPLgXJVfd85fhpP8TjY1bzk/HvIeXw/MKnb6yc65z5CVR9W1TxVzUtLs72ZTP9eKa6krUPJzxleTVBdJo0exfxJyawuOuB2FDMMeFMsWrt/03eajAZNVSuBMhGZ4ZxaChTj2bp1hXNuBfC8c/8F4FpnVNRioLZbc5Uxg7Zq4wEyR48iZ2Lo7V3hrUty0tlyoI49hxvdjmJCnDfF4v+c0VDJIvIFYC3wyBA/9yvA70WkCFgA3Af8CDhXRHbimc/xI+e5a4DdQInzuV8a4mcbw5GGFt7ddYT8nHREhl8TVJeL5jlNURvt6sIMTZ/LfYjnt+gpYCZQB8wAvqOqrw7lQ1V1A5B3goeWnuC5Ctw0lM8zpqeXNlfS0Tl8m6C6ZCTHkpeVwuqiCr6ydJrbcUwI67NYqKqKyBpVnQcMqUAYE0xWFx1galocs9KH0v0WGvJz0rlnVTE7D9Yzbdzw/3mNf3jTDFUgIif7PYkxAXKwrpn39xwlPydjWDdBdbloXjoi2JwLMyTeFItTgX+IyC5nuY1NTl+DMSFpzaYKVAn55ci9NTYxhlMnj2Z10QE8rbrGDJw3+1mc7/cUxgTQqo0HmDk+gZPGjpwmmfycDO78y2a2VtQzOyPR7TgmBHkzg3vfiW6BCGeMr+2vOUZBac2wW96jPxfO9ax9ZXMuzGCF/uL9xgzAi84fy0uG+SionsbER3Pa1DGsLqqwpigzKFYszIiyamMF8ycmkTlmlNtRAi4/J53So01s2l/rdhQTgqxYmBFj7+FGNu2vHfZzK3pz/pzxRIQJL9qoKDMIvRYLEdkjIru73bof7wpkSGN8oau9/uJhsiPeQCWPiuKMaanWFGUGpa8rizzg5G63U4D/BgTY4PdkxvjY6qIK8rJSyEiOdTuKa/JzMthfc4zCshq3o5gQ02uxUNUjqnoEqAbygTeAjwEXq+oVAcpnjE/sPFjPtsr648t2j1TnzhlHVHgYqzdaU5QZmL6aoSJF5N/xrAh7BnCZql6jqsUBS2eMj6wqqiBM4KIRXiwSYyI5c0YaazZV0NlpTVHGe31NytuDZ1e8B4BSIEdEcroeVNVn/RvNGN9QVVZvPMCpk8cwNiHG7Tiuy89J59Xig6zbV80pk0e7HceEiL6KxVo8e1jMd27dKWDFwoSE4oo6dh9u5PNnTHE7SlBYOmsc0RFhrC46YMXCeK3XYqGq1wUwhzF+s7qogvAw4YK5492OEhTioyM4Z+ZY1myq5O5L5hA+DLeUNb7XV5/FA93u39zjsSf8F8kY31FVVm08wOknpTI6LsrtOEEjPyeDww0tvL/7iNtRTIjoa+jsx7vdX9HjsRyMCQEby2sprz424kdB9XTOzLGMigq3ZcuN1/oqFtLLfWNCxqqNB4gKD+O8OdYE1V1sVDhLZ43j5c0VtHV0uh3HhIC+ikWYiKSIyJhu90eLyGggPED5jBm0zk7lxaIKPj49jaTYSLfjBJ38nHSqm9p4d5c1RZn+9VUskoD1wDogEShwjtcDI2cjABOy1u2rprKuecRscjRQZ05PIyE6gtUbbdly07++RkNlBzCHMT63uugAMZFhLJs1zu0oQSkmMpxzZ4/jlS2V/ODyeURF2Lqipnd9jYYKF5H4bseLReTjzs2uLExQa+/oZM2mCs6ZOZa4aG82hByZ8uenU9fczts7q9yOYoJcX18l7ge+1O34j8B/AncBd/ozlDFD9f6eoxxuaB1xmxwN1OknefpzVtuoKNOPvr5yLcWz2myXGlW9REQEeNu/sYwZmlUbDxAXFc7ZM8e6HSWoRUWEcf6ccazZVElzWwcxkTZ2xZxYn6OhVLW92/E3AdSzEH78iV9ijPvaOjp5eUsl584eZ3/8vJCfk0FDSztvbremKNO7vopFVPe+CVX9K4CIJAG2GpsJWu+UHKamqW3E7og3UKdNHcPouKjjm0MZcyJ9FYtHgKdEJLPrhIhk4em7eNTfwYwZrFUbD5AYE8EZ01PdjhISIsLDuGDueF7beoim1vb+X2BGpL42P/op8ALwjogcEZEjwN+AVar6X4EKaMxANLd18OqWg5w/ZzzREdYE5a38nHSOtXXw+rZDbkcxQarPgdWq+mtVzQSygWxVzVLVhwKSzJhBeGtHFfUt7Vwy35qgBuLUyWNIS4i2HfRMr7yahaOq9apa7+8wxgzV6qIKRsdFcdrUMW5HCSnhYcJFc8fzxvZDNLRYU5T5V65N2XQm/RWKyGrneLKIvC8iJSLylIhEOeejneMS5/FstzKb4FbT1Mra4oNcOHc8EeE2G3mg8udn0NLeyavFlW5HMUGorxnc/r6OvxnY2u34fuBnqnoSUA3c4Jy/Aah2zv/MeZ4x/+L+l7fR2tHJZz+W5XaUkLQoM4UpaXE8sHYnzW0dbscxQaavr1+Pish7IvIjETlLRHy2ZoKITAQuxhlV5Uz0Owd42nnKSuAy5/5y5xjn8aXO8405bv2+o/zxn2VcvySbmeMT3Y4TksLChHuXz2XfkSZ+9UaJ23FMkOlrNNRFwFnAm8DlwHsi8qyI3Nh9OO0gPQDcBnQtpD8GzwzxrsbScmCCc38CUOZkagdqned/hJNrnYisq6qyyUUjSVtHJ3c8t5n0pBhuWTbd7Tgh7bSTUrl84QR+/dYuSg5ZN6X5UH+joZpV9WVVvVlV84Cv41ki5Jci8s/BfKCI5AOHVHX9YF7fR9aHVTVPVfPS0tJ8+dYmyD3x971sq6znnkvn2KKBPnDHxbMYFRXBHc9txrNggzED7OBW1T2q+j+qeilw+iA/cwlwqYjsBf6Ep/np50Byt6auicB+5/5+YBKA83gSYLu1GAAO1BzjZ2t3sGzWWM6bbUuR+0JqfDS3XziT9/cc5en15W7HMUFi0ENGVLV1kK/7lqpOdPbLuAp4XVX/DXgD+KTztBXA8879F/hwD/BPOs+3rzsGgO+u2kKnKndfMgfryvKdK/MmkZeVwn1rtnK0cVC/6maYCabxhd8EbhWREjx9Eo855x8DxjjnbwVudymfCTJriw/yypaD3Lx0OpNGj3I7zrASFibc94l51De388M1W/t/gRn2+m3gFZEY4CTnsERVm3314ar6Jp4OdFR1N3DKCZ7TDHzKV59phoem1nbufmEL08fF8/kzJrsdZ1iaPi6BL3x8Cg+9uYsrFk1k8RSb6DiS9TXPIkJEfoxnZNJK4LdAmYj8WEQiAxXQmBP5xesl7K85xr2XzSPSJuD5zVfPmcak0bHc8dwmWtpt7sVI1tdv2U+A0cBkVV2kqrnAVCAZsIUEjWt2HKznkb/t5lOLJnLK5NFuxxnWYqPC+d7yueyqauTht3a7Hce4qK9ikQ98ofuaUKpaB3wRuMjfwYw5kc5O5Y7nNhEfE8G3LprldpwR4ewZY7l4Xjq/eKOEvYcb3Y5jXNJXsdATjTpS1Q7ARiMZVzxdUM4He6v59oWzGB0X5XacEeM7l8wmOjyMu563uRcjVV/FolhEru15UkSuAbb5L5IxJ1bd2MoP12wlLyuFTy6a6HacEWVcYgzfOH8Gb+88zAsbbUe9kaiv0VA3Ac+KyPVA12zrPCAWz/IfxgTUj17aRn1zO/dePpewMJtTEWjXLM7imYJyvr96K2dNH0vSKBvnMpL0tTbUflU9FfgesNe5fU9VT1HV/b29zhh/WLf3KE+tK+OGMybbQoEuCQ8T7rt8HkcbW7j/FWtcGGl6vbIQkZOBVFV9CXi92/kL8cPaTsb0pmuhwAnJsdy8dJrbcUa0uROS+NySyTz2zh6uyJ3IoqwUtyOZAOmrz+J+oPgE54vxDKs1JiAef2cP2w96FgocFWULBbrt1nOnk54Uwx3PbaKto7P/F5hhoa9ikaCq+3qedM6l+i+SMR8qr27igbU7OXf2OM61hQKDQlx0BPdcOodtlfU8/s4et+OYAOmrWPR1fWkL8ZiAuOcFz8XtPZfOcTmJ6e78OeNZNmscD6zdSXl1k9txTAD0VSzWisgPuu9KJx7fo1sfhjH+8tctlazdepBblk1jQnKs23FMD99dPgcRuPv5LTb3YgToq1h8HZgClIjIMyLyDLATmO48ZozfNLa0c88LW5gxLoHrT7eFAoPRhORYvrZsOq9tO8QrWyrdjmP8rNfeQlVtBK4WkSlAVxvAFmd1WGP86sHXdnKgtpmnr15oCwUGsc8tyebZwv3c80IxS05KJSHG5l4MV32tOnsNHF86/KiqruoqFCLy5QDlMyPQtso6Hntnj2cDnmxbKDCYRYSHcd/lczlY38x//3WH23GMH/X1le3Wbvd/0eOx6/2QxRg6O5U7n9tMQkwEt1840+04xgsLM1O45tQsfvuPvWwqr3U7jvGTvoqF9HL/RMfG+MSf15exbl81375oFim2UGDI+M8LZjAmPppvP7eJjk7r7B6O+lx1tpf7Jzo2ZsjKq5v44UvbOCV7tC0UGGISYyL5Tv5sNu2v5X//tsvtOMYP+poOO1NEivBcRUx17uMcT/F7MjOi1DW3ccMT6+joVH54xTy6jdg2ISI/J52XN1fyk1e2kz0mjovmpbsdyfhQX8XCdpYxAdHe0cmX/1DIrqoGVl5/ClPT4t2OZAZBRPjvT8+novYYtzy1gXGJ0SzKsgEKw0Vfq87uO9ENKANOD1xEM5ypKves2sLfdlTxg8vnsuQkW0kmlMVEhvPoipPJSIrh8yvXscd21hs2+ho6mygi3xKRX4rIec7s7a8Au4FPBy6iGc4ee2cPT75Xyn+cOZUrT850O47xgdFxUTzxuVMQEa77zT850tDidiTjA311cP8OmAFsAj4PvAF8ErhMVZcHIJsZ5v66pZIfrNnKhXPHc9v5M9yOY3woOzWOR67No7K2mRtWruNYa4fbkcwQ9VUspqjqdar6v8DVwGzgfFXdEJBkZljbVF7LzX/aQM7EZH766QW2890wtCgrhZ9ftYCN5TXc8lShDakNcX0Vi7auO6raAZSrarP/I5nhrqL2GDes/IDRcVE8cu0iYqPC3Y5k/OSCuencdfFsXtlykHtfPNH2OCZU9DUaar6I1Dn3BYh1jgVQVbW9Lc2ANbS0c/0T62hq7eCZL57K2IQYtyMZP7v+9MmUVTfxm7/vZWLKKG6whSFDUl8LCdrXPeNT7R2dfPWPhew4WM/j153MjPEJbkcyAXLnxbM5UHOMe18sZkJyDBfMtTkYocaW8zQBc++LW3l92yG+e+kczpye5nYcE0DhYcIDVy5kwaRkbv7TBtbvq3Y7khkgKxYmIJ74+x6eeHcvnz99MtcsznI7jnFBbFQ4j16bx/ikGD6/8gObgxFiAl4sRGSSiLwhIsUiskVEbnbOjxaRV0Vkp/NvinNeRORBESkRkSIRyQ10ZjM0r287yPdWF7Ns1ji+dZEtDDCSjYmP5onPnQLA537zT442trqcyHjLjSuLduDrqjobWAzcJCKzgduB11R1GvCacwxwITDNud0IPBT4yGawig/U8ZU/FDI7I5EHr15AuA2RHfEmp8bx6IqTqaht5vMrP6C5zeZghIKAFwtVrVDVAud+PbAVmAAsB1Y6T1sJXObcXw78Vj3eA5JFxHrHQsDBumZuWPkBCTGRPLbiZEZF9TX4zowki7JSeODKBRSW1XDLnzbYHIwQ4GqfhYhkAwuB94FxqlrhPFQJjHPuT8CzHlWXcuecCWJNre3csPIDao+18dh1eYxLtCGy5qMunJfOHRfN4uUtldy3ZqvbcUw/XPuqJyLxwDPALapa131JalVVERnQVw0RuRFPMxWZmbbGkJs6OpWv/nEDxQfqeHRFHnMyktyOZILUDadPprz6GI+9s4cJybFcb3MwgpYrVxYiEomnUPxeVZ91Th/sal5y/j3knN8PTOr28onOuY9Q1YdVNU9V89LSbFimm364Zitrtx7kO/mzOWfmuP5fYEYsEeGu/NmcN3sc33+xmJc3V7odyfTCjdFQAjwGbFXVn3Z76AVghXN/BfB8t/PXOqOiFgO13ZqrTJB58r19PPrOHq47LZvrlti3RNO/8DDh51ctZP7EZG7+UyHvlhx2O5I5ATeuLJYAnwXOEZENzu0i4EfAuSKyE1jmHAOswbMsegnwCPAlFzKbflTVt/DNp4u46/nNnD0jjTsvtiGyxnuxUeE8tiKPCSmxfObR9/nO85tpaGl3O5bpRlSH3yiEvLw8XbdundsxRoSW9g4ef2cvv3qjhOa2Dq47LZtbz5tuI5/MoDS2tPNff93OE+/uJSMplvs+Mc9m+weQiKxX1bwTPmbFwgyGqvLKlkruW7ON0qNNLJs1lm9fNIsptiWq8YH1+45y29NF7Kpq5IrcidyVP4vkUVFuxxr2+ioW9vXPDFjxgTq+t3oL7+0+yvRx8fzuhlM4Y5p9+zO+syhrNC9+9Qx++XoJD721i7d2VPH95XO4cJ5NsXKLXVkYrx1uaOG//7qDpz4oJSk2klvPnc7Vp2QSEW5LjBn/2XKgltueLmLLgTounDue7y6fY0vb+4ldWZghaW3v5Il39/CL10o41tbBdadN5ual00gaFel2NDMCzMlI4vmblvDw27t5YO1O3t11hLvyZ3NF7gS6z88y/mVXFqZXqsrarYf4wYvF7D3SxNkz0rjj4tmcNNb6JYw7dlU18M2ni1i3r5qPT0/jvsvnMjFllNuxhg3r4DYDtq2yjntXb+WdksOcNDaeOy+exVkzxrodyxg6O5XfvbeP+1/ehgDfvHAm15yaZfu4+4AVC+O1o42t/PTV7fzh/VISYiL52rJp/NviLCKtX8IEmfLqJr717Cbe3nmYk7NT+NEVOUy10XhDYsXCeOXtnVV8+Q+FNLS089nFWdyybJoNVzRBTVV5en05319dTHN7J19bNp0vnDHZBl0MknVwm349+d4+7n5hC9PGxvOLqxcybZztj22Cn4jwqbxJnDkjje/8ZQv3v7yN17cd5Ff/lmsjpnzMyu8I19GpfG9VMXf+ZTNnTk/j6S+eZoXChJyxCTH8+rOL+PlVC9i8v478B99h/b6jbscaVqxYjGANLe3c+Nt1PP73PXxuSTaPXJtHfLRdbJrQtXzBBJ676TRio8K56uH3+N0/9jIcm9rdYMVihDpQc4xP/fofvOnMjL37kjm25akZFmaOT+SFL5/OGdPSuOv5LXzjz0W2dasPWLEYgTaW1bD8V3+n/GgTj193Mp/9WLbbkYzxqaTYSB69No9blk3jmYJyrnjoXcqONrkdK6RZsRhhXtpUwZUP/4PoiDCe+dJptqKnGbbCwoRblk3nsRV5lB5t4pJfvsPbO6vcjhWyrFiMEKrK/7xZwhd/X8Ds9ET+ctMSpltHthkBls4ax6ovn864hBhWPP5P/ufNEuvHGAQrFiNAa3sntz1dxI9f3s4l8zP4wxcWkxof7XYsYwImOzWO5246jYtzMvjxy9v5jyfXU9/c5naskGLFYpirbmzls4+9z5/Xl/PVpdN48KoFxESGux3LmIAbFRXBg1ct4K782azdeojlv/o7JYfq3Y4VMqxYDGO7qxr4xEPvUlhawwNXLuDWc6fbKp1mRBMRbjh9Mr///KnUHWtj+S//zsubK9yOFRKsWAxT/9h1hMv/511qj7Xxhy+cymULJ7gdyZigsXjKGFZ95XSmjUvgP54s4P6Xt9HRaf0YfbFiMQz93wdlfPax90lLiOYvX1pCXvZotyMZE3TSk2J56t8X85lTM3nozV2sePyfHG1sdTtW0LJiMYwcrGvmBy8Wc9szRSyeMoZnvngamWNsrX9jehMdEc59l8/j/ivm8c+9R7nkF+/w0qYKGlra3Y4WdGxthxDW3tFJQWkNb24/xJvbqyiuqAPgM6dm8t1L59iy4sZ46cqTM5mVnsgXnyzgi78vIDJcODl7NGfPGMtZM9I4aWz8iO/vsyXKQ8yhumbe3FHFW9ureHtnFXXN7YSHCYuyUjhrRhpnzxjLrPREt2MaE5Ja2ztZv6+aN7cf4o3th9hxsAGACcmxnD3T8/v1saljGBU1PL9n234WIay9o5MNZTW84Vw9bDnguXoYmxDNWTPSOGvGWJaclEpSrO2HbYyv7a855ikc26p4d9dhmlo7iIoI49TJnquOs2eOZXJqnNsxfcaKRYipqm/hrR1VvLH9EG/v+PDqITczmbOcy+LZ6Ykj/rLYmEBqae/ggz3Vzhe3Q+yqagQge8yo47+Xi6eMCel5TFYsgpiqUl59jILSagpLa/hg79HjVw9pCdGcOd1z6Xv6NLt6MCaYlB5p4s0dniv+d3cdprmtk5jIMBZlpbAoazR5WSkszEwmISZ0fm+tWASR5rYONu2vpWBfNQWl1RSU1lBV3wJAbGQ4OROTOGNaKmfNGMvs9ETbhN6YENDc1sF7u4/w5vYq/rnnKNsq6+hUCBOYMT6RRVnJ5GWNZlFWChNTYoO2VcCKhUt6XjUUlFZTfKCOdmfyT9aYUSyclExuVgq5mSnMHJ9gewcbMwzUN7exoayGdXurWb+vmsLSahpbPXtqjEuM/sjVx+yMxKAZuWjFIgBUlaONrZQcamBDWU2vVw1dhWFhZrIt5mfMCNHRqWyrrGP9vurjBWR/zTHA87dh/qQkFmWlkJc1mjkZiaQlRLty9WHFwocaW9rZc7jxI7fdhxvZU9VAXfOHE3nsqsEY05eK2mMfKR7FFXXHlxyJiwonOzWO7NQ4pqTGkT0mjslpcUweE0dKXJTfMg2LYiEiFwA/B8KBR1X1R709d6jForW9k7LqJvZUdSsGhxvYc7iRg3UtH3luRlKM5z9iahyTU+OZkhrHvIlJdtVgjBmQptZ2NpTVsPNgw/EvonuPNFJ2tInuy1YlxUY6f288t+xUTxHJTh015M70kC8WIhIO7ADOBcqBD4CrVbX4RM8fbLHYvL+WL/+hgLLqYx9ZVGx0XNRH/uNMSfVU+azRccRGhe4wOWNM8Ov68rq3R4vG3sONHKht/shzU+OjuWxBBnfmzx7UZ/VVLEJlGuIpQImq7gYQkT8By4ETFovBGhMfxdwJSVw6P8O5Wohn8pg4kkaFztA3Y8zwEhURxtS0eKamxf/LY8daO9h31FM4djsFJCM51i85QqVYTADKuh2XA6d2f4KI3AjcCJCZmTmoD0lPiuWXn8kdZERjjAms2KhwZo5PZOZ4/y/xM2x6XFX1YVXNU9W8tLQ0t+MYY8ywEirFYj8wqdvxROecMcaYAAiVYvEBME1EJotIFHAV8ILLmYwxZsQIiT4LVW0XkS8Dr+AZOvu4qm5xOZYxxowYIVEsAFR1DbDG7RzGGDMShUozlDHGGBdZsTDGGNMvKxbGGGP6FRLLfQyUiFQB+9zO0YdU4LDbIVxkP7/9/PbzB6csVT3hRLVhWSyCnYis6239lZHAfn77+e3nD72f35qhjDHG9MuKhTHGmH5ZsXDHw24HcJn9/COb/fwhyPosjDHG9MuuLIwxxvTLioUxxph+WbFwkYh8XURURFLdzhJIIvITEdkmIkUi8pyIJLudKRBE5AIR2S4iJSJyu9t5AklEJonIGyJSLCJbRORmtzO5QUTCRaRQRFa7nWWgrFi4REQmAecBpW5nccGrwFxVzcGzt/q3XM7jd84+8r8CLgRmA1eLyOA2Sg5N7cDXVXU2sBi4aYT9/F1uBra6HWIwrFi452fAbcCIG2Ggqn9V1Xbn8D08m1kNd8f3kVfVVqBrH/kRQVUrVLXAuV+P5w/mBHdTBZaITAQuBh51O8tgWLFwgYgsB/ar6ka3swSB64GX3A4RACfaR35E/bHsIiLZwELgfZejBNoDeL4gdrqcY1BCZj+LUCMia4HxJ3joDuDbeJqghq2+fn5Vfd55zh14mid+H8hsxj0iEg88A9yiqnVu5wkUEckHDqnqehE5y+U4g2LFwk9UddmJzovIPGAysFFEwNMEUyAip6hqZQAj+lVvP38XEbkOyAeW6siY7DPi95EXkUg8heL3qvqs23kCbAlwqYhcBMQAiSLypKpe43Iur9mkPJeJyF4gT1WDdRVKnxORC4CfAmeqapXbeQJBRCLwdOYvxVMkPgA+M1K2BxbPN6OVwFFVvcXlOK5yriy+oar5LkcZEOuzMG74JZAAvCoiG0Tk124H8jenQ79rH/mtwP+NlELhWAJ8FjjH+W++wfmWbUKEXVkYY4zpl11ZGGOM6ZcVC2OMMf2yYmGMMaZfViyMMcb0y4qFMcaYflmxMMYY0y8rFsYYY/plxcKYABGRb4rI70TkFRHZJyLfcDuTMd6yYmFM4MzDs+LoRcAy4Fp34xjjPVtI0JjAmQd8QlU7RKQDOOp2IGO8ZVcWxgSAs+Jqmqruck7lAJtcjGTMgFixMCYwZvDR7TQXALb5lQkZViyMCYx5QFG34wVYsTAhxFadNcYFIlICzFPVY25nMcYbdmVhTICJSDLQaoXChBK7sjDGGNMvu7IwxhjTLysWxhhj+mXFwhhjTL+sWBhjjOmXFQtjjDH9smJhjDGmX1YsjDHG9Ov/ARxnKtscEk49AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# jet energy response\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "h_jer = output['r_pt_ptveta'].sum('dataset').sum('eta').sum('pt')\n",
    "ax = hist.plot1d(h_jer)\n",
    "ax.set_xlim(0,2)\n",
    "plt.legend(['QCD'])\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# jet energy response vs eta\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "h = output['r_pt_ptveta'].sum('dataset').sum('pt')\n",
    "\n",
    "xaxis='eta'\n",
    "yaxis='r'\n",
    "xaxis = h.axis(xaxis)\n",
    "yaxis = h.axis(yaxis)\n",
    "xoverflow='none'\n",
    "xedges = xaxis.edges(overflow=xoverflow)\n",
    "xcenters = xaxis.centers(overflow=xoverflow)\n",
    "\n",
    "\n",
    "\n",
    "ax.plot(xcenters, np.average(h.values()[()], axis=1))\n",
    "ax.set_xlabel(xaxis.label)\n",
    "ax.set_ylabel(yaxis.label)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# your code here\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## Bonus: Plot  response correction for $p_{T}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## Bonus: Plot JEC for separate $\\eta$ and $p_{T}$ Bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " eta bins:\n",
      "\n",
      "Bin #0: (-inf, -5)\n",
      "Bin #1: [-5, -4.5)\n",
      "Bin #2: [-4.5, -4)\n",
      "Bin #3: [-4, -3.5)\n",
      "Bin #4: [-3.5, -3)\n",
      "Bin #5: [-3, -2.5)\n",
      "Bin #6: [-2.5, -2)\n",
      "Bin #7: [-2, -1.5)\n",
      "Bin #8: [-1.5, -1)\n",
      "Bin #9: [-1, -0.5)\n",
      "Bin #10: [-0.5, 0)\n",
      "Bin #11: [0, 0.5)\n",
      "Bin #12: [0.5, 1)\n",
      "Bin #13: [1, 1.5)\n",
      "Bin #14: [1.5, 2)\n",
      "Bin #15: [2, 2.5)\n",
      "Bin #16: [2.5, 3)\n",
      "Bin #17: [3, 3.5)\n",
      "Bin #18: [3.5, 4)\n",
      "Bin #19: [4, 4.5)\n",
      "\n",
      " pt bins:\n",
      "\n",
      "Bin #0: (-inf, 5)\n",
      "Bin #1: [5, 10)\n",
      "Bin #2: [10, 15)\n",
      "Bin #3: [15, 20)\n",
      "Bin #4: [20, 25)\n",
      "Bin #5: [25, 30)\n",
      "Bin #6: [30, 35)\n",
      "Bin #7: [35, 40)\n",
      "Bin #8: [40, 45)\n",
      "Bin #9: [45, 50)\n",
      "Bin #10: [50, 60)\n",
      "Bin #11: [60, 70)\n",
      "Bin #12: [70, 80)\n",
      "Bin #13: [80, 90)\n",
      "Bin #14: [90, 100)\n",
      "Bin #15: [100, 120)\n",
      "Bin #16: [120, 140)\n",
      "Bin #17: [140, 160)\n",
      "Bin #18: [160, 180)\n",
      "Bin #19: [180, 200)\n",
      "Bin #20: [200, 250)\n",
      "Bin #21: [250, 300)\n",
      "Bin #22: [300, 350)\n",
      "Bin #23: [350, 400)\n",
      "Bin #24: [400, 450)\n",
      "Bin #25: [450, 500)\n",
      "Bin #26: [500, 600)\n",
      "Bin #27: [600, 700)\n",
      "Bin #28: [700, 800)\n",
      "Bin #29: [800, 900)\n",
      "Bin #30: [900, 1000)\n",
      "Bin #31: [1000, 1500)\n",
      "Bin #32: [1500, 2000)\n",
      "Bin #33: [2000, 3000)\n",
      "Bin #34: [3000, 4000)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print('\\n eta bins:\\n')\n",
    "for i in range(EtaBinNums):\n",
    "    print('Bin #' + str(i) + ': '+ str(EtaBins[i]))\n",
    "    \n",
    "print('\\n pt bins:\\n')\n",
    "for i in range(PtBinNums):\n",
    "    print('Bin #' + str(i) + ': '+ str(PtBins[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "coffea-jerc",
   "language": "python",
   "name": "coffea-jerc"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}