02 Comparing simulation yields to data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.22/00\n"
     ]
    }
   ],
   "source": [
    "import ROOT as r\n",
    "import numpy as np\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing simulation yields to data\n",
    "\n",
    "Often we want to compare the event yields of simulation samples to actual data, to see if number of events you get from data matches what your simulation predicts.\n",
    "\n",
    "To relate the two numbers (MC yields v.s. data), we just have to remember the expected number of events from any process at LHC is given by\n",
    "\n",
    "$$ N = L \\times \\sigma \\times Br$$\n",
    "where \n",
    "\n",
    "$\\begin{align}\n",
    "L      &= \\text{integrated luminosity of data collected} \\\\ \n",
    "\\sigma &= \\text{cross section of the process}\\\\\n",
    "Br     &= \\text{Branching ratio of the process into specific final states}\n",
    "\\end{align}\n",
    "$ \n",
    "\n",
    "For example, if a process A has a cross section of 1 fb, and you have collected 10 /fb of data, the expected number of events in data would be 10 events.\n",
    "If you have generated 50 simuated events (more than what you expect in data), you simply need to scale each event with a weight of w,\n",
    "\n",
    "$\\begin{align}\n",
    "w &= (L \\times \\sigma )/ N_{sample} \\\\\n",
    "   &  = 10 \\times 1 / 50\\\\\n",
    "    & = 0.2\n",
    "     \\end{align}\n",
    "$\n",
    "\n",
    "Let's see how to do this in ROOT in practice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filling a histogram with a weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = r.TCanvas('c1','c1',800,600)\n",
    "h = r.TH1F('h_test','Test',10,0,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h.Fill(5,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h.GetEntries()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dUZKjurYtUHhx+wXuzLnNAJpxd2eAlvE+VmwFBztddpawAY3xUZGFSVt2Oq2ZS0Kql2WpAAD4O//v2w0AALgCoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIIP/+XYDgCubpunFM9u23bEdAPsTqoC99H0/DMOLJ9vdATg7w3/A9zVNs9M9t22rBgZ8Ru2vQ+BjUu3qY588dV1/8uGAkqlUAQBkYE4VcETTNMUk9+fjd3HONE0PT3t9mjzA3xOqgMOJMbsQw4VN09wnpLZt53n+6bT1nRgEBD7A8B9wLCkMNU3TdV3MYZ/neVOIqus6EtVPp6WDcc5+c+EBgkoVcCARiTZ1qWmabrdbKkpVq3G9dfFpc1rf99W/Ec04IPABKlXAUUzTFJFok4Hato06U6pCxQmb4lPbtl3XdV33ibYC3BGqgKN4GJVClJ3Wxar47yZ+9X0fZwJ8nuE/4Fjmeb4PRvfhKWam3263qqqapun73iKfwHcJVcBRpOT0yuY2y7Kkq//meZaugK8z/AccS9M0y8/WZ07TtCzLOI5pHlVKVwCfJ1QBR/G7ClPbtn3fR7r6m/sB+EtCFXAUEYY2s9FD3/d1Xae01LZtXdebqVdx9d/ejQT4iVAFHMU6M62PT9MUs6w2Kep+6lUcUakCvqK2bwPwMemqvZ8+eWIBz/g6yk5p8ar1iqD3p1WrjLW+81j8s+u653sIAvw9oQr4nD+Gquq/A1Nyv/dfuquNcRzX4Wm9A6CPO2BXQhVwRNM0pRT1pMj0ymlp0NC6oMCuhCoAgAxMVAcAyECoAgDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyOB/vt0AADir9ZL9XNJby3kKVQDwe9bQvrB3Q7PhPwCADIQqALimtm0NUH6SUAUA1zTP8998+zRNdV3n2on8yc7ol2FOFQCwu79MeKegUgUA1zdN0zRNn7m3vI91IkIVAFxZjOLdbrfb7fb6cF7btrfbraqqYRjWE7Ninla6t014quv6/rGiAemEXOOJByRUAcCV3W63ruvGcey6rqqqYRhe+a6+7+P8pmnGcYyDdV3P8xxH4tbb7ZZyVSSn+8dq2zbdwziOV55ZtQAAv3LwbjQ6+q7r0pGmaaqqGsfxlW+PJJS+PXLS+t7iIZqmWT9cuinOT491xtTxboNVqgDgytbDbX9TJYqK1GbwrmmazQz0VLjq+35ZlivXpe4IVQDAn0V4qv9bHIwglQYE67pu2/bCc6d+IlQBAK9qHolyVJSmIlrN8xwz3Iu6DNA6VQDAn8VI3x9DUt/3UaPq+34Yhr7vy8lVKlUAwJ+lctT6YAwCVv+um7CeQRVnlrDmZ6JSBQD8aD3xfBiGtEpCVVWxkFWM98WReZ77vl/Hr7i1EPXy71WOAMBb6vrQ3WjUkNYtjGD0+mJRadHOdCebHZq7rku1q2maImYlTdOkTNa2bVSt1t9ycO/+fA/9bgCAIzt4qMoiUtE6hMUuND9tkPzk1rjpLImqEqoA4GNKCFUle/fna04VABTnj8N/ZS409ZfKClWbkWAAKJPAtIey6pbqtABkpFu5tnd/vtapAgDIoKzhPwDIy8QSEqEKAH7P8N+FvZuYDf8BAGQgVAEAZGD4DwB2kXZouffiLjGv3H+WuyKLsq4Fde0rABk96VZil72fvvHFzij2e3ny6K/fFb9gRXUAOIqmaX69zOYfNz9umubXDSPUdf369tJ/JFQBwF5+2nU4iyfDi7wi8u7zcuBbipuoXj/17dYBUJC+76Nf7/u+ruv1dntp9DCdE7dGAogQkG5a32Hcel8ei2/cPEqx4oV6Mj77S0tJSnu+AOzqSbfSdV1VVV3XPf/2JI3lxbesh/aaprk/OR15eG9hHMe4Kd3b5lHKNI7j+lV68lK8GxuKq1QBwMcMw9De2RSKIiFN0xS9eJRPpmmKWDaO43qYL528eaCoXaUUFbkhPdA8z1VV3T9Kmdq2jVcpXuGMzKkCgB1FoHniralRP508z3PTNGluUNu2TdP89NAPYxl/T6UKAPbycGgpe6CJO5zneV0Pi0QVN8Wo33pCldWt9nDNUFXXtVl4AJRmXlkfT4OJ8zwPw6CX3MkFQ1VcxKewCUBRYrrVRqpI9X2/nm41DIOOMrurhSrRG4DSRHLaVKfSUkHTNK1LU23bRtVKqMruUqFqmqaSL2cA4GgeXv3Xtu3rgebFMyMnxT2n1SxjKlV8nUpTqa80rSq7S139d7vdqqoaxzG+AICv++PVfz+J1SnD8qcd6KIQNQxD6gGbpkmBrOu69U1xRKjK7jqVqrREx/PTnq+obrF1KFBd//PtJnBBaQ7TQ9FnpRUmk/WRWE6p67ro2p6fnB5xHMdYrWpd4ko3xb0ty2K2zB7e2375sGI5/67rYgn/2+22TujJu9tNAyWo63+W5T/fbgWnpFu5tnd/vleoVMXw8N/sBA4A8JeuMKdqfUVDOhhroMWOid9oFABQliuEqnA/E3Ce53QFBADArq4QqjZzp57MqQIA2MkV5lQBAHydUAUAkMEVhv82YmGPb7cCACjLBUMVAHxdmtf78Hqp57dyUmWtWmaVNuCexT/5tSfdStqK4+EJz28lr7iCbXPwlQvaSlz88y32ogHgk+57bhenX1VxoerJTkz+XAAgu/vdPuz/8WGRYmPTw2SPaFtcqAKAj2ma5uHa1E3T3J8cu4DEdiCbm2It6xhUadt2s1lynN/3fdwqtG3Ey/WJ6WvPKzcXU9rzBV5RVf/37SZwVk+6lehku66Lf9PxcRzT8fW333fQqbKSTl5HsYffmE5YPyLxsozj2HVd13WbktUT78aGskKGUAXcE6r4tT+GqviiaZp0PDr49QnLqteP/0bwSt+1SVGboJbi1P1DE+4D6/rlev6Nbz2Q4T8A2NFmBPDh2F8cTONTbduuv2scx4hZ6db7RzH5/Y9SDI3Xdo/RQOtUAcCO+r6/3W79v6q7ieqRhzbdfCSqmEoVx6d/3U/S4rlxHNev7TRNdV3v8TIKVQCwoxSJqj/NmP6pm+/7fhiG+LppmoeT33ni/gWP1zAya8YHMvwHAPtKMein6/6qH6aWx4V+kajiyDRNLu5718fGRoUqANhXxKBYYvo+EkWxJNWiQlqSOgJBugCwMn3qTbGc+qYiFRk3+7QqoQoA9rXuvB925JGZoi6VxqSiprWeUFWthgJFqxfFCzjPc9/365d3nVOzeetawbPzagD3LKnArz3pODbdSiSkJwsfbPr4+1UYkvWyVff38/AIm+7+xXW83n0Zy9pg2IbKwD0bKvNr2buVn2ayp7pUuulzq4RfRSpTvf6ivfvzLStkCFXAPaGKX9OtXNu7P19zqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAy+J9vN+DTYn/Kn1jDDYC80rZ9b63lnb63+nnZ9M32f1ZXf65t29hHuWmanXZOLGspWEvfAvesqM6vPe9W0ubHa+M4rtNP2t/3ocgBDx/i4Z1XVdV1Xd/3f2h3ee5LKq/kASuqA8D3TdMUoWccx9htNzZCvt1uqUzS9/36v7/QNM04juM4dl0X9z8Mg5LVRrwgaRPleKH2iJ5CFQDkF332ui7V9/04jlXW7rz9V9/3fd9HWWWe552Gt04qRv3Syx5fPKzz/SWhCgDyi458UzRq27Zpmvg6jd9FHoqDMRpY13Xbtr8LRg9zW9/3KXttDm6+fXPOBcTLmF72XRU3UR0APqBpmnmeI6NsJlFtvoj4FUdut1s6mL5+SzxWus/qv6cTzfM8DEOqn83zvI50fd+vv/HC4qeT/W5VqgAgv8hMkY3quu77flN5mqYpJveM4xg3RYrazMH6SxGe0n2u61hx07pV8fXFKlXh4Tyz7IOkQhUA7CKCUQw8DcOQ0tXDk+N413XrOVh/8+gp1TVNk+4zxh/XQ5Prgs3DIctr+MwkM6EKAPYSBaooEaV09XwNhVwPnWZlxShkEskpblpfB5dSXa4GHF/2+ChUAUB+ac3PEBEnXZ338Pz7g7mmV88r6+MRpOKhrzr2d1+Qu/9vLkIVAGQWU84fBpTISfcR6n6GU/Wrvj8edJ3GmqZZ7qyLNLEEQwwUvvtwp7B5XvtdD1hcqKqf+nbrALiCh9WR8NO8pYeh6l1pxdH1/jabZmz6u7QkaXXFMlXYTMzf78kWF6ru0/rat1sHwEVEWEnX/U3T1Pd9pJnNvKUUgNIqDHH+K3/qxwytUNd1xIX1/cfX6T4jWKyLNOtscckp6lVVtW3bdd08z+kl2mwWlM3zkHExpT1f4BVV9X/fbgJn9bxbeTjpO22Wsvy7wMG6O96MScV/X7/z2LXm+Zn3o4HxKOuGXVVs6fP6+e/GhrI2GLahMnDPhsr82ivdyma6+v0Jabnzzfl5SynrAUFe9G5sKCtkCFXAPaGKX9OtXNu7P9/jzqlK+x/FFkjPz7yfci6MAwCfdNC9/9Lk/Bjojcll/hoA4GhcOU5y0FAViSqlqNjK+6ctu+PgXjP5AeAH/tp/rrSCyEGfbV3XTdOsI1T8KfCwtbHo/itPpLSfLvAKc6pgJ6V1uwetVK2vMq1eW/w0bQiwvoYCAOAzjh4hY1GyWB/2pwG+h+PZD59XaZEZeIVKFeyktG73uFf/hWEYIlE1TfO8/pRWLYuC1k8nP9+mxg42AMDvnCBCxtL+sdHjTxPVNxHqpwlYpUVm4BUqVbCT0rrd0zzbJxPV78XU9fvhwtJ+usArhCrYSWnd7hGH/2LLyU1RKgb1niypAADwRUeMkLHy54tLKrx1cmmRGXiFShXspLRu94iVqhizm+e57/tYKCGOpCUVopQVB386+eH23QAAOzlohEzb1CTrWlQssL4+srk6r+u6WItho7TIDLxCpQp2Ulq3e+hnmzLTK4t5pjLVk5NL++kCrxCqYCeldbuFPdvCfrrAK4Qq2Elp3e4R51QBAJyOUAUAkMFBN1Tez/MNZ4qqUgIAGRUXqsQmAGAPhv8AADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyKG5FddvUAAB7KC5UiU0AwB4M/wEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGRQ3IrqtqkBAPZQXKgSmwCAPRj+AwDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyKC4FdVtUwMA7KG4UCU2AQB7MPwHAJCBUAUAkMF1QtU0TW3b1nVd13Xbtt9uDgBQlovMqZqm6Xa7VVXVNE1VVfM813Vt+hQA8DEXqVRFolqWZZqmaZq6rquqSr0KAPiYi5Rz6rpummaapvWR6u5aP+Ur4F5d/7Ms//l2K+CCSut2LzL8N47j+r+RrmIoEADgAy4SqtJIX9/3VVUNw5C+BgD4gKvV5dKC6ZvRwM2tv3CxFwpIDP/BTkob/rvIRPVkWZZxHJummef54UT15bc+/lQAgDO5bIQ0UR14kUoV7KS0bvcKlaq+7+u63gz2xSz1+xFAAIA9XCFUxTDfZlr6PM+VpaoAgE+5Tqia57nv+1j8M45YUgEA+JiLDHambWqSn67+u8bzBTIypwp2Ulq3e511qmKPmvTfb7YGACjPRUJVkKUAgG+5wpwqAICvE6oAADK41PDfK57vVFPUfDoAIKPiQpXYBADswfAfAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABsWtqG6bGgBgD8WFKrEJANiD4T8AgAyEKgCADIQqAIAMhCoAgAyEKgCADIQqAIAMhCoAgAyEKgCADIQqAIAMiltR3TY1AMAeigtVYhMAsAfDfwAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABkUt6K6bWoAgD0UF6rEJgBgD4b/AAAyEKoAADIQqgAAMjhuqJqmqW3buq7rum7bdpqmJ2fWd9q2/VxbAYDiHTRUTdN0u93meW6apmmaeZ5vt1vf999uFwDAYwcNVbfbraqqcRynaZqmKS7ZG4bh4clRxBrHcVl5UtkCAMjuoKEqrIfwmqap/s1PG3HQeB8A8EUHXaeqaZp3Q1LUtKqqattWwAIAPqw+xWKYMcWq+mHpzoeLpP905imeL/BJdf3Psvzn262ACyqt2z308F/o+z5NsXpyWtd1MZsqBgp/KlbdXyf4oh2eGQBwHYeOkKlAVVXVOI4/5aRYfGF9JDLQ/VMrLTIDr1Cpgp2U1u0et1KVClRRgnoyTer+piez2gEA9nDQiep93w/D0DTNK8HovlIFAPBhB63L/TR+dy+GCDfxy/Af8DrDf7CT0rrdI1aqUjy6rz/1fd+27bqOFefM8xw3xTlVVXVd97EGAwAcOlTN8/zK+cuy1HU9DENacr3rOnvaAACfdJ26XCz++Xzlz9LqkMArDP/BTkrrdgt7toX9dIFXCFWwk9K63eMuqQAAcCJCFQBABkecqL6r5xvOFFWlBAAyKi5UiU0AwB4M/wEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGRQ3IrqtqkBAPZQXKgSmwCAPRj+AwDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyKC4FdVtUwMA7KG4UCU2AQB7MPwHAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQQXErqtumBgDYQ3GhSmwCAPZg+A8AIAOhCgAgA6EKACCD64SqaZratq3ruq7rtm2nafp2iwCAglwkVE3TdLvd5nlumqZpmnmeb7db3/ffbhcAUIr6GlfDxUIJ4zi2bbs+snl2dX2R5wtkVNf/LMt/vt0KuKDSut2LVKpCSlRVVTVNU1WVQUAA4DMusk5V0zTrRAUA8GHXrMvFFKvq0fDfr+/zki8UnEVd/7Pn3f9vVf3fnvdfGV6kTKUN/12kUrXW9/0wDFVVjeN4f2tRP124jF1DSV0LPUAGlwpVqUBV/fekdQCAvV1nonrf95Gouq5blkWiAgA+6SKVqhjya5rG5X4AwFdcZAbZw1WpHp52jecLZGSdKthJad3uFSpVqTp1P+TX971xQADgAy4VquZ5/mpDAIBylVWXK60OCbzC8B/spLRu9zpX/wEAfJFQBQCQwRXmVL3l+U41RVUpAYCMigtVYhMAsAfDfwAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABkUt6K6bWoAgD0UF6rEJgBgD4b/AAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMihuRXXb1AAAeyguVIlNAMAeDP8BAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkUNyK6rapAQD2UFyoEpsAgD0Y/gMAyECoAgDIQKgCAMjgBKGqrutpmp6cME1Tfadt2w+1DwDg+KGq7/uqqp6HKgCArzvu1X9930/TNM/zH8+MyDWOo+oUAPAtBw1V0zQNw/D6yVVVSVQAwBcdNFS1bRsLSvV9/2K6mqYppSsBCwD4sIOGqrfEEOHtdov/RgizyCcA8ElHn6j+uq7rlmVZlqVpmurn0cD76wRf9NEnAwCczRUqVZsp6rHCwk8z3FWwAIA9XKFSdV+UimKVhRgAgI+5QqgSngCArzt9qJqm6Xa7bYpVMfbnGkAA4GNOGar6vk8b0cS/8zzHYqHTNMWRruu+2UQAoDBXmKi+LEtd18MwpBWtuq6L/W0AAD6jvszVcKlM9WTUr66v83yBXOr6n2X5z7dbARdUWrdb2LMt7KcLvEKogp2U1u2eck4VAMDRCFUAABlcYaL6W55vOFNUlRIAyKi4UCU2AQB7MPwHAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQQXErqtumBgDYQ3GhSmwCAPZg+A8AIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACCD4lZUt00NALCH4kKV2AQA7MHwHwBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGxS3+aUV1AGAPxYUqsQkA2IPhPwCADIQqAIAMhCoAgAyEKgCADK4cquq6nqbp260AAIpw2VDV931VVUIVAPAZF1xSoe/7aZrmef52QwCAglwtVE3TNAzDt1sBABTnasN/bdsuy7IsS9d1324LAFCQq1Wq/uj5NjVPWIodAHiiuFAlGwEAe7ja8B8AwFcIVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZ1EUtMVDXZT1f4BV1/c+y/OfbrYALKq3bVakCAMhAqAIAyKC4FdWfb1NTVJUSAMiouFAlNgEAezD8BwCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkEFxK6rbpgYA2ENxoUpsAgD2YPgPACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgg+JWVLdNDQCwh+JCldgEAOzB8B8AQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGxa2obpsaAGAPxYUqsQkA2IPhPwCADIQqAIAMhCoAgAwOHaratq3ruq7rtm2fnDZNU33n+bcAAOR13Inq68v05nmu69occwDgsA5aqYo6U9d1y7Isy9J1XVVVfd8/PHmapqqqxnFcVuIgAMBnHLT8E2WqddvujyRt287z/MoTUe4C7tX1P8vyn2+3Ai6otG73iMN/UWRqmubd74pvbNvWhCoA4MOOGKoeappmnueHN8Xx2+0W/x2GobLIJwDwWQedU1X9O61q48lMqTQBK0pcPxWr7q8TfFGWJwUAXNVxK1WvzzQfx3EdoWKFhZ/KWipYAMAejlupeuhh/en+YBSrXAAIAHzMEUNVhKRNqemnylMlPAEAB3DEUFXdXfr35HrAaZput9umWBUJzDWAAMDHHDRUxTqfkYoiNlWrxT/7vk8b0aSyVt/3sapCWjj0880GAIp10Inqbdt2XTcMQ7rsbjMbfW1Zlrquh2GIxRSqquq67qfl1wEA9nD0pU7Tep6vnBllqicnl7a0K/AKK6rDTkrrdgt7toX9dIFXCFWwk9K63YPOqQIAOBehCgAgg4NOVN/P8w1niqpSAgAZFReqxCYAYA+G/wAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIobkV129QAAHsoLlSJTQDAHgz/AQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZFDciuq2qQEA9lBcqBKbAIA9GP4DAMhAqAIAyECoAgDIQKgCAMhAqAIAyECoAgDIQKg6k+eLbB3cqRtfaf+3nbr9p258pf1fderGF0ioAgDIoLjFP62oDgDsobhQJTYBAHsw/AcAkIFQBQCQgVD1wItXW7xyWsa7et3nG5ax/Z9v2Odf/BdP8+Z5lzfPHnf1Om+et07z5rmkS4Wqtm3ruq7rum3bb7cFACjLdSaqr1PzPM91XZuTDgB8zEUqVVGa6rpuWZZlWbquq6qq7/vvtgoAKMdFyjlRplo/l/sjcfCV55vxtM8/4mEb5qXQsMM2rK7/WZb/fPIR9z5Nw47fMC/FJV2hUjVNU1VVTdN8uyEAQLmuEKoekrEAgE+6zkT1h1f8TdO0OX7Y62ZLaJiX4t3TNOzd0359V3X9vx9+xL1P07B3TyvhEV88zaoKv3adUBWDgM8VNbILAHzSZYf/ggWrAIDPuMi0/Bev/gMA2MlFKlWbaemuBwQAPuwioSrW+YzBvmmabrdbZfFPAOCDLhKq2rbtui52p4lENY7jekLVZbYFrOv6lSn5xxEXYKYX/1yNr+7a/+3m/JW6rs/1l8Y0TfWdc/0U+r6P98+53vwPX/lwrrfQeX95T/rJ81MPdcbn8kvLtYzjOI7j5uBlnnVsv5N24zm+cRzjBW+aJo3Gnr39327UL6Un8u2GvCG9/msnegrpPZO+uP90OqaHr/y5fn9P/cu7bvyJuq2feqjLdMGvuPJzC/Gm3GwLeJbPhaTrujOGkvuO5Fy/UZvWxpvnRJ16Ei0/XeOj2WcJIhubj5rUTX61UX/lXO3ffPjE63+W9//DT54jf/I/6aGu0QW/7jS/Ib92/0Fwro+G5e6vxhO9He9f6vgFO0s3ef8pfLo3z/Lf75+zdCoh3i3fbsUv3b9Vuq470S/vxrl+c5czf/I/zH9HbvzzHuq8P4jfuc7inw9d4zLAtm2XZamqqu/7YRi+3Zw3NE1z6hH0zYfFSd9OaZZhfHE60zTFK9+27bneTpu3yrlmI61N0zTPc9d153r91070y5ve7euDTdPM8/yV9vzRH3uoU7zsuVxkovpbivoBf9c0TeuOJD6aq/Msypp68b7v+74/41Wl0f4nU2SOLN4tt9ttGIZhGG6321neOWmubpqfe7op3mtnfOfHez5e9rZt4ymc4lqBdBn7+mD8Lpyi/WsPG3ztLriIUPXTtoCfbkfZUig5YwcfnXp1ttpb3/dnLzBUq9GE+GP9RM9lnud5npumiXkkwzCcK5eE+KhM0/LOom3b6LyHYYhEcpa+PN7h8zyvo/n3mrOXq3bBRYSqq/7wziIuz45Qslnq4iyWZRnH8Vyd+jRNwzA0TXPGjjzEDJ7U/vhFPuwgyEPLskS9NsZHzjV8H85Ypqqqqq7r+IsiEnmsuXOWX974y/N2u0WN81zv+cIVEaoeOstv19mlAlV8urO3DCQAAAJKSURBVJ33ZU/rDJ3lAy71gjGIGf+NfuUsf2bcv1ui2HCW9m9KI+dqfDjRVKS11Oz0W9D3ffxRdIrXP2Ypdf9a/r1Q4Lyfn/eu9FzWLj5RPdVR1wfP0ileQMxbbJrmFB9kG9H4TWktfS6f5RPh/t1+ovafpZ33TtrsexFKTlemejjX+3TWL/tJu60Cu+DrV6psC/hFMdhxxkRV/ftxsOlOTjTRfpqm9YW+6+u0T9FHxn5Tm5f6RK9/9UNfcpbGhzO2uXo61/sUzyVG/dJ/T91tPWz2SZ/LK64fqmwL+C3rWZYbp4hZ6W+svu/jqv44cuGPg0P56fU/y4zpaGe820/65jlvX372X9402BcT8k506eK9TRccYfHKXfDeC2EdweZT+ETr122cay3aJ53fWX4E9xcqnmvxzLVzrSidbF7/s7z5w+ZX4HQv/rk+cO7dz2n7dovesGn8KT4zf3rDXKYLfkW9PNqX55KuMcrO5137wubjS5WGk77+Pnm+67yv/9nf+Rvn/UG8paBQBQCwn+vPqQIA+AChCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIIP/D/OfKXJaxQ9OAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h.Draw()\n",
    "c.Draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing Z-peak from simulation with data\n",
    "\n",
    "Let's do a concrete exercise with the dataset we have at hand.\n",
    "\n",
    "The cross-section for the simulation dataset (`DYToMuMu_M-20_CT10_8TeV-powheg-pythia6`) is **1870** pb, since the Z-boson is set to decay to muons during the event generation, the branching ratio has been incoporated in the cross section value.\n",
    "\n",
    "As for the data, the full luminosity of 2012 data-taking runs is around 20 fb$^{-1}$, but we are only processing a fraction of the data. If we sum up the luminosity collected in the released data, that corresponds to a luminosity of 4.42 fb$^{-1}$ = **4420** pb$^{-1}$\n",
    "\n",
    "Therefore, the expected total number of Z-to-mumu events in the data is \n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "   N = 1870 * 4420 = 8265400\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "[1] Open data record for MC dataset : http://opendata.cern.ch/record/7741\n",
    "\n",
    "[2] Open data record for parked data: http://opendata.cern.ch/record/7741"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Getting the data-files as before\n",
    "dy_mc      = r.TFile('dy_mc.root')\n",
    "data       = r.TFile(\"Run2012B_DoubleMuParked_nanoaod.root\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "239370L"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Let's check the number of events in each file\n",
    "## The simulation dataset is small enough to be processed in this notebook\n",
    "dy_mc.Get('aod2nanoaod/Events').GetEntries()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For the datafile, we are processing only a fraction of all the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16922L"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.Get('aod2nanoaod/Events').GetEntries()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How much luminosity do we have in this datafile?\n",
    "\n",
    "From (http://opendata.cern.ch/record/7741), the full luminosity 4.42 fb$^{-1}$ corresponds to 29,308,627 events.\n",
    "\n",
    "Therefore, by simple proportion, 16,922 corresponds to 0.002551 fb$^{-1}$ or **2.551** pb$^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Effective luminosity of in this file is 2.5520 pb \n",
      "\n",
      "Expected number of Z-to-mu-mu events in this file is 4773.2368 \n",
      "\n",
      "We should weight each event in the simulation by a factor of 0.0955\n"
     ]
    }
   ],
   "source": [
    "N_mc     = 50000   #For the sake of processing time, let's process only 50k MC events\n",
    "xSection = 1870.4\n",
    "\n",
    "fulldata = 29308627\n",
    "tinydata = data.Get('aod2nanoaod/Events').GetEntries()\n",
    "\n",
    "lumi = (4.42 *1000)*tinydata/fulldata\n",
    "print('Effective luminosity of in this file is %.4f pb \\n'%lumi)\n",
    "\n",
    "print('Expected number of Z-to-mu-mu events in this file is %.4f \\n'%(lumi* xSection))\n",
    "\n",
    "weight = lumi * xSection / (N_mc) \n",
    "print('We should weight each event in the simulation by a factor of %.4f'%weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the event-loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Same functions as before\n",
    "def invmass(pt,phi,eta,mass):\n",
    "    m = 0\n",
    "    px = pt[0]*np.cos(phi[0])+pt[1]*np.cos(phi[1])\n",
    "    py = pt[0]*np.sin(phi[0])+pt[1]*np.sin(phi[1])\n",
    "    theta0 = 2.*np.arctan(np.exp(-eta[0]))\n",
    "    theta1 = 2.*np.arctan(np.exp(-eta[1]))\n",
    "    pz = pt[0]/np.tan(theta0) + pt[1]/np.tan(theta1)\n",
    "    e0 = np.sqrt(pt[0]**2 + (pt[0]/np.tan(theta0))**2 + mass[0]**2)\n",
    "    e1 = np.sqrt(pt[1]**2 + (pt[1]/np.tan(theta1))**2 + mass[1]**2)\n",
    "    \n",
    "    m = np.sqrt((e0+e1)**2 - px**2 - py**2 - pz**2)\n",
    "    return m\n",
    "\n",
    "def deltaPhi(phi1, phi2):\n",
    "    x = np.abs(phi1 - phi2)\n",
    "    if x >=  math.pi: \n",
    "        x = 2.*math.pi - x\n",
    "    return x\n",
    "\n",
    "def selectMuons(muons):\n",
    "    select = False\n",
    "    if (muons['n']==2 and np.all(muons['isTight'])\n",
    "        and muons['charge'][0]*muons['charge'][1]<0\n",
    "        and muons['pt'][0]>25 and muons['pt'][1]>10\n",
    "        and abs(muons['eta'][0])<2.1 and abs(muons['eta'][1])<2.4\n",
    "        and muons['iso'][0]<0.12 and muons['eta'][1]<0.2):\n",
    "            select = True\n",
    "    return select"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the histograms to contain the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_dimuon_mass_mc   = r.TH1F('DiMuon_mass_mc','DiMuon_mass_mc',50,80,100)\n",
    "h_dimuon_mass_data = r.TH1F('DiMuon_mass_data','DiMuon_mass_data',50,80,100)\n",
    "\n",
    "h_dimuon_mass_mc.Sumw2()\n",
    "h_dimuon_mass_data.Sumw2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_dimuon_mass_data.Reset()\n",
    "h_dimuon_mass_mc.Reset()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function to fill the histogram\n",
    "\n",
    "We want to process the simulation and data in the exact same way, in order for them to match.\n",
    "\n",
    "Therefore, we can define a function, taking a `TTree` as input, with `weight` and number of events to process `N` as arguments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FillHistogram(t,h_dimuon_mass,weight,N=-1):\n",
    "    if N==-1:\n",
    "        #Process all the events by default\n",
    "        nEntries = t.GetEntries()\n",
    "    else:\n",
    "        nEntries = N\n",
    "    print('Processing %s events in the tree'%nEntries)\n",
    "    for i in range(0,nEntries):\n",
    "        t.GetEntry(i)\n",
    "        nMuon   = t.nMuon\n",
    "        pt      = np.array(t.Muon_pt)          # note: this is an array\n",
    "        phi     = np.array(t.Muon_phi)         # note: this is an array\n",
    "        eta     = np.array(t.Muon_eta)         # note: this is an array \n",
    "        mass    = np.array(t.Muon_mass)        # note: this is an array\n",
    "        isTight = np.array(t.Muon_tightId)     # note: this is an array\n",
    "        charge  = np.array(t.Muon_charge)      # note: this is an array\n",
    "        iso     = np.array(t.Muon_pfRelIso04_all)\n",
    "        muons   = {\n",
    "            'n':nMuon,\n",
    "            'pt':pt,\n",
    "            'phi':phi,\n",
    "            'eta':eta,\n",
    "            'isTight':isTight,\n",
    "            'charge':charge,\n",
    "            'iso':iso,\n",
    "        }\n",
    "        if selectMuons(muons):\n",
    "            h_dimuon_mass.Fill(invmass(pt,phi,eta,mass),weight)\n",
    "\n",
    "    h_dimuon_mass.GetXaxis().SetTitle('m_{#mu#mu}[GeV]')\n",
    "    return "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the tree from the file\n",
    "t_mc = dy_mc.Get('aod2nanoaod/Events')\n",
    "t_data = data.Get('aod2nanoaod/Events')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 50000 events in the tree\n"
     ]
    }
   ],
   "source": [
    "# For simulation, we want to fill the histogram with a weight we have computed before\n",
    "FillHistogram(t_mc,h_dimuon_mass_mc,weight,N_mc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 16922 events in the tree\n"
     ]
    }
   ],
   "source": [
    "# For data, we want to set the weight to 1.0\n",
    "FillHistogram(t_data,h_dimuon_mass_data,1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finally, let's draw both histogram together for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<cppyy.gbl.TLegend object at 0x7f808a21f510>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dXdKjyJk2YPjCu7LDoPCpvY3qJfScOEAzq3Bto/twHIDDsy6+g5zKoQAhXikl8XNdUdHxNkJSCoF0KzN5yPu+zwAAeM7/+3QDAACOQKgCAEhAqAIASECoAgBIQKgCAEhAqAIASECoAgBIQKiCBMqyzG8oy7Ku69l71XUdVhgtCZafsW3buGbbtsleCW8UdptbuwewO0IVvFbXddfr9YHos7y+IEVd1wIZbMofPt0AOJSmaUZL6rruui7Lssvl0jTNsF8qmC6Jd1xITtfr9YlmcgRhH5CrYDuEKkhpmpBCMArDeZfLZXhhqOWehhDFZummAtggw3/wDrEH60v9CrfCU1heFMWTrQIgIaEK3iH2YA2H7dq2vTWNvaqq7HYCWxj3Kctydjxx4bmGs+NvzZsePuxwVv6tscuVYpOG8+6H88+GbSvLcjZlhpd29yXMrrnyAZ/vGhw95kK2nj779IWP3o7Ru7B+gwCJ9cDTYqfRwjohJw3XCUuKophd59YDxk6vuE7TNPHW6ZLpvaYL734yxIdN+zES7h5f9VDTNLNdcaOXttBdN3qu2WfJft7+61f7kluNDMurqoprLrwj03d59tb1GwRITk8VvEnsTljZ7RG+Had9DAnH/i6XS/gjfiXHSDFb0yGsP135yY6Q6/Uag0VMFZfLpeu66fLhc7VtG2aeDXPJrVaF7r1phO26bviOrFxtvXimQlEUYdPFvDidNjd9R4YbJK7WT1Jy2LtmN0hc88luReC+53MZsKanqh90zIT/Xe6pmt46fJDwpTv8e3rr0LSnatpzNlpztmtk9oU/3IsTHvNWX9EwGcw2+Naz3+oUnLZ/9h25u9qXzL7GfrDPTIPjrQe5u3D9BgFeQU8VbFToaBl1ZsTOkud7HUKXzHS0qyzLW51ksytPG/lVoyeKL+3W8uEdq6p6pp+s7/umae4+wsrVpoaTw0Y3TZeUZVlV1cII4N2uslsbRB8VvIdQBds1DTfha/XWvJ8HzH7d3opKL/pufvhhwzz36TTtafgYTuueTvr+6mrrLYTg2QecvpzsK6OrtzaIierwHkIVvMkDnUxxokxcErqXEoabhVC1fvnHhSQRznQLk7FGK8Tut67rLpfLrcsHrVxtveUJcLeWh9P34it6oNDr3Q0CvIJQBe+2PpqMRgATjv0dpnxojB1xKxVFMRtW2rYdnlEYLx807ddZs9qLhCQUYlCc3v6ljsnZDZKwaxNYIFTBmzzW4TEcAUw49rfZPqcviacohqlIYaLowmhXGNQLE6TiZuy6brQ1Vq62xvKcs9Hy6XmCyy9n6tYGOcbbDdsnVME7DM91/9IdhyOAz4z9LRdnX7nypsRGhinkX9osYVCvH5RLeGa1VGI1h8di0DMbBEhCqIJ3iEWGvtpfFUcAnxz7uzUvZ7Y98dv9gSd6vzUzu28VFh8tWbnaevGO6x/zVmWy9aYbZBcpGQ5AqIJXCaehhQuthCULZ8svCOEmxLKVHV2jL+aFK+F0XbeyosFmjRJD7BQcCttwmixHr3Hlal8SH3PYzrZtb8XcaQAalv1cY7pBHpjqDjzirVWx4KDWdOpMC3IuF/+cLpw+yHRhXDlMT66qKhaEnD5yXBhXjk80KnR56+PiyaqSs/e9VQNzunzY2jD/afRGTOtqLr/Ylas99jLjYw43ezZXejTMqZq+nOFcq35QVjQuXL9BgFcQqiCB5VB165tsZajqbxc0DwvvXg4vPMXsI8y2fNrg2fsuNHil2fuuD1WzPX+h8dPXcqubcPRiV672JbfaGbbe8JFvvXfD5aOtF8QL4KzfIEByeT93bU5g19ofyh+W1x8Ob+1l1C+K1T6Hr3R2YbZ6y3x1A65vZ3jMhY08LF46fOpwl1Fj6rqePuCXNgiQkFAFAJCAieoAAAn84dMNAA7igUpRr2rKy3z1OnpqGcCp6KkCPsC0HuB4zKkCAEhATxUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACf/h0A94qz/NPNwEATqTv+0834X3OFaqyk727ANuX57lP5qM6W1+G4T8AgASEKgCABIQqAIAEhCoAgARON1F9edKcyZIAwGNOF6rEJgDgFQz/AQAkIFQBACQgVAEAJCBUAQAkIFQBACQgVAEAJCBUAQAkcLo6VYp/AgCvcLpQJTYBAK9g+A8AIAGhCgAgAaEKACABoQoAIAGhCjiyPP/+6SYAZyFUAQAkIFQBcFLtwGh5vHW6/tual9yuG78L+anqNuX5uV4vkOff+/7bp1vBkg9+Mo/KQTdNU5ZlXdchPIVbh22bLtmXW5v6n//8Z9d1//mf//m2Zzyq0/VU5Ys+3ToA3qppmr7v+76vqupyuWRZFkJVXCH+Xdf1B9q3wmP9T8N7dV3Xdd2TD0h2wlDVL/p06wD4jJCZ2rat67osy9Hy7Oeoked5HCKMP8jj7/Nwl7Is7/5oD3efCg8+fISwflwSWhjufrlc4l1mX9foQYb3quv6n//853/913/961//KopidNMXNyHnC1UAcMsomlRVFbtwuq6rqurWHcuyLIqi7/umaa7Xa1wef7EvdP+E3rLwdHHltm27rot3D/1ncUnXdSH/hSetqupWBrper+Hxi6IISy6XS3iiqqqu1+tf/vKXv//973/+85+7rosPOHoVrLTLa//leT7cgdq2DX22QyFxv7lhADzvK4Uwfsmyf6xc9YHZdWVZXq/XNfPTu65rmibcJfvRvxVzzN1nGf4Rvr/quu77fjgWGW4N34AhaYXmhXlgs6Eq9niFJoXvyr7vw+NPY1N4pbM3scb+QlXowBSYAI7qK+kn8VkIZVlOO6tCXlnopnqF0F8Qnj12HMSYled50zRhpDKssL4rIc/zoijKsgw9VcObyrIMHXLTm1hjZ8N/C0k8TjYMpC4A1ojn+oWoMbo15IwwNDa6Vzb4hV8URZyVlQ06n55pVZZlwwleoZGx7yr8b2h80zTDmebDxg8fKi4PPVvTL8rQ32Y21eOWJ25vSuhZDcKgbxD6V9c8wr5eL/C8LPvHp5vAHR/8ZB5+G8avlaqqwt/ZYLZTaGTorIp/DL9Jh99QYXiuKIr4mNnkl38Q7jVaJ9xx+IDxWe4+6dS0qaNByaqq/vu//zvLsj//+c+jlYdftQ9v4ScfYV/2VEAiDPw1TTPq5ww/I5qmib8Pbv1EOFvBDECdqu3b7ydz27ajr5vpklRPMfwj+7knbM2TjtZZvkvCV7HfN/cxu3m1MTllWTYKVbOnqs6+rrO9u4BQtX0n+WSenf0dyo2meorpt+HHz9k6yZsb7WNOVV3XYercws4Xez5Dx+ZCZ9VjXvPKADiFcDbfSNpurenjm178Zjs4+69t2+v1GucATo2Sfpi7NztlL9vz5QUAgC3bQU9VzFJhslT4367r4omv06QfOqskdADgbXbQUxVMe55CPdkQrdL2oAIAfNX+ZpCFemhx8t3of4P8xoXEzzZjDjBRfft8Mh/Y2d7cHQz/LQt9VKEsW6iBFpa8ufQtAHByuxn+W9D3fZ7n1+s1nq26cGlJAIBX2F+oKsty2pcYThwN3VTmVwEA73euwc6zDe4C5lRtn0/mAzvbm7u/nqonLZfxPNV7DwAkdLpQJTYBAK+w+7P/AAC2QKgCAEhAqAIASECoAgBIQKgCAEhAqAIASECoAgBIQKgCzi7Pv3+6CcARnK74p4rqAMArnC5UiU0AwCsY/gMASECoAgBIQKgCAEhAqAIASECoAgBIQKgCAEhAqAIASOB0daoU/wQAXuF0oUpsAgBewfAfAEACQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACQhUAQAKnK/6pojoA8AqnC1ViEwDwCob/AAAS2GWoyvO8ruvRwrIs8zzP87wsyw+0CQA4t/0N/4VJUW3bThcGXdfleW6YDwB4p531VE07qLIsC11TVVX1fd/3fVVVt9YEAHiRPfXotG17uVzC30VRxM6q0E01fCHTJXH5jl4v8Lw8/973355fh9fxyXxgZ3tz99RTFRJV0zTTm4qieHtzAAD+z25CVRjjmyaq0eSqQMYCAN5sHxPV67ruuq6qqrIsZ1PUrLZtp2cCLhf/XHCqDkzYozz/Prf4l+lBb7CPYPiFMvy+CF8f4dbR8tESGNpBqGrb9nq9FkWRZO65bARHNRuV8lyE4qY4TzdomqYsy7qu27aNs3iH3xrTJTC0g+G/mKXKsgy7e5ZlXdfd7bXyYwKAZU3TxDPHQ2YKoSquEP/e7Enl6wdweLUdhKqgG4hL4gBfXBhven8LAdivkJnatq3revibPGapYXbJ8zz8b9u2cVZJ/kO4SyxJHcw+abj7VBxnHN19VOY63P1yucS7TA0fdvbuw2fkSTsIVW3b9gNhrnpRFH3fhx13dlq6ueoAfNUoW1RVNfwlH+ogzirLMnwxNU1zvV7j8vDNNX3kodBblv0ouBhWbtu267p499B/FpeEboW6rsOTVlW10JEW1imKItw9vqjL5RKeMdUEG3YQqu6KvwmywY8G+wfATuX5+n/5+pUfEL9Z7n6ndF03/DJa+M1/61niH6EQY1mWoe9gdGvoCev7PkyJCZNhlufDxMaH9sQXlQ365/RUJXGEUFWWZcjdoRc0+zHZ8NPtAuARfb/+X79+5TWm3x2hE6ht24VuqleY7SMI/UzhphD1wujN5XJ58ltPqEpif6EqhPfR2x9ie9M0oRNVogJgpdBPk+d5URTTr4/QGxQ7oob3ygZZJI6gpSq7EHuS4kOFRsZ59DFatW3bNM1XJxMPe9TiSWA8aQclFdaTpQD4ktmrnw3NfrNUVXW9Xodzp+q6DhPGw63Pfx+VZXm9XuMU9ev1GnJPXBLy1uVyCc14oCOtaZp499mrlfBV57ooz93Kn6faGnAGrv23ffu9PNy0xPRs0ekkTzH8I5srVfr8U7zCft/cx5zs1Z7s3QWEqu07ySdzXdfDnq0g7QzgacfBre63tznJmxud7NWe7N0FhKrt88l8YGd7c/c3UR0AYIOEKgCABIQqAIAEhCoAgASEKgCABIQqAIAEhCoAgAQOdZmaNZaLqp+qnAYAkNDpQpXYBAC8guE/AIAEhCpgl/L8+6ebAPAToQoAIAGhCgAgAaEKACABoQoAIAGhCgAggdPVqVL8E05u7rTBX2Y/GPr+2+ubAxzH6UKV2AQnN41KeS4/AQkY/gMASECoAgBIQKgCWEUNd2DZ6eZUAbA1y6cQwV4IVQB8kvOHDuxscdnwHwBAAkIVAEACQhUAQAKnm1OlojoA8AqnC1ViEwDwCob/AAAS2E2oatu2LMs8z/M8L8tydFM+MVoHAOCl9jH817bt5XLJsqwoiizLuq7L89xAHgCwHfvoqQqJqu/7tm3btq2qKsuy2BfVtm2WZU3T9ANhIQDAe+wjVGU/+qiCuq6zLOu6LvxvyE/G+wCAD9rH8F/TNMP/DSlqGLPCwpiuBCwA4M12NjMp9FFdr9csy5qmCeFptvTU7OsyEwsOI8+/9/23JKslfChg6Gxfu7sZ/guu12tIVEVRjLqjqqoKs6lCD9atzqrpeYIrvfy1AQB7trNQ1fd90zRFUXRdF2NTmKIeOrGyH4ODccbV9BEe84ZXBwDs185CVZZlZVmOYtO0Uyp0VjkBEAB4mx2Eqrqu8zwfJaRhbBKeAICP20GoCh1RcXQvCN1UodfqcrmMOqvirW9qIgBwersJVV3X1XUd6iaEJcMJ6dNbQ4FQAID32Me5jvEyNVFRFMNRv9HZeVVVjXq24mq7eL3AXUoqwPad7Wt3H8U/y7IcXnlmOq4Xr2Cj8icA8BH7CFXBcloSpwCAD9pTqEpiuYznqXopAYCETheqxCYA4BV2cPYfAMD2CVUAAAkIVcCRKYIAvI1QBQCQgFAFAJCAUAUAkIBQBQCQgFAFAJDA6Yp/qqgOALzC6UKV2AQAvILhPwCABIQqAIAEhCoAgASEKgCABIQqAIAEhCoAgASEKgCABE5Xp0rxT9ijPP8+WfbL9Gju+29vaQ7AjNOFKrEJ9mialvJchAK2xfAfAEACQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkIVAEACQhUAQAKnK/6pojqwhhruwFedLlSJTcDIbDBSwx34KsN/AAAJ7CZUtW1blmWe53mel2U5XWH5VgCAl9rH8F/btpfLJcuyoiiyLOu6Ls/z4UDecKbU9FYAgFfbR09VSFR937dt27ZtVVVZlsUeqfBHVVV93/d9H26t6/pDjQUAzmgfPTp5nhdF0bbtcEn2Y9b58O/praPH2cXrBe7K8+9vnjb+/meEvTvb1+4+hv+aphn+b0hXYSgwGP4NAPB++whVcaQvDOpdr9f497D7KiqKouu6NzUOAGAvoSoKcSrLsqIo7p7lF04YHC1cLv654FQdmMDDVo4SGkyE49nHRPWo7/umaUJH1GOlE/pHpX4pAMCh7CxUZVlWlmUY8rs7wKdgFQDwNjsIVXVd53k+mjsVZqbHAb5RwDKhCgB4sx2EqhCbRnWnQmwKN82e+ud8QADgnXYTqrquq+s6FP8cZamQt8LCtm3DVHTFPwGAd9pHVa54mZpoVAu0rut4YmCWZU3TzE6oOlsVMjiwzRb/dPYfRGf72t3Tq40p6tYM9LDCwvz0s727cGBCFWzf2b5291Sn6u7ZfE73AwA+ZU+hKonl4p+nCtQAQEKnC1ViEwDwCjs4+w8AYPuEKgCABIQqAIAEhCoAgASEKgCABIQqAIAEhCoAgASEKgCABE5X/FNFdQDgFU4XqsQmAOAVDP8Bb5Ln3z/dBIAXEqoAABIQqgAAEhCqAAASEKoAABIQqgAAEhCqAAASOF2dKsU/4agWD+6fONCBVzhdqBKb4MDWHN/rsxfAlxj+AwBIQKgCtuUNhdfzfNW/j1ODHvbldMN/wMmtnAKwhVAF7IueKgCABIQqAIAEhCpgl/r+26ebAPAToQoAIAGhCgAggdOd/aeiOgDwCqcLVWIT7JECB8D2nS5UATvlBxGwcbuZU9W2bVmWeZ7neV6WZdu2w5vyibIsP9ZWAOB89hGq2ra9XC5d1xVFURRF13WXy6Wu60+3C/hf57ygyo2L23zb4BVvgDfYR6i6XC5ZljVN07Zt27ZhXtT1eg23hl6rpmn6gWFXFsCL9P34X5Z9nywBTmEfoSoYjugVRZH9iFPhv8b7gJdSbhRYto+J6kVR3M1MoRMry7KyLAUsAODN9hGqRmN5bdt2XZf96J0Kf4chwuzHsKDSCQDAO+0jVA3VdR1iU9M0w+VVVYWp62VZdl03OkMwWi7+uUBKA6bmZuj/MvmY+bZuNSOMsG97ClXhHMDwd9M0cYxv+Hf2o8JC6L6ako2AhKYxKM/HC6dLbi0Edm03E9Xrug6Jqqqqvu+HKWo6g2o4jR0A4A320VMVhvyKopjNSaEu6LvbBAAwsI+eqjCJ6laiulwuo1A1nMYOAPAGO+ipillqGpLquo4nAMa/w3T1qqre1kKAD8rz76ZnwRbsKVQtzD3P8/x6vcYa6/FMQACA99hBqKrr+m5CCtelCZOrjPoBAO+3g1C1kjgFAHzQPiaqAwBs3HF6qlZarqiuNCgQPXr9BeCkTheqxCZgDR8VwFcZ/gMASECoAgBIQKgCAEhAqAIASECoAgBIQKgCAEhAqAIASOB0daoU/4T3yPPvk2W/zB5/ff/t9c0BeLnThSqxCd5jGpXyXH4CjszwHwBAAkIVAEACQhUAQAJCFQBAAqebqA7wfnOnPX6bLnQiDeyaniqA1+r7mX9Z9n2y5LXmilwAKQlVAAAJnG74T/FP2JTVNUJPUeBqsjWUS4U9OV2oEptgU1bWCF38NXQcoxeuXCrsi+E/AIAEhCoAgASEKgCABIQqAIAEhCoAgASEKgCABIQqAIAEhCoAgAROV/xTRXV40o1LyM3U/la4EjiV04UqsQmeNBuV1P4GMPwHAJDAbkJV27ZlWeZ5nud5WZZt245WGN76gfYBAOe2j+G/tm0vl0uWZUVRZFnWdd3lcqmqqq7rsMJwplTXdXmeG+YDAN5pHz1VIVE1TdO2bdu2ITBdr9dwa+iaqqqq7/u+76uqyrIs5i0AgDfYR6gKhuN6ocsqDAJ2XZcNUlT4I0YuAIA32MfwX1EUyzOlQsYCAPiUfYSq0bT0tm1D79TsjPUsy4qiCCsAALzHPkLVUF3XYWivaZrlNcMJg6OFy8U/F5j5zsHk+XeVpfZorviqyquwCXsKVfEcwCzLmqZ5rHSCbASbMvcz59ujv302YWWaeTj0TO+o8ipsxG4mqtd1HRJVOMtvTaJSsAp2oe9/+pdl30dLfiwH2LR99FSFIb+iKGZrfmY/TgCMTKgCAN5sHz1VYRLV7Jz07Mapf84HhM/K8/G/MK43WQhwEDsIVTFLlRPhplCYKnRZtW0bpqIr/gkfZ1wPOJUdDP/FUHVrUK8sy6qqrtdrPLPv4WnsAACP2UFPVV3X/Q0xOYV1mqZpmmblNHYAgIR20FO1niwFAHzKDnqqAAC271A9VWssV1RXGhQAeMzpeqpuTc8KPt062Jy5i6IAMON0oQoA4BWEKgCABIQqAIAEhCoAgASEKgCABIQqAIAEhCrgTfr+26ebAPBCin/+RKkq4ABuVBf7Zfr5J+lCQqcLVWITcHizUSnPRSh4LcN/AAAJCFUAAAkIVQAACQhVAAAJCFUAAAkIVcC2OEMtsilgX4QqAIAETlenSvFPAOAVTheqxCZgsxZ/9P0fH2OwTacLVQDbtDIqrQxewPuZUwUAkIBQBUdz42K6kIx9DGYJVQAACQhVAAAJCFVAAspUAghVAAAJCFUAAAmcrk6Viuqc09zpWr/MHg0G8gAec7pQJTZxTtOolOfy015N0vC3aT72UQfvZ/gPYE/6fvwvy75PlgAfsL9Qled527bDJW3b5hNlWX6mfQDAKe0sVNV1nWXZKFQBAHzcbuZU1XXdtm3XddObQsZqmkbvFHBC5sbBRuwjVLVte71eF27NskyiAgA+aB+hqizLcNZeXde30lXbtjFdCVgAwJvtI1QtC2OCl8sl/G9IXUonAADvdIRQFVRVFaaxl2XZdV1ZlrPz2ZeLfy6Q0jiDlTVCTeIBmDpCqBpNUQ8VFmantGeyESxSIxTgYTsrqTBrOoOqKIpM5QUA4I2OEKqEJwDg43Yfqtq2vVwuo86qMPbnHEAA4G12H6pCcuq6LlQHbds2LKmq6rMNAwBO5QgT1fu+z/P8er3GElbxTEAAgPfYWaiq63o2LfV9H7upjPqxR3n+3Ul2ALu2s1C1QJwCAD5o93OqAAC24Dg9VSstV1RXGhQAeMzpQpXYBAC8guE/AIAEhCqAU3B6KbyaUAUAkMDp5lQBT1o82QPgvIQq4Muc7wEwZfgPXivPv2/20WDEDgbPEKoAABI43fCf4p8czFzXwi+zu7mTv3iMfQxWOl2oEps4mOnXWJ77biMl+xisZPgPACABoQoAIAGhCgAgAaEKACABoQoAIAGhCgAgAaEKACCB04WqfNGnW8fLuQoHJ5Hnq/4BCSn+CSeleOOBrfycE6ogrdP1VAEAvMLpeqqABXNdF9/0ZwCsoacK+Enf//Qvy76PlhhCB5glVAEAJCBUAQAkIFQBACQgVAEAJCBUAQAkcLqSCstl05UG5T3mCrv/Mt03t1CfcwttANiF04UqsYktmCaVPBdfOJo8/55qr175UGlXg68y/AcAkMD+QlWe523bTpeXZRkuilyW5bvbBACc3s5CVV3XWZZNQ1We513Xhb+7rlueOAUAkNxuQlVd12VZXq/X6U2ha6qqqr7v+76vqir7Eb8AAN5jH6Gqbdvr9Rr7okbC8piiwh+z8QsA4EX2EarKshz2Qk0VRfHmJgEADO0jVC2YnbQuY8FInq/6B8DDjlynqm3b6ZmAD89hV+CKNVZW9cxWFNWce6j5R1tZcccuzNRkN3twdwWyY4eqWbIRL5WwqufsvdQIJa3R7mQHg2fsfvhvgYJVAMDb7D5UheQ0OjHw1nmCAAAvsvtQld2Ylm6u+tncmIEEAG9yhFAVClOFLqu2bcNUdMU/AYB3OsJE9bIsq6q6Xq/xzL6maUyoAgDeaWc9VXVd930/7YUKy5umaZqm73uJijNz9havZh+DWUfoqYpkKQDgU3bWUwUAsE2H6qlaY7miutKgB7OyvrmxDACed7pQJTadSsL65gCwzPAfAEACQhUAQAJCFQBAAkIVzHDRG8DnAF8lVAEAJHC6s/8AiCYVRr7Nlp1x2jSsoacK4KT6fvwvy77PLQRWOV1PleKfAK+QsNbuyoea8/hqKtjxvNOFKrEJ4BUS1tp9+KHSrgZfdbpQBQez7nc5AC8nVMHu6X4F2AIT1WETDEYA7J1QBQCQgOE/2Kgbk6XmywgB8HF6qngJl3dIQhkhgB0RqgAAEjjd8J/in8ewvjagCeAAvMfpQpXYdAwJywwCd62cyefzlZM7XagC4EtWRqU8X3V5ZsGLAzOnCl4rbf+Z3jhe6pkdbM15FXBsQhUAQAJCFQBAAkIVAEACQhUAQAJCFQBAAkoqwAe4fh/A8ZwuVO2lonqef19zbvPK1Vhw4zKFM/XZ5yqOPr79N7OvAU99DkB0ulC1ndjER0w/EGc/ItVnh1PxOUAS5lQBACRwkFDVtm0+UZblp9sFAJzFQUIVAMBnHWROVdu2WZY1TaN3CgD4iEOFKomKj1MrAeC0DhKqgrZtY7oSsPgUJ5gCnNNBQlXXdVmWXS6X8L/X67UoihCwAADe4CChKqiqqq7rLMvKsuy6rizLaWK+A6AAAAzOSURBVK5aLv65QIGrBXN18xTNg7NLeMivfKi0q8FXHSRUjaaohwoLoftqRDZ6hblS4z62gBkrf9j6qGaPDlJSYTqDqiiK7McEdgC2oO9X/YOdOkio2kh4unH1qM9L2LC0r3GzWwx4m9nPgTxf9e/jVn6IbfazbrMN26kjhKq2bS+Xy6izKoz9OQcQYHd0aLFTRwhVITl1XVfXdaiqEJZUVfXZhgEA53GQiep93+d5fr1er9drWBLPBAQAeIODhKosy/q+j91URv1OyMmGAHzWcUJVppA6APA5R5hTBQDwcYfqqVpjuaL6p0qDrqxIPmd+NWNhyd14O75t4aRuALbgdKFqmxXVH65IrnD5O033nTz/bvsDEJwuVHEqLogBwNsIVRzc3cBk/A6AJExUBwBIQKj6ANda+qpXb7G5C4p929olxoA3M2OSrzL8x9mZUAVAEkIVAEfmhBXeRqhil4zHAes5YYX3OF2oSlj8c2XFzvePyids2PqqpHcf7ca8qMe3mJ+VsDvJPwfeLG2h5l1/OzDrdKEqYfHPhyt2vlrChr30oZ55NGCP9v45sPdCzZv92joMZ/8BACRwup4qvmSuW3v+anejHkA/fYA3MBeKTdFTxR19/9O/LPs+WmJuE/AR08+iW//gPYQqAIAEDP+RxppOeL8XATgwoYoE1qQlUx8AODah6oBWzi5/f7/R+mnvj64G8CCF13ne6UJVwuKf657uM5XuRq8jz7+PHv/9GWV20841bH6LZdk/JgudYAiksfKzfyO/7qafnI+tQ3KnC1XJY9O9p3uw0t3qquJHO2b2VRtwm60CXmeSqx7vbtfpdTzO/gOAVeaKNcxUmVmzGockVAEAJHC64T8A2AJT449HqJqxfirim/f1tHMkNzLjEuCE9jU1npWEqnmpCi/dWOfxGYupYtwn6inkbz5LIC3t/6xdt3/Xjc+0H1Yzp+rlVs5YzPP7/wA4m7nvgm+v+4JY82XkK+mWc+X3tON6eX5/tdl1Hi4fsuvfW880fgsFV3a98TPt/6hdNz7T/o+abfwzpQcf+9p6YrUdb/wHnK2nKndJcwAOZmWtB17tdHOq3lxRHQCeMXediZnLcqxbR73i1zpdqBKbANiR6eSHNReZmF3HRKhXO9vw39csd2tlKyb0Pfn4T/rq4796/Zc++KYa/8Dja39Cu278A4+v/QntuvEPPP7W2r93hwpVf/rTn/I8z/P8T3/60xuebm64embOFgBwBscZ/hvG5//5n/852xkHAHDXpKep1/eU0EF6qkLX1K+//tr3fd/3v/76a5Zl//Ef//Hpds37eHWA3bHFAJ6XdoBlXS2rc/VuHKQ7J3RTDV/LdEn29YIZ1k+4/qYaY/1jr7+pxlh/X+tvqjG31l8o3RfXX1neb2uvd++O0FP1+++/Z1n2xz/+8dMNAQDO6wihapaMBQC803Emqs+e8ff777//9a9/HS7Z2ummp1p/U42x/rHX31RjrL+v9TfVmFvr5/kvd9dfWCd5ewiOE6r+/e9/313nVCO7AMA7HXb4Lxh1UwEAvMhBpuWvPPsPAOBFDtJTNZqW7nxAAODNDhKqQp3PMFf9999//9vf/pZtuPgnAHA8BwlVf/3rX3/99ddwdZqQqH777be//e1vZVlOVy7LMlwicPbWLVhoYdu2w1vbtv1A+xat3Lyb3f7L7a/rOqywzY2f3Wv/xnf+0e49XWHL7V9u/PaP3LsbP9rj9s+2ffDud8/P83x2Y+705STQH8tvv/3222+/ZVlWFEVRFNPXuPEt0DRNaNVs+2dvrarqc+39yXLjh+Jq723gHXfbHxfGP5qm+VBjZzyw83ywtVPDFs4enls+eJcbv/Ejt1+x8aPtH7yz7d/ywfvAzvPB1g5VVTW7Jy8fqls+kJ93qBfT/3i34tESdsd4/Ic9Mu4Bt3aID1pu//SzYFN75HLjo7DZN/i5vNz+0Oz4v/GT7gMNveFLO8+td+eDRttzdHhu/OBdbvzGj9z+XvtHy7e25/T32r/xg/dLO89Gjtyqqm79PFg+VJf/9wC2slelcisUr7l1C77a/rD7buQn15rNGz/OtvC5MPLVjV9V1aY+C3a9889+Veyl/Q80flNH7t32D1fb4MH7wPbfzsG7xz1/uCdkk0i06w+i5x1kTtUtYax32Kc6/Hv7Ru0viiL+Uty+6cbPsuxyuWRZNjomt+nuzlPXdV3X723UF8xu/5W3vl9oz2iCxah522ntyN3Gb/zIXbPxsw0fvA/sPNs5eFdu/NH6nz0WyrIMAeLWXr3cvM0eyGl8NtMlFw/4Yefk8KbRD4KwzgcaesNC+xdWfmcLF9xtfPx1vpEe7JE1O8/w42Ajv3Sj5e0/e+tGekr6xd/rt3aY7Ry8y42/tf5GGt+va/+WD96VO882D967G3/LR+508G75UF24dSOv6HlbOaoTmv5ACcv38nbeav9I/Imwl8YPD78Nfi4HyztPXBg3/nY+moPlnWflrvUps9/iC1/kmzp4Fxo/WnObR+5y+7d/8N7debZ88N7deTZ75H4pVO3iQH7e0ULV6GgZzk/cxdu50P5o+BmxnZb39xqf3Qu4H3d358m2PRXg7va/u2t91q1xpV18Fi80fnadjTQ7Wm7/9g/euzvPlg/euxt/s0euUDW1lb0qiYffznc2csGaHW6DP7OC5caHP4qB8CqKotjIsbT3z4K9tz+qfujXjRp8oo03zTY+3rTNI3dotv3bP3ijr+48m2r/Hhv/1eG//ueAPr31AP4wG5B3anbGXxSWd103XDj6389abn+WZXVdX6/Xoii2VrkuW9H4bG5rd10Xqt69qlmrrdl5tmy5/WvenY0YTh+OO8z2D95gtvHZto/coVvtn/5vtqWDN1rYebZvtvE7OnKDNYfq9g/kp3w61aU0m5GHL3OaiKfrf9Dd9m/5Lbvb+Lsrf9YDG39Tb8dy+7/07nzKqD2jNm/84F1u/NY29dRy+4c2ePD2X9/+m3pHFhq/8SN3tsrU8qG68QP5eZt4YxIKe1tVVWHII3Zch1unO2u2jU7UaKH9scHFxEZewvLGH9ry5/LyzpP9GE1beHWfstz+4kdFvs22PzapaZo4WBZv3fjBu9D47R+5/b2NP7TNg3fNzrPZg3e58Vs+cmdD1fKhuvED+XlHC1X9vRMlqp/ramzwvbzV/up2nZvtvIrljR9t83O5/+LOs7v2r3x3PmjUwtGOvfGD91bjd3Hk9vc2frSXg3d559la+5cbv9kjdzZU9fcO1Y0fyE/K+7kL8RzA3qeYbL+FC3bd+OzoO8/G29+2bZips8f23238xh2+/bveebbc+Fm7/iB62GFDFQDAOx38MjUAAO8hVAEAJHCoOlUAcHLDcmijSUth5lZY/qX5TGsmSA0db7LUSuZUAcBx5Hke/igG9WZDBdrRmlVVDYuOLlgoYFuWZdd1VVW1bRsreZ42WghVAHAcIVQ1TRO7i2KiqqoqLGzbNiyZzUkLDzvNDMPlbdteLpfZ1U5CqAKA4xiln4WgM41fC0KP1Gjl8ODDZHYre52EieoA8AF1Xdd1HcpT5XlelmWIJnVdhyUrx+buPkuWZbF8+VBYOOqpCs9eluXo2cP/zi487QyqKT1VAPABw8lP2c+XFi6KIvzv+mlPo4eNX+5f6jqKTYqGXVPTh1qz5FT0VAHAx4Qp3m3bxsvR9H3ftm3oRprOLv+SlfOlghCeRtfpG0a60MK4JDz46Co6JydUAcDHxIwSMk28NF7CMbVR7sknwvKu64qiiM9blmXsMBs2NQa1o15q5hnqVAHAMc0mnuEljWNPWEhIXdcN7xISVZj1FR8txqxw3yQTvw5DqAKAIxv2NmU/x6Dp8OJo5ZGqqq7Xa5jMnhn7mzD8BwCHFXLP7OSq6cKqqvqJYd9V+Pt6vcYTFV/S6N0SqgBg09q2jTOf8jwPgWZ24ex9sywLpapGhgtjWhquMJxxNVwtrmlC1YhQBQBHFiZRhcJX4UzDuq6n1RPCaqFcVpxHNR3gi1OyjP1NmVMFALsxG2WW800YpLv+EJc3TROvVzNcLfZgzV7EJl70xtjflOKfALBp4Wow4fs61DoP/UnThdli+c1hF9TCyN0ztRJOXvxTTxUA7Mbs2XnLp+xF4foza1b7YqP4X+ZUAcAOhOsDFkVxuVyGFzCeLsx+9F29s3mxq+zM9FQBwKaFq9aMIktd19OF0aiM5xu0bRs6zM48gV2oAoCti3Oh7i781HymL11n8KiEKgDYtNk+J2NtG+TsPwCABExUBwBIQKgCAEhAqAIASECoAgBIQKgCAEhAqAIASECoAgBIQKgCAEhAqAIASECoAgBIQKgCAEhAqAIASECoAgBIQKgCAEjg/wNG902YLOJa5AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h_dimuon_mass_data.SetStats(0)\n",
    "h_dimuon_mass_mc.SetLineColor(r.kBlue)\n",
    "h_dimuon_mass_mc.SetFillColor(0)\n",
    "h_dimuon_mass_data.Draw()\n",
    "h_dimuon_mass_mc.Draw('same hist')\n",
    "c.Draw()\n",
    "c.BuildLegend(0.7,0.7,0.9,0.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data-MC agreement is not bad with this very simple example.\n",
    "\n",
    "\n",
    "As the data sample is very tiny, fluctuation in data could make the comparison with simulation difficult."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset for H to 4l analysis\n",
    "\n",
    "With the above knowledge, you can process the simulation files in the H to 4l analysis and compare them with data.\n",
    "\n",
    "Here are the relevant data that you will need to compute the weights for respective simulation samples.\n",
    "\n",
    "Can you compute the weights below?\n",
    "\n",
    "| Sample Name  | cross-section(pb)  | Nevents  |\n",
    "|---       | ---                | --- |\n",
    "|ZZto4mu   |   7.6865000E-02    | 1499064  |\n",
    "|ZZto4e    |   7.7241300E-02    | 1499093  |\n",
    "|ZZto2e2mu |   1.7656500E-01    | 1497445  |\n",
    "|DYJets    |   2.7903452E+03    | 30458871 |\n",
    "|SMHiggs   |   1.2932700E+01    | 299973   |\n",
    "\n",
    "|Dataset | Luminosity(/pb)|\n",
    "|---| ---|\n",
    "|2012B |4410|\n",
    "|2012C |7060|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Compute the weights here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extension: To handle more than 1 histograms, you may use a more generic function like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FillHistogram_extended(t,histograms,weight,N=-1):\n",
    "    if N==-1:\n",
    "        #Process all the events by default\n",
    "        nEntries = t.GetEntries()\n",
    "    else:\n",
    "        nEntries = N\n",
    "    print('Processing %s events in the tree'%nEntries)\n",
    "    for i in range(0,nEntries):\n",
    "        t.GetEntry(i)\n",
    "        nMuon   = t.nMuon\n",
    "        pt      = np.array(t.Muon_pt)          # note: this is an array\n",
    "        phi     = np.array(t.Muon_phi)         # note: this is an array\n",
    "        eta     = np.array(t.Muon_eta)         # note: this is an array \n",
    "        mass    = np.array(t.Muon_mass)        # note: this is an array\n",
    "        isTight = np.array(t.Muon_tightId)     # note: this is an array\n",
    "        charge  = np.array(t.Muon_charge)      # note: this is an array\n",
    "        iso     = np.array(t.Muon_pfRelIso04_all)\n",
    "        muons   = {\n",
    "            'n':nMuon,\n",
    "            'pt':pt,\n",
    "            'phi':phi,\n",
    "            'eta':eta,\n",
    "            'isTight':isTight,\n",
    "            'charge':charge,\n",
    "            'iso':iso\n",
    "        }\n",
    "        quantity={}\n",
    "        if selectMuons(muons):\n",
    "            quantity['mass']=invmass(pt,phi,eta,mass)\n",
    "            quantity['pt0'] =muons['pt'][0]\n",
    "            for histogram in histograms:\n",
    "                histogram['h'].Fill(quantity[histogram['fill_var']],weight)\n",
    "            #h_dimuon_mass.Fill(invmass(pt,phi,eta,mass),weight)\n",
    "    return "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning in <TFile::Append>: Replacing existing TH1: DiMuon_mass_mc (Potential memory leak).\n",
      "Warning in <TFile::Append>: Replacing existing TH1: DiMuon_mass_data (Potential memory leak).\n"
     ]
    }
   ],
   "source": [
    "h_dimuon_mass_mc   = r.TH1F('DiMuon_mass_mc','DiMuon_mass_mc',50,80,100)\n",
    "h_dimuon_pt0_mc   = r.TH1F('DiMuon_pt_mc','DiMuon_pt_mc',50,0,200)\n",
    "\n",
    "h_dimuon_mass_data  = r.TH1F('DiMuon_mass_data','DiMuon_mass_data',50,80,100)\n",
    "h_dimuon_pt0_data   = r.TH1F('DiMuon_pt_data','DiMuon_pt_data',50,0,200)\n",
    "\n",
    "\n",
    "h_mc_list = [\n",
    "    {'h':h_dimuon_mass_mc, 'fill_var':\"mass\"},\n",
    "    {'h':h_dimuon_pt0_mc, 'fill_var':\"pt0\"},\n",
    "]\n",
    "h_data_list = [\n",
    "    {'h':h_dimuon_mass_data, 'fill_var':\"mass\"},\n",
    "    {'h':h_dimuon_pt0_data, 'fill_var':\"pt0\"},\n",
    "]\n",
    "for hist in h_mc_list+h_data_list:\n",
    "    hist['h'].Sumw2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 50000 events in the tree\n"
     ]
    }
   ],
   "source": [
    "FillHistogram_extended(t_mc,h_mc_list,weight,N_mc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 16922 events in the tree\n"
     ]
    }
   ],
   "source": [
    "FillHistogram_extended(t_data,h_data_list,1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<cppyy.gbl.TLegend object at 0x7f808b6bae00>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3df5Kjxv0/YPpbvlYMylblKOtDfFILciWHyNwkrvIKnIPx/aOzHQyIYTQtCcTzlMs1Awx6o1+8trtpQt/3BQAAn/P/nl0AAMArEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAy+OnZBcBBtW3bNM3squqH2bVN07RtW1VV+vPhrtq2XX7c+OcrN2al9BJ4SuHIhCp4jrZtu66bXdV13fl8LoricrlMo1X8w+Hy4a6aprmW1aK4Z7K79moCxyFUwZPVdT38dZiQTqdTXdfLIWlkuaXkyO0oTdOcz+eyLLf5JIQQiisxGtgLoQqebDYztW17Op2Kojifz6OuwOVMsNxess08AfAaDFSHLaqqKt2XM6ar9RZatmLfX1mWnygNgHlaqmC7LpdLTFRxZPq729d1fT6frzVHpeVVVeUaANS27XC3a4pMmW/l9jdUcq3xr8g08undx7q2cfHeUV97oW94noEn6IFnSEOpljeL25RlmZbEdqa6rmd3tbDPuFlZlpfLZbrZwh/Orko7GblcLtf+fDR6bHRcH5UOJz3E0Ozzk+Wrb7qra0/ptadodNTLG8wWP/s8A0+n+w82LUaoj7avzDafxL6/Dw17vyYN+YrKskxdiqfTabapLI4TL4qirusUFEaXMd6g67o4xHtURhyLFn+uqmrY4znc7KPSYw3303XdbBdtWlj+kGoeHvWonrIs09r0pE13cu15Bp7p2akODmplS1U6iY6WXGupGrbfDA2bUj7fUpW+QIbtJcOGmdmNhzXPHtqHDJtwrrX9TLf/fNvY9EBmv1HTxqMmpWuv++zGs0eXlo/KAJ5OSxVs2kcbcmJD1LRlK7ZqZBminhpIRtf/V1WVctW1xqqFX7OUFPU/YkfesUdpmP+o8n6uCy894aMa1k/TujBHa3wdtVTB1ghV8JpGJ/68fX/xh2lkSUum5/vp2KBciWd26FJ8uIwTci5PQz9Nq23b9n3/mdwTrwCdTWzANglVsGk3nJVjnhjOnL4Qg24u6Vqj17VGlPtdsDa75+wPt/xCrHm4eCub4ais28rIePEmkJdQBa9m2hyVse+v+NH8cy1GxOVPP+svtJnd5rZQ1TRNVVUhhBDC6XS64R5BKYqlnTz9uQWuEapg0+IZ9KN5KG6f0lXGvj9Wats2hHA+n1MGKsuyrutr8yzMqqpqGMXSHkzfCtskVMF23dxtF7ePJ+NP9v19dJR03laxm+Xt8Xx3P9NnI82nkK7Ri21O6+tpmiZF6nhVYNqDyT9hm4Qq2K7UtvTRRqbh9tlTznIH3+M7p2bjXfYr4z4UqtKvfd/f3ECYLjacTrPuuj/YJqEKNiqNR742p/ay1AP4ob6/6dl6+ocLw5XW37wlo4WZTrNHyeJKuhpFyYXQ89FnZs3DARshVMGTtRNxYHI6cd4WUIY9gMXqXrDpVFLT83fa1WhS7+E064/sn+q6blT2cEby6cY3P1CaP3202+kFfddyZ9u268eqx4ebbq/vD7brCROOAtfv6TY0O/338ozqQ9f2MzujemrUiaOh67qOS2bnPR+Oto7bD9uEZqcFn71X3We+hWbvJDi638vs9vHobnvQaw83fZbSZvUPo2azNFJqun1cPjzAtJPRIZtUHTZFqILnWA5VozPu0PpQlc7io13Nhqr++p19Zzf+6A2V7xeqZitZuFHPZx60X/0sLVQ1TFdp+9k4OPsmuVwuw53LVbAdoTddL/BD6oKsfsi7fUZpuFj8EltTyXCbz4z6Wn/U17ZsmmZaRlo43PLaHoZXAuoQhI0QqoBdGoUqgKczUB0AIIOfnl0AwH87s9Zv//kOr4/uwdRQwLu0VAHPJ7IAL8CYKmCX0rRYvsSAjRCqAAAy0P0HAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJDBT88u4KFCCM8uAQAOpO/7Z5fwOMcKVcXBXl2A7Qsh+GZ+VUdry9D9BwCQgVAFAJCBUAUAkIFQBQCQgVAFAJDBFkNVCKFt2+UNqqoaLayqKoQwuwoA4N42F6qapimKYiFUzV6fGULoui7+3HXd0a7hBACebkOhqmmaqqrO5/PyNtcW1nXd933f93VdX9sSAOBOtjLlWtu2p9Mp/VrX9TQVDbcpyzK1ZsV2qeGBTJek5Rs5XgAi38wv7Ggv7lZaqqqqGrYzzYqJ6nK5TFeVZXnH4gAA3rOVUPWuOPx8mqhmR1/FjLU82h0AIKN9hKqmabquq+v681f2hVvlOA4A4GXt44bK5/O5LMssY88P1bkLADzMDkJVap0aNlN1XVdVVbxg8N0/BAC4tx2EqihNQzVc0rZtTE6jtdONAQDuagdjqtq27QfiWPWyLPu+jx2Cs5f+uR4QAHikHYSqd8VoFZus2raNg8pN/gkAPNJuuv8WVFVV1/X5fE7X6F0uFwOqAIBHeqmpTuPEVAtx6mhTuwJsn2/mF3a0F/dgR3uwV3f7Qnjr+6/PrgJ4Jt/ML+xoL+4rdP8BwA1Go2+rqhr2dcRZe+L0PdON05IdDeFNR/TsQl7WsSLk0SLz9mmpAp74zTy9W0ZZlukWZyGEOEI3blbX9TA/NU1zPp+Lp84p/fPPP//973//29/+tnL7EMLwAKeyx8SjnXZf4eo/ALhNnKAnTdkTZ5aOq0bXPMUIlWzh9rL/+c9/8u7wfD5v4bj2S6gCgKL4cS15mj66aZqFhDGcZTrNRJ3+cNjYE9u6qqpKe2t+iJ1xyznm+/fvVVXF/8cK4/Kff/65KIp//OMfv/3228Kfx9qmj9I0TbyzbSo1/tB1Xfwh/mGsfEddnM8lVAHAf6UpD4s/x6Y4ofQofwxnmR5lrNSsFULouq4sy67rTqdT3HPc4Hw+V1U1XH5N13VfvnyJP//666/rB0U1TXM6neLP6Yd4mPGmumVZxjKmf3s6nWLlhRas9fojOdrxbl9R/OvZJQBP9sRv5uLP3X9pYV3X8YfL5ZI2i/EibhN/jY1G/Y9bfaQ9pC2HfzJ8uOny+Iizfv/996Iovn37Fn/99u1bURS///57/MN///vfKw9weD+S4SMOi0kbxI3j4V97otY42mnX1X8AbEsIb6u3/aUo/rVy009eFhNbfWKnWNd1l8vl3cab1NITxfaq2x497ed8Pv/6669d1/31r39d/pNYXmpdGzZH9X1fDPr7pn9bVVXf923bLveBMnK4UDW91mOoP9JFCgDb9JH0k/ny4YVJpFPPYNrmo2kjprH483NvUJsqiU1us7kqni7LshyWzbLDjalabrh7dnUAPNPwZrJTcQRS27azkWiasUZ5JUuTz/fv34t1mWw4PmyobdvY0hbbohYSZGqs+lzJB3K4UAUASbzYLYqDyuNwolmj6+NG0iD04fWD6f9xebp276O+fPkS49Svv/5aFEXq+/vjjz8W/iqlwGKSFIdD5pcfWqj6gHsP2tqUox3v9hmoDjzxm3l0QizLMg3N7icD1Yd/En9OA9X7H8O9036uLU8Lh+O+ixUD1YfiKPVUzP/93/996Bhnqx0tn02WC0UuPPpH/2TXjjXV6dGmdt0+M6oDL/PNvDAe61ov2xrfv3//8uVL3/expWo0Pv233357d0b1hZaqYRfhuz/fcAgv8+KudLCjPdiru31CFeCbObrWy/bzzz/HUPXRP1xe9RhHe3GNqQKAHXOD5O04VoQ8WmTePi1VgG/mF3a0F1dLFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABj89u4BHCyEsrD3UdBoAQEaHC1ViEwBwD7r/AAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADI43L3/eJYQ3uYW/zK9w3Xff71/OQCQmVDFg8xGpRBEKABehO4/AIAMDtdSFaa9TQN93z+sEgDglRwuVIlNAMA96P4DAMhAqAIAyECoAgDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyGCLoSqE0LbtaGHbtlVVhRBCCFVVTTcYrn1ImQAA/7O5UNU0TVEUo8zUtu3pdOq6rizLsiy7rjudTnHLKITQdV38ueu65XvRAABkt6FQ1TRNVVXn83m66nQ6FUVxuVzatm3bNt5qJm0Z01Vd133f931f13VaCADwGGEj98KLbVHp17quRw1RxZ9v21dVVdd1l8sl9vqN1k6XpOUbOV6iEN76/uuzqwCeyTfzCzvai7uVlqqqqobtTCNlWc4uH25wt9IAAN63lVC1rG3bYcNV27ZxBNXsiPXiR8aaXcVjhPD27BIA4KH2EaqGmqZJQ6xu+PNwq9zHAQC8lJ+eXcAHDMddxdFUN+zkUJ27AMDD7KalKjVQxav8UqJaiFYmrAIAHmYfLVVN05zP57Isrw2TSpNUzf4KAHBv+whVcUqqa4kqTgc6XXjvqgDYtdGMhlVVDbs44uyJVVXFzabTH15bznH1GxOnTkgzefZ9nwaklxOXyyVtUJblcOO4amSDx/uqiuJfGTcDXtgTv5mn58R4Kklr46kkrhqemPrBBEAPrfjP/vKXv/z73/9+YgHvOtppdwctVamB6lqnXlVVdV2fz+d0jd7Nw9gBOJThwJJ4OVSarGd0Kjmfz6PJfR5X5RX/+c9/nl0Cf/bsVJfT5XKZbaBKXux4t0xLFbDSE7+Ziz83TfU/2p/iz6k/ZPaMOVx4uVyG+6nretisFYejpL2lDeq6Hva6XPP777+XZRn/X5blt2/f4vK//OUvRVEsN1aNHigumT5oPOp3K7nB0U67Bzvag726TyRUASttKlTFMSQpS6UfYjAa3WQ2Lkx/lXaSlvc/slca5ht3mH5NfYgLaeb3339PO0n5rF8XqtIDDccZjzou0z7freQGRzvt7mZKBQC4t9jfd+1eHfGqqbhBWZbvjjOJG/R937Zt3/fFn0e1933fNE1c/m5n4rdv39q2bdv227dvXdd9//79jz/+KIri73//+9/+9rflv40FxCx1uVyapokpMO4w3kh3tkI+SqgCYFtCWP9fWL/xJ8W0EdNP13VrwkfXdcMmotlr1VdK+4nB7ub9jIJgPJy0MMav2/ZMIVQBsDV9v/6/fv3Ga4xCxlBqxIpx6obLoYZ/8rB5f5YfSITKS6gCgP9aDkyxBzD2/U3XTgPKqGkqS4L5/v17kS+TjY50NFMXHyVUAXBcsSMvCiHEAUbXNo6R61rfX8xMcZTScPvUb9h1XRok/lFfvnyJcerXX38tiuKvf/1rXB5HVt1sGCJjhULVZ+xgnioAuJ80/LwsyziL+rUtr912tqqq2IgVd5UaqKbLPzMM/MuXL/GHdD1gURT//Oc/i6L4xz/+cfNuL5fL6XSKEz1+skJCv7Kf+SWEcKzjfaIQ3vr+a67NgBf2Mt/MC+Ox2ra9uQXo+/fvX7586fs+tlSlNqrot99+e/fqvzUWiv+Ml3lxVzrY0b53+cehno27EqqAlY523r3mWhPRzz//HEPVR/9wedVjHO3FPdyYquVpu55dHQB8jFFQ23GsCHm0yPxEmqCAlXwzv7CjvbiHa6kCALgHoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACCDn55dwKOFEBbWHupm2hmF8DZZ9sv0me77rw8pBwCe4HChSmy6h2laCkGEAuBYdP8BAGQgVAEAZCBUAQBkIFQBAGRwuIHqAGzN8nXZsBdCFQDP5KLsF3a0uKz7DwAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACCDw03+uTwRmTnoAIDbHC5UiU0AwD3o/gMAyECoAgDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyECoAgDIYIuhKoTQtu10eVVVIYQQQlVVH10LAHBXm5tRvWmaoijath1lo+HtZbquCyEM50ZfXgsAcG8baqlqmqaqqvP5PF0VA1Zd133f931f13XxI36lH66tBQB4gK206LRtezqd0q91XQ9TUWyImjZNxSXLa4e0YD1MCG99//XZVQDwTEc77W6lpaqqqmE701RZlgt/vrwWAODethKqFswOWk8pamHt7CoAgHvY3ED1D7khNg2HtH/IoRowAYCP2neouoFsBADcw75D1fKUVCaseiKj1AE4mh2MqYrZqOu64cIP/QoAcG87CFXFlYv70sLltQAAD7CPUBXnrIpNVm3bxsHmo8k/r63lBYTw9uwSAOAd+xhTVVVVXdfn8zldu3e5XNKQqeW1AAAPsLOpTuMcCtcC0/La4nhTu74M87MD7NHRTrsHO9qDvbovQ6gC2KOjnXb3MaYKAGDjhCoAgAyEKgCADIQqAIAMhCoAgAyEKgCADIQqAIAM9jGjekZp1vVZh5pOAwDI6HChSmwCAO5B9x8AQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFR8TwtuzSwCALRKqAAAyEKoAADIQqgAAMhCqAAAyEKoAADL46dkFPFoIYWFt3/cPqwQAeCWHC1ViEwBwD7r/AAAyEKoAADI4XPcfG3dlxvZfpmPh+v7r/csBgLWEKrZlNiqFIEIBsHW6/wAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMjjcbWrC9B5yA33fP6wSAOCVHC5UiU0AwD3o/gMAyECoAgDIQKgCAMhAqAIAyECoAgDI4HBX//EhIbxNlv0yOytF3399b1drH9QFmgDskVDFkmlUCuH9/HR9b+9vsz57AcCm6P4DAMhgT6GqqqoQQgihqqqPrgUAuKt9dP+1bXs6nYqiKMuyKIqu60IIw7nRhzefma4FALi3fbRUxUR1uVzatm3b9nK5FEWRWqTiD3Vd933f931d10VRNE3zpGIBgCPaR4tObIiaNk3FJctrR/vZxfFuWQhvtw1UD2HtQPXpZjc/KABPdLTT7j5aqkbati1+dAVGw58BAB5vH6Eq9veFEJqmqaoq9gbGaBX/PyJjAQAPto9QVVVVzEnn87nrumJdbJrNW+FWuY8JAHgp+whVIYSu64ZD0buuu23qhP5WuY8JAHgpOwhVaQRVuqCvaZqyLLuum22LSkxYBQA8zG5C1bWEFJfHPsFk9Ct759I/ALZvB6EqxqZRo1SMTXHV7PgqY9UBgEfaTajquq5pmjj55yhLxW7BlL3ioHKTfwIAj7SbWbmqqhp26pVlOWy7aprmfD6nXy+Xy2x34dFmIbuHp0z+CcAeHe20u7OjXR5ftby2ON6rew9CFQArHe20e7CjPdirew9CFQArHe20u4MxVQAA2ydUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZPDTswt4tHgTm2sONZ0GAJDR4UKV2AQA3IPuPwCADIQqAIAMhCoAgAyOdafDo93ZcVPcUBngaI522tVSBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJDB4WZUZ/sW7yT0P0e6oASAHRCqyGBlDFpjZVTK+IgAkIVQRR7ajQA4OGOqAAAyEKoAADI4XPdfWByMc6jZ9AGAjA4XqsQmAOAedP8BAGQgVAEAZCBU8T8hvD27BADYK6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIIOfnl3Ao4UQFtb2ff+wSgCAV3K4UCU2AQD3oPsPACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIIM9haqmaaqqCiFUVdW27WhtXBXXPqE4AODYwl5mGK+qquu6oijKsow/XC6XlJ+mN5+ZPa4QdnO8TxHCW99/vekPiwc/r49/RAA+6min3X20VDVN03VdXdd937dte7lciqI4nU5xbYxWcW3f93Vdxz95WrkAwPHsI0LGhqhhqTEzxf9P106XpOW7ON4HCOFtbvEvRfGv0aI1bVdaqgCYOtppdzc3VC7LcvjrqCFqtJZ3zUalEMbLQygmPasAwIwddP+lMelpKHoIIYWq6Yj1QsbKqu9X/QcAB7eDUBV1Xdd1XVmWccjU+Xx+d9TUbN4Kt7rDMQEAr2M3oaooijhKvWma2EF7Pp9v28ltch8NAPBSdhOqRj168dfZtqjEhFUAwMPsIFQtZ6O4Ns5clYx+BQC4tx2Eqmg2NsVENTss3Vh1AOCR9hGq4uD0eHeatm1HWSqOWI8L27aNg8pN/gkAPNJuZuVqmmY4Mr0sy+GAqtHa4R1sho42C9lHTW9Ts9k5NjdbGADJ0U67OzvaGKSujbJaXlsc79X9KKEKgIyOdto92NEe7NX9KKEKgIyOdtrdx5gqAICNE6oAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAy+OnZBTxavInNNYeaTgMAyOhwoUpsAgDuQfcfAEAGQhUAQAZCFQBABkIVAEAGx7p99NFul/15IRTbfMI2WxgAydFOu1qqAAAyEKoAADIQqgAAMhCqeBEhvD27BAAOTagCAMhAqAIAyECoAgDIQKgCAMhAqAIAyOCnZxfwaCGEhbWHmvgVAMjocKFKbAIA7kH3HwBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABj89u4BHCyEsrO37/mGV8EkhvP15wS+zr23ff31IOQAc3eFCldj0MkZpKQT5CYBn0v0HAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkMEuQ1UIoWma0cKqqkIIIYSqqp5QEwBwbPubUT3eZ6Zt2+nCqOu6EIKZ0wGAR9pZS9W0gaooitg0Vdd13/d939d1fW1LAIA72VOLTtu2p9Mp/lyWZWqsis1UwwOZLknLd3S8WxBCsc0nbFpYCG/u/QewKUc77e6ppSomqsvlMl1VluXDywEA+J/dhKrYxzdNVKPBVZGMBQA82D5CVdM0XdfVdf2hK/tm81a4Va5j2YgQVv0HAKy0g6v/2rY9n89lWWYZe36ozt1lngkAyGgHoSplqWEzVdd1VVUtxywTVgEAD7ODUBV1XTdd0rZtzFWjtdONAQDuagdjqtq27QfiWPWyLPu+j4lqdli6seoAwCPtIFS9K0ar2NnXtm0cVG7yTwDgkXbT/begqqq6rs/nc7pG73K5GFAFADzSS011GudQWIhTR5vadcFmp0pfae/1AxzB0U67Bzvag726yfQWLnsPJXuvH+AIjnbafYUxVQAATydUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZPAKM6p/SJp1fdahptMAADI6XKgSmwCAe9D9BwCQgVAFAJCBUAUAkIFQBQCQweEGqvMyFq/j/B9XJgDwGEIVu7QyKq0MXgDwebr/AAAyEKoAADLQ/feCQnibLPtl0hH29THFAMBBCFUvqO/HgSmE8UKDjQAgL91/AAAZCFUAABkIVQAAGQhVAAAZHG6gelgcod2bfhsAuMnhQpXYBADcg+4/AIAMhCoAgAyEKgCADIQqAIAMhCoAgAyEKgCADIQqAIAMhCoAgAyEKgCADIQqAIAMhCoAgAzCoe6FF8KxjndBCMURnomDHCbANh3ttKulCgAgA6EKACCDn55dwKOFEBbWHqqVEgDI6HChSmwCAO5B9x/HEsLbs0sA4DUJVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZHG5KhSNYnIoLALgLoeo1mY0LAB5sN91/bdtWVRVCCCFUVdW27WiD4don1AcAHNs+bh/dtu3pdCqKoizLoii6riuKoq7rpmniBtObz8we10Fulx2Clqr/mj4VIbz1/dcnlQNwLAc57Sb7aKmKiepyubRt27ZtfIXO53NcG5um6rru+77v+7qui6JIeQsA4AH2ESFjQ9Sw1Kqquq67XC6x12+0drokLd/F8X6SlqpESxXAEx3ktJvsY6B6WZbLI6VityAAwLPsI1SNhqW3bRuHVc2OWC+KoizLuAFMhtt9nU45caR/RwFwL/sIVUNN08TRVJfLZXnLeMHgaOF0SPtKh2rAfCFvkyW/FMW//rxEbyAAGewpVKVrAIuiiKOpbtiJbHQo0+FTIYwXmisVgCz2cfVfURRN08REFa/yW5OoTFgFADzMPoblxy6/sixnR1C5+m/E1X8Lplf/eboA7uQgp91kHy1VcRDVbKIqrlz653pAAOCRdhCqUpaqJuKqOM9n7Oxr2zY2U5n8EwB4pB0MVE+h6tosCVVV1XV9Pp/TlX03D2MHALjNS3V2xvi1EKcO0rlrkNACY6oAHuYgp91kBy1V62mdAgCeZQdjqgAAtk+oAgDIQKgCAMhAqAIAyECoAgDIQKiCGSG8PbsEAHZGqOJYRpNUAUAuLzVP1Rpp1vVZh5qjDADI6HChSmwCAO5B9x8AQAZCFQBABkIVAEAGQhXcl9kZAA5CqAIAyECoAgDIQKgCAMjgcPNUwdTcsKdfptPE3nU29hDezPYOsGtCFcykpRDc0AaAjxGqoJi7d9HXuZaqRxQDwE4ZU8XR9f3Mf0XxNlkCAEuEKgCADIQqAIAMhCoAgAwON1A9zI1JTnpjZ/iclbMzFHeeoAGAxztcqBKbuCuzMwAclu4/AIAMDtdSBVuwhTncAchLqILb3XxvGb2EAK9H9x/MkG8A+CihCgAgA6EKACADoQpezdwoeADuTqgCAMhAqAIAyECogvtyISHAQZinCtYyYycAC4QqWOuuM3aKYgB7p/sPACCDw7VUhWlvzUDf9w+rhN2Ze+98nev+e0QxycpOyUJjGMCdHa+IZlEAAAgySURBVC5UiU3cZvaNM73332Joz2D6iG4jCLARhwtVsAUrs5d/AgDsiFC1JyGE/ba07br44g71v7uzEHJ2OHr+n2jXxRfqf6pdF39ABqrDRvX9zH9F8TZZAsAmCFWQWWxhWv7v9bjh4PZ5jeDedP/B7abjwTUdARyWUAUbdaVdYTpdggv9ADZBqIKNmp0WYTpdwkt2JgLskTFVS5ZnCn389h+16/p3XfwN+8+4fZZJqj7+/Pxy5/3fcefePHntuv5dF3/D/rdW/969VEtVVVVd1xVFUZZl27bPLgc2bTqP6Po/XLl80qh24yNee7g1e8v7oA+W+zX6X99x2mC/Tw5s0OuEqmF87rrO3B4cx7p/OvbTwVhX7mbz3o7m+yV/cXrejuW+Yy8W3MmLdP9VVVUURV3Xfd/3fV/XdVEUTdM8typ4gNnprOYmuArvTnnlnyEAn/EizTmxmWp4LNMlxcenpt3p9iH89+y4kXoesHPbr9m+KP7152W/TJYUK68lnD7ytJ7ZOz0PHjH9vL7JZLTD2foXHmjB/DbXmnPu9slNB5jqWVP8/GZzN4V8+9FS9bQ357V+yTX136OeB+/c9i/vFY62bdvT6TQaRyVUbaeeB+zc9jdsvzxeZ/G8uP5hV1l3m50P7GpY/JphSZt9sVaPG1sb2p4eqra//aaKOeD2e/ci3X9TZVk+uwR4TSu7F+e2eZvdfp33uy+LovgxYX0/mLz+64oJ7vt3t3nK9muKv7b/uRfOOCq4r9cZqB6HVY20bfvn5dOxusv2un0atr+py2s3VYzti6IIYWkehHvUM3zEe9d/bMvfFVv5ptrk9psq5hW2/9DWe/c6oWrNHApHaoM80KECwBa8bPdfNNt8BQCQ3YuMIFt59R8AwJ28SEvVaFh67Ao0Vh0AeJgXCVVxns/Y2RdnWChM/gkAPNCLhKqqquq6jneniYnqcrkMB1RVVRVCCCFsf5RVvGIxVTsdgL+XYwkhTHPt9otvmiYWudMnf7nCzdYfQpi91mQXhzNb/I4+yNee/OEG0yK3X/8uPss7fed/8u399PrvqH8tl8vlcrmMFu7oqC+XS6ywLMvUfZluv9Pv51jSUUwXbrn49JynH4Zvp43XP/vmGW6w2frjfaWG7/NoueCNHM5s8Tv6IF978pONf5av1b+Lz/L6N89wg6cX/8m399Prv6uXOphZ8SUf3RZw4Rvkua59+OPPezmWWNjoi3j7xY9KSl8c8dft1z9688T600uwzfrrup79Uu7fK3j518dYKH75g7yF4vvF+ofbTD/L268/lpRqHn2Wt1D/+jfP6IO8heI/c57aQv139fqh6lpMflY9y6a1xTdofPvu4ljS99foi3j7xU/rqes6fdr3WP9wyQbrH75Vrp1ari15+uF8tPitfZCX659us7XP8kef/019lnf9zp99xPVv7y3Uf1cvMqbqmt1dBliWZfqn4bUNHlbMbdKYtumq7Rc/qrBpmuGwsO3XPzR982+t/qqq4tfQtff8csHPPZzl4rf/QX73yS+2/Vn+6JtnU5/lNU9+ssEP8iff3k+v/65ePFTN2vIr2rbt8JPftm3XdUVRzI4ELLZ3LHHU4fRbePvFpwrTCMrhQPvt11/8eNpj2VVVxTNirHwX9Q8tF7ywds2dFR7gtg/yRoqPbvgsb6T+mz/LG6l/9oMc699I8TefpzZS/10dIlRduy3go+v4oKZpFv6lOLSRY2mapuu6uq4/dEHHRoqPuq7rui79O+x8Pr87Mcd26q+qKn49nc/n+B23JjZtp/6Vdlfw+g/ydtz2Wd6UGz7LGzH7Qd7sC7G789S9HSJU7e61bNs2hHA+n4vJ3BCb1bbt+Xwuy3Iv31zX9H0f/x3W931RFPFV2IUQQjwRpp6Frut28eZ5VXv8IEc+y0+0lw/yft/ed3WIUDVrs++AFPzjh2pNnVs4lvT9W1VVVVXx1/hdsBxqt1B8MmraWdMuvZH608CL9EI0TVOWZdd1u6h/vfgGW1j7uFIWXfsgb7/4YanTz/L264+ufZY3Xv/yB3k7xd92ntpO/Xd0vzHwGzE9zC0f+OhK4JEtH8tCT1P81G25+Gj6zG/tiq0Fs1cm76X+2eKXC97O4cwWv6MP8rT+fX2Wr715dvFZnha/iw/yZ97eW6j/rl6/pWr0BbHx6wFjU+q1poXZsjdyLG3bDt9Yw7lV4j+5tlx8EkcwjH6N/4TaeP3Vj3s0DRfuqP6p5YI3fjj7/SAXPstPtYsP8mfe3luo/65eP1Tt6LaAw4tWRuKq0bGEEIqtHsvU9ouP//yKz3ZqaU+f9o3XHwvruq5pmj3WP7Vc8JYP57U/yMUe6k8X0O3us1BVVezs2+wH+ZNv76fXf3ePaxR7ntGMGtP72GzEwswfqea9HMtoFuBo+8WPKtxd/dNxJMO1m63/2qzKywVv5HCu9eDs5YP87pTWG/8sr3zzbLP+a8Vv+YP8+bf3Rp78Own93I14XlIM0a8xGm7Xx7L94pcrVP+DvdjhDO26+GIP9e/6zbPr4ov913+bA4UqAID7ef0xVQAADyBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZPD/AQtKvJLQuMREAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h_dimuon_pt0_data.SetStats(0)\n",
    "h_dimuon_pt0_mc.SetLineColor(r.kBlue)\n",
    "h_dimuon_pt0_mc.SetFillColor(0)\n",
    "h_dimuon_pt0_data.Draw()\n",
    "h_dimuon_pt0_mc.Draw('same hist')\n",
    "c.Draw()\n",
    "c.BuildLegend(0.7,0.7,0.9,0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<cppyy.gbl.TLegend object at 0x7f808b709e40>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dUbKjSnouULjhWbmjQeHX7mnUGcLxiwNkj8J7Guc8ugNw2OPiPuQtLgUIsaWURMJasaNChZCUQiB9ykx+5X3fZwAAPOf/fLoBAABHIFQBAEQgVAEARCBUAQBEIFQBAEQgVAEARCBUAQBEIFRBBGVZ5jeUZVnX9eKt6roOK0yWBOuP2LbtsGbbttGeCW8UdptbuweQHKEKXqvruuv1+kD0WV9fkKKua4EMduWfPt0AOJSmaSZL6rruui7Lssvl0jTNuF8qmC8ZbriSnK7X6xPN5AjCPiBXwX4IVRDTPCGFYBSG8y6Xy/iHodZ7GkIUW6SbCmCHDP/BOww9WN/qV7gVnsLyoiiebBUAEQlV8A5DD9Z42K5t21vT2Kuqym4nsJVxn7IsF8cTVx5rPDv+1rzp8d2OZ+XfGrvcaGjSeN79eP7ZuG1lWS6mzPDU7j6FxTU33uHzXYOT+1zJ1vNHnz/xycsxeRW2bxAgsh542tBptLJOyEnjdcKSoigW17l1h0On17BO0zTDtfMl81vNF959ZxjuNu7bSLj58KzHmqZZ7IqbPLWV7rrJYy0+Svbr9t++2rfcamRYXlXVsObKKzJ/lRev3b5BgOj0VMGbDN0JG7s9wqfjvI8h4tjf5XIJF4aP5CFSLNZ0COvPV36yI+R6vQ7BYkgVl8ul67r58vFjtW0bZp6Nc8mtVoXuvXmE7bpu/IpsXG274UyFoijCphvy4nza3PwVGW+QYbV+lpLD3rW4QYY1n+xWBO57PpcBW3qq+lHHTPjvek/V/NrxnYQP3fHl+bVj856qec/ZZM3FrpHFJ/5wL064z1t9ReNksNjgW49+q1Nw3v7FV+Tuat+y+Bz70T4zD4637uTuwu0bBHgFPVWwU6GjZdKZMXSWPN/rELpk5qNdZVne6iRbXHneyO+aPNDw1G4tH9+wqqpn+sn6vm+a5u49bFxtbjw5bHLVfElZllVVrYwA3u0qu7VB9FHBewhVsF/zcBM+Vm/N+3nA4sftraj0os/mh+82zHOfT9Oeh4/xtO75pO/vrrbdSghevMP508m+M7p6a4OYqA7vIVTBmzzQyTRMlBmWhO6liOFmJVRtX/5xIUmEM93CZKzJCkP3W9d1l8vl1s8HbVxtu/UJcLeWh9P3hmf0QKHXuxsEeAWhCt5tezSZjABGHPs7TPnQIXYMW6koisWw0rbt+IzC4eeD5v06W1Z7kZCEQgwaprd/q2NycYNE7NoEVghV8CaPdXiMRwAjjv3tts/pW4ZTFMNUpDBRdGW0KwzqhQlSw2bsum6yNTautsX6nLPJ8vl5gutPZ+7WBjnGyw37J1TBO4zPdf/WDccjgM+M/a0XZ9+48q4MjQxTyL+1WcKgXj8ql/DMarEM1Rwei0HPbBAgCqEK3mEoMvTd/qphBPDJsb9b83IW2zN8uj/wQO+3ZWb3rcLikyUbV9tuuOH2+7xVmWy7+QZJIiXDAQhV8CrhNLTwQythycrZ8itCuAmxbGNH1+SDeeWXcLqu21jRYLcmiWHoFBwL23CeLCfPceNq3zLc57idbdveirnzADQu+7nFfIM8MNUdeMRbq2LBQW3p1JkX5Fwv/jlfOL+T+cJh5TA9uaqqoSDk/J6HhcPKwwNNCl3eert4sqrk4m1v1cCcLx+3Nsx/mrwQ87qa609242qPPc3hPsebPVsqPRrmVM2fzniuVT8qKzos3L5BgFcQqiCC9VB165NsY6jqbxc0Dwvv/hxeeIjFe1hs+bzBi7ddafBGi7fdHqoWe/5C4+fP5VY34eTJblztW261M2y98T3feu3GyydbLxh+AGf7BgGiy/ul3+YEktb+VP60vv54eCuVUb/BUO1z/EwXF2abt8x3N+D2dob7XNnI4+Kl44cON5k0pq7r+R1+a4MAEQlVAAARmKgOABDBP326AcBBPFAp6lVNeZnv/o6eWgZwKnqqgA8wrQc4HnOqAAAi0FMFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQwT99ugFvlef5p5sAACfS9/2nm/A+5wpV2cleXYD9y/PcO/NRna0vw/AfAEAEQhUAQARCFQBABEIVAEAEp5uovj5pzmRJAOAxpwtVYhMA8AqG/wAAIhCqAAAiEKoAACIQqgAAIhCqAAAiEKoAACIQqgAAIjhdnSrFPwGAVzhdqBKbAIBXMPwHABCBUAUAEIFQBQAQgVAFABCBUAUcWZ5/fboJwFkIVQAAEQhVAJxUOzJZPlw7X/9tzYsu6cYnIT9V3aY8P9fzBfL8q+9/fLoVrPngO/OkHHTTNGVZ1nUdwlO4dty2+ZK03NrU//jHP7qu+/d///e3PeJRna6nKl/16dYB8FZN0/R93/d9VVWXyyXLshCqhhWGy3Vdf6B9GzzW/zS+Vdd1Xdc9eYdkJwxV/apPtw6AzwiZqW3buq7Lspwsz36NGnmeD0OEwxfy4ft5uElZlne/tIebz4U7H99DWH9YEloYbn65XIabLD6vyZ2Mb1XX9T/+8Y//+I//+O///u+iKCZXfXMTcr5QBQC3TKJJVVVDF07XdVVV3bphWZZFUfR93zTN9Xodlg/f2Fe6f0JvWXi4YeW2bbuuG24e+s+GJV3XhfwXHrSqqlsZ6Hq9hvsviiIsuVwu4YGqqrper//yL//yb//2b3/961+7rhvucPIs2CjJ3/7L83y8A7VtG/psx0LifnPDAHjedwph/JZl/7lx1Qdm15Vleb1et8xP77quaZpwk+xn/9aQY+4+yvhC+Pyq67rv+/FYZLg2fAKGpBWaF+aBLYaqoccrNCl8VvZ9H+5/HpvCM128ii3SC1WhA1NgAjiq76SfyGchlGU576wKeWWlm+oVQn9BePSh42CIWXmeN00TRirDCtu7EvI8L4qiLMvQUzW+qizL0CE3v4otEhv+W0niw2TDQOoCYIvhXL8QNSbXhpwRhsYmt8pG3/CLohhmZWWjzqdnWpVl2XiCV2jk0HcV/hsa3zTNeKb5uPHjuxqWh56t+Qdl6G8zm+px6xO3dyX0rAZh0DcI/atb7iGt5ws8L8v+89NN4I4PvjOPPw2Hj5WqqsLlbDTbKTQydFYNF8afpONPqDA8VxTFcJ/Z7Jt/EG41WSfccHyHw6PcfdC5eVMng5JVVf3Xf/1XlmV//etfJyuPP2of3sJP3kNaUiogEQb+mqaZ9HOGrxFN0wzfD259RThbwQxAnar9S/eduW3bycfNfEmshxhfyH7tCdvyoJN11m8S8Vmk++I+JplnOySnLMsmoWrxVNXF53W2VxcQqvbvJO/Mi7O/Q7nRWA8x/zT8+DlbJ3lxB2nMqarrOkydW9n5hp7P0LG50ln1mNc8MwBOIZzNNxG3W2t+/6YXv1kCZ/+1bXu9Xoc5gHOTpB/m7i1O2ctS/nkBAGDPEuipGrJUmCwV/tt13XDi6zzph84qCR0AeJsEeqqCec9TqCcbolXcHlQAgO9KbwZZqIc2TL6b/DfIb/yQ+NlmzAEmqu+fd+YDO9uLm8Dw37rQRxXKsoUaaGHJm0vfAgAnl8zw34q+7/M8v16vw9mqKz8tCQDwCumFqrIs532J4cTR0E1lfhUA8H7nGuw82+AuYE7V/nlnPrCzvbjp9VQ9ab2M56leewAgotOFKrEJAHiF5M/+AwDYA6EKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAqEKOLs8//p0E4AjOF3xTxXVAYBXOF2oEpsAgFcw/AcAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEMHp6lQp/gkAvMLpQpXYBAC8guE/AIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAhOV/xTRXUA4BVOF6rEJgDgFQz/AQBEkGSoyvO8ruvJwrIs8zzP87wsyw+0CQA4t/SG/8KkqLZt5wuDruvyPDfMBwC8U2I9VfMOqizLQtdUVVV93/d9X1XVrTUBAF4kpR6dtm0vl0u4XBTF0FkVuqnGT2S+ZFie0PMFnpfnX33/4/l1eB3vzAd2thc3pZ6qkKiapplfVRTF25sDAPD/JROqwhjfPFFNJlcFMhYA8GZpTFSv67rruqqqyrJcTFGL2radnwm4Xvxzxak6MCFFef61tPi3+UFvsI9g/IEy/rwIHx/h2snyyRIYSyBUtW17vV6Loogy91w2gqNajEp5LkJx0zBPN2iapizLuq7bth1m8Y4/NeZLYCyB4b8hS5VlGXb3LMu6rrvba+XLBADrmqYZzhwPmSmEqmGF4fJuTyrfPoDDqyUQqoJuZFgyDPANC4er3t9CANIVMlPbtnVdj7+TD1lqnF3yPA//bdt2mFWS/xRuMpSkDhYfNNx8bhhnnNx8UuY63PxyuQw3mRvf7eLNx4/IkxIIVW3b9iNhrnpRFH3fhx13cVq6ueoAfNckW1RVNf4mH+ogLirLMnwwNU1zvV6H5eGTa37PY6G3LPtZcDGs3LZt13XDzUP/2bAkdCvUdR0etKqqlY60sE5RFOHmw5O6XC7hEWNNsCGBUHXX8J0gG31psH8AJCrPt//l21d+wPDJcvczpeu68YfRynf+W48yXAiFGMuyDH0Hk2tDT1jf92FKTJgMsz4fZmh8aM/wpLJR/5yeqiiOEKrKsgy5O/SCZj8nG366XQA8ou+3//XbV95i/tkROoHatl3ppnqFxT6C0M8UrgpRL4zeXC6XJz/1hKoo0gtVIbxPXv4Q25umCZ2oEhUAG4V+mjzPi6KYf3yE3qChI2p8q2yURYYRtFhlF4aepOGuQiOHefRDtGrbtmma704mHveoDSeB8aQESipsJ0sB8C2Lv342tvjJUlXV9Xodz52q6zpMGA/XPv95VJbl9Xodpqhfr9eQe4YlIW9dLpfQjAc60pqmGW6++GslfNe5fpTnbuXPU20NOAO//bd/6f483LzE9GLR6SgPMb6QLZUqff4hXiHdF/cxJ3u2J3t1AaFq/07yzlzX9bhnK4g7A3jecXCr++1tTvLiDk72bE/26gJC1f55Zz6ws7246U1UBwDYIaEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAggkP9TM0W60XVT1VOAwCI6HShSmwCAF7B8B8AQARCFZCkPP/6dBMAfiFUAQBEIFQBAEQgVAEARCBUAQBEIFQBAERwujpVin/CyS2dNvjb4htD3/94fXOA4zhdqBKb4OTmUSnP5ScgAsN/AAARCFUAABEIVQCbqOEOrDvdnCoA9mb9FCJIhVAFwCc5f+jAzhaXDf8BAEQgVAEARCBUAQBEcLo5VSqqAwCvcLpQJTYBAK9g+A8AIIJkQlXbtmVZ5nme53lZlpOr8pnJOgAAL5XG8F/btpfLJcuyoiiyLOu6Ls9zA3kAwH6k0VMVElXf923btm1bVVWWZUNfVNu2WZY1TdOPhIUAAO+RRqjKfvZRBXVdZ1nWdV34b8hPxvsAgA9KY/ivaZrxf0OKGsessHBIVwIWAPBmic1MCn1U1+s1y7KmaUJ4Wiw9tfi8zMSCw8jzr77/EWW1iHcFjJ3tYzeZ4b/ger2GRFUUxaQ7qqqqMJsq9GDd6qyanye40cufGwCQssRCVd/3TdMURdF13RCbwhT10ImV/RwcHGZcze/hMW94dgBAuhILVVmWlWU5iU3zTqnQWeUEQADgbRIIVXVd53k+SUjj2CQ8AQAfl0CoCh1Rw+heELqpQq/V5XKZdFYN176piQDA6SUTqrquq+s61E0IS8YT0ufXhgKhAADvkca5jsPP1AyKohiP+k3OzquqatKzNayWxPMF7lJSAfbvbB+7aRT/LMty/Msz83G94RdsVP4EAD4ijVAVrKclcQoA+KCUQlUU62U8T9VLCQBEdLpQJTYBAK+QwNl/AAD7J1QBAEQgVAFHpggC8DZCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEpyv+qaI6APAKpwtVYhMA8AqG/wAAIhCqAAAiEKoAACIQqgAAIhCqAAAiEKoAACIQqgAAIjhdnSrFPyFFef41W/bb/Gju+x9vaQ7AgtOFKrEJUjRPS3kuQgH7YvgPACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAqEKACACoQoAIILTFf9UUR3YQg134LtOF6rEJmBiMRip4Q58l+E/AIAIkglVbduWZZnneZ7nZVnOV1i/FgDgpdIY/mvb9nK5ZFlWFEWWZV3X5Xk+Hsgbz5SaXwsA8Gpp9FSFRNX3fdu2bdtWVZVl2dAjFS5UVdX3fd/34dq6rj/UWADgjNLo0cnzvCiKtm3HS7Kfs87Hl+fXTu4niecL3JXnX2+eNv7+R4TUne1jN43hv6Zpxv8N6SoMBQbjywAA75dGqBpG+sKg3vV6HS6Pu68GRVF0XfemxgEApBKqBiFOZVlWFMXds/zCCYOThevFP1ecqgMTeNjGUUKDiXA8aUxUH/R93zRN6Ih6rHRC/6jYTwUAOJTEQlWWZWVZhiG/uwN8ClYBAG+TQKiq6zrP88ncqTAzfRjgmwQsE6oAgDdLIFSF2DSpOxViU7hq8dQ/5wMCAO+UTKjquq6u61D8c5KlQt4KC9u2DVPRFf8EAN4pjapcw8/UDCa1QOu6Hk4MzLKsaZrFCVVnq0IGB7bb4p/O/oPB2T52U3q2Q4q6NQM9rLAyP/1sry4cmFAF+3e2j92U6lTdPZvP6X4AwKekFKqiWC/+eapADQBEdLpQJTYBAK+QwNl/AAD7J1QBAEQgVAEARCBUAQBEIFQBAEQgVAEARCBUAQBEIFQBAERwuuKfKqoDAK9wulAlNgEAr2D4D3iTPP/6dBMAXkioAgCIQKgCAIhAqAIAiECoAgCIQKgCAIhAqAIAiOB0daoU/4SjWj24f+FAB17hdKFKbIID23J8b89eAN9i+A8AIAKhCtiXNxRez/NNfx+nBj2k5XTDf8DJbZwCsIdQBaRFTxUAQARCFQBABEIVkKS+//HpJgD8QqgCAIhAqAIAiOB0Z/+pqA4AvMLpQpXYBClS4ADYv9OFKiBRvhABO5fMnKq2bcuyzPM8z/OyLNu2HV+Vz5Rl+bG2AgDnk0aoatv2crl0XVcURVEUXdddLpe6rj/dLuD/OecPqtz4cZsfO/zFG+AN0ghVl8sly7Kmadq2bds2zIu6Xq/h2tBr1TRNPzLuygJ4kb6f/mXZ12wJcApphKpgPKJXFEX2M06Ff433AS+l3CiwLo2J6kVR3M1MoRMry7KyLAUsAODN0ghVk7G8tm27rst+9k6Fy2GIMPs5LKh0AgDwTmmEqrG6rkNsappmvLyqqjB1vSzLrusmZwgO1ot/rpDSgLmlGfq/zd5mfmxbzQgjpC2lUBXOAQyXm6YZxvjGl7OfFRZC99WcbARENI9BeT5dOF9yayGQtGQmqtd1HRJVVVV9349T1HwG1XgaOwDAG6TRUxWG/IqiWMxJoS7ou9sEADCSRk9VmER1K1FdLpdJqBpPYwcAeIMEeqqGLDUPSXVdDycADpfDdPWqqt7WQoAPyvMv07NgD1IKVStzz/M8v16vQ4314UxAAID3SCBU1XV9NyGF36UJk6uM+gEA75dAqNpInAIAPiiNieoAADt3nJ6qjdYrqisNCgwe/f0F4KROF6rEJmALbxXAdxn+AwCIQKgCAIhAqAIAiECoAgCIQKgCAIhAqAIAiECoAgCI4HR1qhT/hPfI86/Zst8Wj7++//H65gC83OlCldgE7zGPSnkuPwFHZvgPACACoQoAIAKhCgAgAqEKACCC001UB3i/pdMef8wXOpEGkqanCuC1+n7hL8u+Zktea6nIBRCTUAUAEMHphv8U/4Rd2Vwj9BQFrmZbQ7lUSMnpQpXYBLuysUbo6reh45g8ceVSIS2G/wAAIhCqAAAiEKoAACIQqgAAIhCqAAAiEKoAACIQqgAAIhCqAAAiOF3xTxXV4Uk3fkJuofa3wpXAqZwuVIlN8KTFqKT2N4DhPwCACJIJVW3blmWZ53me52VZtm07WWF87QfaBwCcWxrDf23bXi6XLMuKosiyrOu6y+VSVVVd12GF8UypruvyPDfMBwC8Uxo9VSFRNU3Ttm3btiEwXa/XcG3omqqqqu/7vu+rqsqybMhbAABvkEaoCsbjeqHLKgwCdl2XjVJUuDBELgCAN0hj+K8oivWZUiFjAQB8ShqhajItvW3b0Du1OGM9y7KiKMIKAADvkUaoGqvrOgztNU2zvmY4YXCycL345woz3zmYPP9SWSpFS8VXVV6FXUgpVA3nAGZZ1jTNY6UTZCPYlaWvOT8e/e6zCxvTzMOhZ35DlVdhJ5KZqF7XdUhU4Sy/LYlKwSpIQt//8pdlX5MlP5cD7FoaPVVhyK8oisWan9nPEwAHJlQBAG+WRk9VmES1OCc9u3Hqn/MB4bPyfPoXxvVmCwEOIoFQNWSpciZcFQpThS6rtm3DVHTFP+HjjOsBp5LA8N8Qqm4N6pVlWVXV9Xodzux7eBo7AMBjEuipquu6v2FITmGdpmmaptk4jR0AIKIEeqq2k6UAgE9JoKcKAGD/DtVTtcV6RXWlQQGAx5yup+rW9Kzg062D3Vn6URQAFpwuVAEAvIJQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAW8Sd//+HQTAF5I8c9fKFUFHMCN6mK/zd//JF2I6HShSmwCDm8xKuW5CAWvZfgPACACoQoAIAKhCgAgAqEKACACoQoAIAKhCtgXZ6gNbApIi1AFABDB6epUKf4JALzC6UKV2ATs1uqXvv/P2xjs0+lCFcA+bYxKG4MX8H7mVAEARCBUwdHc+DFdiMY+BouEKgCACIQqAIAIhCogAmUqAYQqAIAIhCoAgAhOV6dKRXXOael0rd8WjwYDeQCPOV2oEps4p3lUynP5KVWzNPxjno+91cH7Gf4DSEnfT/+y7Gu2BPiA9EJVnudt246XtG2bz5Rl+Zn2AQCnlFioqus6y7JJqAIA+Lhk5lTVdd22bdd186tCxmqaRu8UcELmxsFOpBGq2ra9Xq8r12ZZJlEBAB+URqgqyzKctVfX9a101bbtkK4ELADgzdIIVevCmODlcgn/DalL6QQA4J2OEKqCqqrCNPayLLuuK8tycT77evHPFVIaZ7CxRqhJPABzRwhVkynqocLC4pT2TDaCVWqEAjwssZIKi+YzqIqiyFReAADe6AihSngCAD4u+VDVtu3lcpl0VoWxP+cAAgBvk3yoCsmp67pQHbRt27CkqqrPNgwAOJUjTFTv+z7P8+v1OpSwGs4EBAB4j8RCVV3Xi2mp7/uhm8qoHynK8y8n2QEkLbFQtUKcAgA+KPk5VQAAe3CcnqqN1iuqKw0KADzmdKFKbAIAXsHwHwBABEIVwCk4vRReTagCAIjgdHOqgCetnuwBcF5CFfBtzvcAmDP8B6+V51+7vTeYsIPBM4QqAIAITjf8p/gnB7PUtfDb4m7u5C8eYx+DjU4XqsQmDmb+MZbnPtuIyT4GGxn+AwCIQKgCAIhAqAIAiECoAgCIQKgCAIhAqAIAiECoAgCI4HShKl/16dbxcn6Fg5PI801/QESKf8JJKd54YBvf54QqiOt0PVUAAK9wup4qYMVS18UP/RkAW+ipAn7R97/8ZdnXZIkhdIBFQhUAQARCFQBABEIVAEAEQhUAQARCFQBABKcrqbBeNl1pUN5jqbD7b/N9cw/1OffQBoAknC5UiU3swTyp5Ln4wtHk+VesvXrjXcVdDb7L8B8AQATphao8z9u2nS8vyzL8KHJZlu9uEwBweomFqrqusyybh6o8z7uuC5e7rlufOAUAEF0yoaqu67Isr9fr/KrQNVVVVd/3fd9XVZX9jF8AAO+RRqhq2/Z6vQ59URNh+ZCiwoXF+AUA8CJphKqyLMe9UHNFUby5SQAAY2mEqhWLk9ZlLJjI801/ADzsyHWq2radnwn48Bx2Ba7YYmNVz2xDUc2lu1q+t40Vd+zCzM12swd3VyA7dqhaJBvxUhGrei7eSo1Q4prsTnYweEbyw38rFKwCAN4m+VAVktPkxMBb5wkCALxI8qEquzEt3Vz1s7kxAwkA3uQIoSoUpgpdVm3bhqnoin8CAO90hInqZVlWVXW9Xocz+5qmMaEKAHinxHqq6rru+37eCxWWN03TNE3f9xIVZ+bsLV7NPgaLjtBTNZClAIBPSaynCgBgnw7VU7XFekV1pUEPZmN9c2MZADzvdKFKbDqViPXNAWCd4T8AgAiEKgCACIQqAIAIhCpY4EdvAO8DfJdQBQAQwenO/gNgMKsw8mOx7IzTpmELPVUAJ9X3078s+1paCGxyup4qxT8BXiFird2Nd7Xk8dVUsON5pwtVYhPAK0SstfvwXcVdDb7rdKEKDmbb93IAXk6oguTpfgXYAxPVYRcMRgCkTqgCAIjA8B/s1I3JUstlhAD4OD1VvISfd4hCGSGAhAhVAAARnG74T/HPY9heG9AEcADe43ShSmw6hohlBoG7Ns7k8/7KyZ0uVAHwLRujUp5v+nlmwYsDM6cKXitu/5neOF7qmR1sy3kVcGxCFQBABEIVAEAEQhUAQARCFQBABEIVAEAESirAB/j9PoDjOV2oSqWiep5/bTm3eeNqrLjxM4UL9dmXKo4+vv13s68BT70PwOB0oWo/sYmPmL8hLr5Fqs8Op+J9gCjMqQIAiOAgoapt23ymLMtPtwsAOIuDhCoAgM86yJyqtm2zLGuaRu8UAPARhwpVEhUfp1YCwGkdJFQFbdsO6UrA4lOcYApwTgcJVV3XZVl2uVzCf6/Xa1EUIWABALzBQUJVUFVVXddZlpVl2XVdWZbzXOqgmbEAAAqJSURBVLVe/HOFAlcrlurmKZoHZxfxkN94V3FXg+86SKiaTFEPFRZC99WEbPQKS6XGvW0BCzZ+sfVWTYoOUlJhPoOqKIrs5wR2APag7zf9QaIOEqp2Ep5u/HrU50VsWNznuNstBrzN4vtAnm/6+7iNb2K7fa/bbcMSdYRQ1bbt5XKZdFaFsT/nAAIkR4cWiTpCqArJqeu6uq5DVYWwpKqqzzYMADiPg0xU7/s+z/Pr9Xq9XsOS4UxAAIA3OEioyrKs7/uhm8qo3wk52RCAzzpOqMoUUgcAPucIc6oAAD7uUD1VW6xXVP9UadCNFcmXLK9mLCy6Gy/Hjz2c1A3AHpwuVO2zovrDFckVLn+n+b6T51+2PwDB6UIVp+IHMQB4G6GKg7sbmIzfARCFieoAABEIVR/gt5a+69VbbOkHxX7s7SfGgDczY5LvMvzH2ZlQBUAUQhUAR+aEFd5GqCJJxuOA7ZywwnucLlRFLP65sWLn+0flIzZse1XSu/d2Y17U41vM10pITvT3gTeLW6g56U8HFp0uVEUs/vlwxc5Xi9iwl97VM/cGpCj194HUCzXv9mPrMJz9BwAQwel6qviWpW7t5V+7m/QA+uoDvIG5UOyKniru6Ptf/rLsa7LE3CbgI+bvRbf+4D2EKgCACAz/EceWTnjfFwE4MKGKCLakJVMfADg2oeqANs4uf3+/0fZp74+uBvAghdd53ulCVcTin9se7jOV7ibPI8+/Jvf//oyyuGmXGra8xbLsP2cLnWAIxLHxvX8n3+7m75yPrUN0pwtV0WPTvYd7sNLd5qriRztm0qoNuM9WAa8zy1WPd7fr9DoeZ/8BwCZLxRoWqsxsWY1DEqoAACI43fAfAOyBqfHHI1Qt2D4V8c37etw5kjuZcQlwQmlNjWcjoWpZrMJLN9Z5fMZirBj3iXoK+ZvPEohL+z8r6fYn3fhM+2Ezc6pebuOMxTy//wfA2Sx9Fvx43QfElg8jH0m3nCu/xx3Xy/P7qy2u83D5kKS/bz3T+D0UXEl642fa/1FJNz7T/o9abPwzpQcf+9h6YrWEN/4DztZTlftJcwAOZmOtB17tdHOq3lxRHQCesfQ7Ews/y7FtHfWKX+t0oUpsAiAh88kPW35kYnEdE6Fe7WzDf9+z3q2VbZjQ9+T9P+m79//q9V9657tq/AP3r/0RJd34B+5f+yNKuvEP3P/e2p+6Q4Wqv/zlL3me53n+l7/85Q0PtzRcvTBnCwA4g+MM/43j8//+7/+e7YwDALhr1tPU63uK6CA9VaFr6vfff+/7vu/733//Pcuyf/3Xf/10u5Z9vDpAcmwxgOfFHWDZVsvqXL0bB+nOCd1U4+cyX5J9v2CG9SOuv6vGWP/Y6++qMdZPa/1dNebW+iul+4b1N5b329vzTd0Reqr+/PPPLMv++Z//+dMNAQDO6wihapGMBQC803Emqi+e8ffnn3/+7W9/Gy/Z2+mmp1p/V42x/rHX31VjrJ/W+rtqzK318/y3u+uvrBO9PQTHCVX/8z//c3edU43sAgDvdNjhv2DSTQUA8CIHmZa/8ew/AIAXOUhP1WRauvMBAYA3O0ioCnU+w1z1P//88+9//3u24+KfAMDxHCRU/e1vf/v999/Dr9OERPXHH3/8/e9/L8tyvnJZluEnAhev3YOVFrZtO762bdsPtG/Vxs272+2/3v66rsMK+9z42b3273znn+ze8xX23P71xu//yL278Qcpbv9s3wdvunt+nueLGzPRpxNBfyx//PHHH3/8kWVZURRFUcyf4863QNM0oVWL7V+8tqqqz7X3F+uNHxtWe28D77jb/mHhcKFpmg81dsEDO88HWzs3buHi4bnng3e98Ts/cvsNG3+w/4N3sf17Pngf2Hk+2NqxqqoW9+T1Q3XPB/LzDvVk+p+v1nC0hN1xOP7DHjnsAbd2iA9ab//8vWBXe+R64wdhs+/wfXm9/aHZw3+Hd7oPNPSGb+08t16dD5psz8nhufODd73xOz9y+3vtnyzf257T32v/zg/eb+08Ozlyq6q69fVg/VBd/+8B7GWviuVWKN5y7R58t/1h993JV64tm3d4O9vD+8LEdzd+VVW7ei9Ieudf/KhIpf0PNH5XR+7d9o9X2+HB+8D238/Bm+KeP94TslkkSvqN6HkHmVN1SxjrHfepji/v36T9RVEM3xT3b77xsyy7XC5Zlk2OyX26u/PUdV3X9Xsb9Q2L23/jte8X2jOZYDFp3n5aO3G38Ts/crds/GzHB+8DO89+Dt6NG3+y/mePhbIsQ4C4tVevN2+3B3Icn8100Q0H/LhzcnzV5AtBWOcDDb1hpf0rK7+zhSvuNn74dr6THuyJLTvP+O1gJ990B+vbf/HanfSU9Kvf12/tMPs5eNcbf2v9nTS+39b+PR+8G3eefR68dzf+no/c+eDd+qG6cu1OntHz9nJURzT/ghKWp/Jy3mr/xPAVIZXGjw+/Hb4vB+s7z7Bw2Pj7eWsO1neejbvWpyx+iq98kO/q4F1p/GTNfR656+3f/8F7d+fZ88F7d+fZ7ZH7rVCVxIH8vKOFqsnRMp6fmMTLudL+wfg9Yj8t7+81PrsXcD/u7s6T7XsqwN3tf3fX+qxb40pJvBevNH5xnZ00e7De/v0fvHd3nj0fvHc3/m6PXKFqbi97VRQPv5zvbOSKLTvcDr9mBeuNDxeKkfAsiqLYybGU+ntB6u0fVD/120YNPtHGmxYbP1y1zyN3bLH9+z94B9/deXbV/hQb/93hv/7XgD6/9gD+aTEgJ2pxxt8gLO+6brxw8t/PWm9/lmV1XV+v16Io9la5LtvQ+Gxpa3ddF6revapZm23ZefZsvf1bXp2dGE8fHnaY/R+8wWLjs30fuWO32j//b7ang3ewsvPs32LjEzpygy2H6v4P5Kd8OtXFtJiRx09znojn63/Q3fbv+SW72/i7K3/WAxt/Vy/Hevu/9ep8yqQ9kzbv/OBdb/zeNvXcevvHdnjw9t/f/rt6RVYav/Mjd7HK1PqhuvMD+Xm7eGEiCntbVVVhyGPouA7XznfWbB+dqIOV9g8NLmZ28hTWN/7Ynt+X13ee7Odo2sqz+5T19hc/K/Lttv1Dk5qmGQbLhmt3fvCuNH7/R25/b+OP7fPg3bLz7PbgXW/8no/cxVC1fqju/EB+3tFCVX/vRInq17oaO3wtb7W/ul3nZj/PYn3jD/b5vtx/c+dJrv0bX50PmrRwsmPv/OC91fgkjtz+3sYfpHLwru88e2v/euN3e+Quhqr+3qG68wP5SXm/9EM8B5D6FJP9t3BF0o3Pjr7z7Lz9bduGmToptv9u43fu8O1PeufZc+MXJf1G9LDDhioAgHc6+M/UAAC8h1AFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQwf8F75CWsGFhJwkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h_dimuon_mass_data.SetStats(0)\n",
    "h_dimuon_mass_mc.SetLineColor(r.kBlue)\n",
    "h_dimuon_mass_mc.SetFillColor(0)\n",
    "h_dimuon_mass_data.Draw()\n",
    "h_dimuon_mass_mc.Draw('same hist')\n",
    "c.Draw()\n",
    "c.BuildLegend(0.7,0.7,0.9,0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "for hist in h_mc_list+h_data_list:\n",
    "    hist['h'].Reset()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}