Poisson

Poisson

{
 "metadata": {
  "name": "poisson"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import numpy as np\nimport ROOT\nimport rootnotes",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "def compute_poisson_interval(n):\n    a = ROOT.Double()\n    b = ROOT.Double()\n    ROOT.RooHistError.instance().getPoissonInterval(n, a, b)\n    return a, b\n\ndef sum_quadrature(numbers):\n    return np.sqrt((numbers**2).sum())",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Definitions\n============"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "NEVENTS_PER_CAT = [1, 2, 3] # first category 1 event on average, second category 2 events on average, ...\nNCAT = len(NEVENTS_PER_CAT) # numbers of categories\nNTOYS = 200000\nweights = np.random.random(NCAT) # weights for every categories",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "One experiment weighted\n=============="
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# number of events per category\nn = np.array([np.random.poisson(l) for l in NEVENTS_PER_CAT])\n# sum of events\ns = n.sum()\n# sum of weights per category\nnw = n * weights\n# sum of weights\nsw = nw.sum()\n# average weight\naverage_weight = sw / s\n# average weight ** 2\naverage_weight_squared = (weights ** 2 * n).sum() / s\nprint \"entries: \", n\nprint \"weights: \", weights\nprint \"weighted contings: \", nw\nprint \"sum of entries: \", s\nprint \"sum of weights: \", sw\nprint \"average weight: \", average_weight\nprint \"average weight squared: \", average_weight_squared",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "entries:  [0 2 1]\nweights:  [ 0.51122862  0.06631955  0.46076829]\nweighted contings:  [ 0.          0.1326391   0.46076829]\nsum of entries:  3\nsum of weights:  0.593407387243\naverage weight:  0.197802462414\naverage weight squared:  0.0737013276373\n"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Computation of intervals without weights\n-------------------------"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# per-category poisson interval\nintervals = np.array([compute_poisson_interval(nn) for nn in n])\n# low and up error\nelow = n - intervals[:, 0]\nehi = intervals[:, 1] - n\n# quadrature error (across categories) of the errors (lo and hi)\nerror = (sum_quadrature(elow), sum_quadrature(ehi))\nprint \"entries: \", n\nprint \"unweighted intervals:\\n\", intervals\nprint \"-> error of the sum (summing the errors): \", error # this is what I am doing for the paper\ninterval_sum = compute_poisson_interval(s)\nerror2 = abs(np.array(interval_sum) - s)\nprint \"sum of entries: \", s\nprint \"interval of the sum: \", interval_sum\nprint \"-> error of the sum (from the sum): \", error2",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "entries:  [0 2 1]\nunweighted intervals:\n[[ 0.          1.14787446]\n [ 0.70818544  4.63785962]\n [ 0.17275378  3.29952656]]\n-> error of the sum (summing the errors):  (1.533988646638458, 3.6828985290487513)\nsum of entries:  3\ninterval of the sum:  (1.367295313890434, 5.918185832883396)\n-> error of the sum (from the sum):  [ 1.63270469  2.91818583]\n"
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "I have computed the error (up/down) in two ways: for the first one I compute the error for every category as a poisson interval, and then I sum in quadrature them (without weights). In the second case I compute directly the error, as poisson interval, from the sum of events in the various categories. The results is different.\n\nSo if the error is defined as\n\n    error = poisson interval bondary - central value\n\nthe sum in quadrature of the errors from the various categories is different from the error of the sum of the events\n\nso what I am doing in the paper (weighted squared sum of errors in the categories) is wrong."
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Toys\n====\nI generate `NTOYS`. I fill the histogram with the sum of the weights."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "histo = ROOT.TH1F(\"histo\", \"histo\", 100, 0, 20)\nfor itoy in xrange(NTOYS):\n    n = [np.random.poisson(n) for n in NEVENTS_PER_CAT]\n    b = (n * weights).sum()\n    histo.Fill(b)\n    \nf = ROOT.TF1(\"f\", \"TMath::Poisson(x, [0]) * [1]\", 0, 100)\nf.SetParameter(0, histo.GetMean())\nf.SetParameter(1, NTOYS / 5.)\ng = f.Clone()\ng.SetLineColor(ROOT.kBlue)\nhisto.Fit(f, \"\", \"goff\")\n\nprint \"expected sum of weights\", (NEVENTS_PER_CAT * weights).sum()\nprint \"average sum of weights\", histo.GetMean()\nprint \"poisson fit\", f.GetParameter(0)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "expected sum of weights 2.02617259145\naverage sum of weights 2.02793785993\npoisson fit 1.89355350524\n"
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "canvas = rootnotes.default_canvas()\nhisto.Draw()\ng.Draw(\"same\")\nf.Draw(\"same\")\ncanvas",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElE\nQVR4nO3dfbajxoE3YHiP9xJ3sgihFaRnkjW0lzDOiScSyeQ4swT3GjyT9gqEFjFxezW8f1Tfcl1A\n6OMWEojnOT4+EipQwZXEr6uKomzbtgAA4G3+36MrAADwDIQqAIAMhCoAgAyEKgCADIQqAIAMhCoA\ngAyEKgCADIQqAIAMhCoAgAyEKmAq3333XVmWZVmeLfm73/3uwpIAsyVUAQBkIFQBj/d///d/bdve\ncCvSP/7xj9q3gJkQqoAF+5//+Z9HVwHgi68eXQFgLb777ruff/65KIp3797913/9V+fVqqqapumv\n9ac//enz589ff/11URTff/99+tJf/vKX9MFf//rXzrrffvvt58+fi6L4+uuv//GPf2TZC4BTyhva\n2wEu8d133/39738viqJt29/+9rchUUV//vOfY7T63e9+969//SuUjAUG+/W+/fbbEK36G/zDH/7w\n448/hsdff/31L7/80ln397///adPn964UwCnaKkCJhcC0L//+7+/e/fu559//t///d+iKP7+97/3\n26uimKjev38fmqk+ffr0+fPn0OD0/fff//GPf/zXv/4Vuv/+8Ic/hHcJq8RE9fvf/z6s+/nz559+\n+umnn356//69XAVMpQWYxp///OfBn5q4PC6JeSg8/fbbbwd/oN6/f18Uxddffx2X9Iv9x3/8x+C6\np5YD5GKgOjC5NF0VRREbqL777rurtvPPf/6zbdtOr19HaIj6zW9+01kex1SFZAaQnVAFTG6km29Q\nHJD+7t27P/3pT1etGzr+Qq9fR0haYeg6QHZCFTBHoQcwDKIqy/Ldu3f/9m//dvnqg6EKYFJCFTBH\n33//ffsyiKoois+fP3/69Kksy6uiVYekBUxKqALmKwyiStPVp0+fbs5VOv6ASQlVwAKEdBWi1SVz\nIshPwP0JVcAcDTZH/fOf/zy74sho9DCG3dV/wESEKmB23r179+nTp3fv3nWWn+r4i/erKV4y0y+/\n/BInu0qXF8ncCgB5mVEdmJ2ff/65LMvPnz+XZZmOVQ/tT520VBTF3/72tx9//PG3v/3tjz/++I9/\n/CPMn/7f//3fnz59ijOqh2aqOAUoQHZCFTBH3377bbgvTWcEVbz3X/Cf//mff/vb34qiSGcE/fTp\n0/v373/66adffvkl3gHwN7/5zfv37zVTAdNxQ2Vg7sL8n2mW6vjLX/7y17/+dfCl0KwlSwF3IFQB\nAGRgoDoAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZmVAeAG5Vl+egqMK2rpvMUqgDg\ndubQfmLXhmbdfwAAGQhVAAAZCFUAsCJN0+z3+7cU4BShCgBWZL/f13U9UqBpmrqum6a5V42eh1AF\nAPxqv98fDoeqqi4s7BLISKgCAF65MFEVRaFBKyVUAcAaNU1TvUizUVieltzv97FkOtyqqqrj8Rge\npMvTLa9qeJZQBQCrs9/vt9tteHw8HrfbbcxVTdMcj8f4tCzLdAxWXdeD/X2xfFVV2+02hK1YfiUN\nWkIVAKxOXddt2zZN0zTN4XAoimKwSSmEod1u17xIC6c5LDze7/fH43Gz2XQ2HgPcczOjOgCszm63\ni49DZ19sW0r1W5iqqhqZRD60aaVrVVW12WwGN/58tFQBwOpcfnFfURR1XXfGXd3wXmvoARSqAGB1\nLr++73A4hKam7XZbluXN6UqoAgBWLaSotm3TdHVDQro8xi2XUAUAnBfS1cio9lPW0EYVCFUAwLAw\nQ1W6ZGRUe1EUm82m6KWoOJfVFDWcFaEKADjpeDzGcVRxXtD04sEimVIh/D/2D54q/6xMqQAADItX\n/6UTTe12u9j9dzgctttteDVMtZAu6Zd/buXIbBMAwIiyXNFptH/7mrzlZ+jav++KPg1FUbiTNgB5\nreo0ujbXhqrVdf/59AOQi3+rkzJQHQAgg9W1VAFARhqriIQqALidUSVP7NrErPsPACADoQoAIAOh\nCgAWrKqqGU4HNX6/vzgD+20Fpl79ZkIVAOS33+/L0y6fYXxZdyNumibs4Ha7HdzTcFjClOtlWfbj\nYFVVaYHO7oftX756593Prv5W7ZqsbX8BmNTIaSXc7W6z2eyGHA6Hy9/i8sIPF6JFqPDhcEifxiWb\nzSYsiYeos/put0tXj6+eXT3czjmuHp6mR2+weuO7c93uX1V66YQqADI6G6quykP9wuHEH1LCtcbf\neoqg1t/lGIPC036ICbknLdzZ2XT1tPDgO3Yy1rWr9wlVY7TbAZDR20NViBGhcJS2xKQLY+xIm2QG\nm3OiTsjobPPa2HdWf2udUNXJTGmQCnXrrJ4mof7u9FfvVKCzen/7/W12Xj310qDVjakaPxyPrh0A\nq1PXdV3XaZAKI4H2+314GroL4wCgMMwoXRI1TbPdbmMH2eFwOB6PcbKlpmnquo4ZIiSS7XabcV86\nVQpVDQsHB4el5S8ZPdbZfrrl9L0GCxRFETJWx/F4PPu+FzL5JwBMZb/fD46G7gygjv+q3+/3dV3H\n03xVVXVdd67vOx6Pp1oBwmZjhqiq6nA4bLfb/X6/3+/D8vjWVVW1bXv5kPkbhMSW1mrwaISXBsNN\nVVXH4zHXaP2pL5MUqgBgKsfjcTArpFGm03yy2WzG204Gm1vi23VeTZtqQkTbbrdpK9d0oSq0kHW6\nI6cwcrjCLk9dgWh13X8AcDenRiy9ZZvjzS2hvy+Vrhj6E+OEAhPN1xBmLihe91FO10o0kjL7MzJM\nVIdASxUAPJWRxqHYD7jf74/HYxiAlTdqhHFdxYm7IjZN0x9HFZYMNtGNdBrOkFAFAM9js9mcjSCx\njSqOWMqVWuJI+X5QOztcPVSm81JnyWDLU7r6SGjrby0Yaei6lu4/AJi1y1uSBht74qip0EDVfylj\nS9V401e/eqECMRUVr8d4he3EaSBOrd7/f5QOMgvbSes2uMpbaKkCgKmcuvqvuOZcHu5Vd0lj0n6/\nD+OlwmCm2BMXrh8MUyqkT/OmipHB7/GNttttVVXx6fF4jJmpqqrNZhNrWLy+eLBI7oET+jebpklX\nL4pit9vVdR2PeVg9pqhwZWUYp18kE0zk7Fu8alarpVvb/gIwqZHTSn+azcHzb9GbfLIz8Xcsn07+\n2alDcXryz+L1YPlJJ/+8JGx0qtefLL7TGdep3tnVB6dRPVXJkWk/Y+GL975t27Zss854+cYBZWdX\nv6TAyKtlmXl/b1OWH/sL2/bD/WsCwFvM5LQyKONIqdW69u+bbUxV577TnTtLD96juzOUrLP6+Pb7\nTYud22LP/Lbebfsh/e/R1QHg2UhU95cnVMVu2nSW/c7M92HW/FQajELh9MbRaa4K/aDpPRFDp2ks\nEOfDaKeZdx8A4IyrOgtHOh2LK28cPVL4ktXTdzx1W+x+V2uu/X2jovjh7BIA5m8mpxUmcu3fN09L\nVRw9F/VbHUfaIdOh/kF6PcLgOKr0wsh+gfD4njPTAwArlydUde71WLyeeSIVrgvtb2HkxtFn7zs9\n9VxeAABnTTL5Z5g6IsaaEH3quk5Hml87gG68vAgFADxW/lAVpxrrt0jFYU9hUtSH5KTB6xAvMUVl\nAICnkTlUhclSw6DyuLCqqpClYopqmmZwKv3U+Ks3e8toNQCAU3KGqjjxwSVzRI3cWHHQeMmzt2AE\nAObp1Hjr28qf3dq1b3e5nJN/1nW92+1uq+jIYPPB+JWOXj/VUWigFQCPEuasPjV6JN7GbuZTVU+t\nqqqrxlt3yneO3vg85P3VM95K+Yubu8NSp+aFCvrzSLXnJprqz1PVqWq6pP/u5qkC4A5GTivxJnSD\nJ8f4z/6Mt95bnHAQdrtdmJgpPj1VPj2ecSKn+Go44GGm8TY5/p3V4wEPBUberr0+NuSc/HPX00lF\nIwcu7urhcIjHJb6aHqm4evpBTLffP45psSz7+0ZCFcBzOBuqOndHTlfsn8uCcKYbedPnyGEhFY3M\n7D1Yvn8z6biF/rppA80l96Lue2So6ou7GpNQNH5j6v7062fvOz2+eqznG/c0C6EK4DmcDVWDfTXx\nBF8M9dKcOpd1JtnurBsj2vi5clb61RtJOeHgdBammbUfRpcaqi439R94fPtCFQAZnQ1Vg+0xaQDq\n9EbFkv1OmzQn9Tu/Ov2J/QKLMNIsMtjmd6ohcPDV/sZH3i4WOFPjTvmrSi/dTD5eQhXAczgbqtqh\nE3/xMl6lGB093FkxjhYKOqv3OxMHG8nm7FR/aDCYEUf2sR9nY9Dc7XbhvcYTVfuoe/8BAIPSu9n2\nnwbhIrVO91+nZNM0gzd/S59ee7eS+aiq6ng87na7LLsQ5iHvz/EUslRd12HOgeyHS6gCgAmFM3dd\n1+FpmNBxsGS4n1sUbk8SxVkYBl9dtJCoDodDljkO4jzknUQVjljacFXXdd5c9VXGbQEAfbvdrq7r\neI4/FR02m82pl+It4A6HQ5y+8TlyVZjKK+7XKYM3YkknrQzirJmdI9lvIKyq6uzNXa4lVAHAtMKZ\nPk7kfSo9VFV16qWQBtrTl9svVEhUl+xXaM3qLOwsOZWoiqH4Fbc52K96G91/ADC5zWYThvJ0Bk4F\nnS7CoGmake6wJ2imCnt9YVIMhdMD0hmI1jTNqURVXHB3liy0VAHA5MJAn+J031/oIqyqKpzpY+9e\nKB/aVPb7/X6/D2ErdF0t9y43TdMcj8fBHs+4pCzLODQq9NaF3BliUOd4xpTZ2WBo/wuNWNvtNvYz\n7vf7UIGce3XVtYJLN5P9NaUCwHMYOa30J6sshq7wv3zyz/T0HyaxjEvaoYkbZj6lwmCLXZBOMzFy\nEIqhe6v0jcxDnn1KhXKkHs+nLGexv2X5sW0/jC8BYP4mOq1kHOVD3+WH99q/7+q6/07dMDyYQ+QC\nYOUkqklNd3hXF6rEJgBgCq7+AwDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyECoAoClijdpnrmr\nKvn2nRpffbqDJlQBQH77/b48oX9zurB8fDuDC7fb7Xa7HdzsTMR6hgdn4066U51ZOgcPZlomrD7y\nduG4ndp+Blfd1GbpZrK/7v0H8BxGTivh3na73e6QiDe8S+86F29IN3ifvv75Omxks9nE8uE2gsXL\nPQHnI92vWMlThWOBUD7uZn9rHeNvFwukxy39W4zX/7r9var00glVAGR0NlSN5KT4NIaq/v19QzLo\n3Cz5VBQY3MID9Y9A2J1Tye/sPaHHd/Ds2/WPW/8dO66NDbr/AOCuwum/0zO12WyOx2On5FU9em3b\nzmp8VahM2sUWHtd1PVi+v/uh/IU7td/vD4dD/+3i6qGB6ubtX0KoAoC7GjyLh/zUeel4PMZeqii0\nr8xzBFWqH5KKl8pf6FToaV4Mlk+LpQurqhov8HZCFQDcT9M0IW10zuXhaRqVwuN+eApRoK7rMNR6\nv9/PqoEqdVWE6gv7lYaz4/EYR5qHweYjq2+32+J0+tzv94OZ9S2EKgCYSrwMLQpn+jYZgR7tdrs0\nQNR1fSqUhCFEocewruvwLvNvuxoX8k2aNcOx6ugMPD+Vq8LyweNcFMV+vw+HN+9BE6oAYCqbRFgS\nxk0PFu43TY30TIUGqjAcO8SLuq7zzxFwR/v9PsTEGEDj9XqhQMhScR9D+aLXPxgmViguSFT5W/iu\nGta+dDM5Gq7+A3gOIyeOwav/Bs816TVo8XFYvV9gRCg2eL3hQ5zd2UEhI4apKNpzV/z1LyeM10ue\nWiVU4MK5J64NBqtrqTp7+ABgIiEqjXQ5xeHqI31/p9pX5tZMNVj/4/E4PtAqjBLb7/dVVZ0dSN5v\no9putyNNUFVVHY/Hw+EwVVfpVRFs6Wayv1qqAJ7DyGnl1DxV4eSbLu9PQ9Vpc0oLjLTEXNigdTen\nmpFOtRKFUWLpknSP4mSqpwq0109kdda1x/OrSZIaADDkcDiEK9faE90ju90uzOQ02EJTVVUYeBRa\ndMLCML/A2UagO0tnpQrNTuFxrHYY27Tb7cKSqqrCsLA4tiy9Oi8UTo9MKBB3OSxMD0tcntakPxdD\nf6qF210VwZZuJvurpQrgOYycVkbaRTqNLv0WpuJ1k8ypqcY75naPmiDNeZ1mpH7jU2e/+nvUSY1p\ngfGo05n288Ljdm1sKMfr8WTKchb7W5Yf2/bD+BIA5u/hp5WmaeY2lGpQ3nreba+v/fvOImTczcM/\n/S/VEKoAnsFMTitM5Nq/7+qu/gMAmIJQBQCQgav/AOB247efY1WEKgC4kQFV49Y25kz3HwBABkIV\nAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABqubp2p8lrZVTacBAGS0ulAlNgEAU9D9BwCQgVAF\nAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQ\ngVAFAJCBUAUAkMFXj67AvZVlOfJq27Z3qwkA8ExWF6rEJgBgCrr/AAAyEKoAADIQqgAAMljdmKo1\nKMuP/YVt++H+NQGA9RCqnlMnQg3GLAAgI91/AAAZCFUAABkIVQAAGWQOVU3TNE1zc4E3rh4KjNUP\nAGAa2Qaq7/f7uq7TJYfDoaqq+LSqquPxGJ/udrv9fh+fNk2z3W7T1TtTn3e231m96N1/pvPuz8Fl\nfQAwW3laqpqmCYnncDi0bbvb7YqiSENSSFS73a5t28PhsNls6rpOW5VC4bD64XAoXoekkKg2m03b\ntm3bhtXTUBUKp6t3ItrTaNsP6X+Prg4A8KLNob+pkKtCygkFYiTqr9Ip/PbVQ64KGa6z1vU7l19R\n/HB2yVtWvHn7AJDRTE67d5OnpepwOIQcE6Vdb6FJqdNbt9ls4uPQynVqldCg1enLC0EqvNQvEB53\nuiMBAKaTJ1RVVdUJPSEPhYWDqSh9qXidsaIwBuvs6ulQrWhwgwAAE5lkSoWmaY7H43isuXYU+Xh5\nEQoAeKz8oSpex3fJ7AYjZSbKSeWtpqgMAPA0Moeq/X6/3W7DoPJLyo+0Pw126r3dzaPPpqgMAPA0\ncoaqOPFBp/1pMDldO0vnePl+ApsokwEADMoWqkKi2u12/fTTGZMedIafjww2P7v6qY5CA60AgLvJ\nOflnf5bzIESfzkvpSPZ0foQgnVJhcH6EdBaG/pQNgxcMAgBMp8wyWiiM4w7ZKBWnWog9gyH6xPnT\nY+6JU6IXLxEt7UbsrL7f74/HY3/13W5XVVWc3r2/a2WZZ3/fqCw/diZD7y95y4o3bx8AMprJafd+\nskwhemrj6TTonciVToAepL11nfnT3756rOdbdjMXM6oDsAYzOe3ezbNFyKZpRnr9ZhKZtVQBsAYz\nOe3ezSSTfz6QcVQAwEM8W6gCAHgIoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIIOv\nHl2Bewt3CTxlVRO/AgAZrS5UiU0AwBR0/wEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEA\nZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGTw1aMrcG9lWY682rbt2zb+\ncWibH96yTQBgEVYXqt4Ymy7Y/qsINRizAIDno/sPACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAg\nA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACCDrx5dAYaV5cdHVwEAuIJQNV9t++HR\nVQAALrW6UFWW5cirbdverSYAwDNZXagSmwCAKRioDgCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUA\nkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQwVePrsC9\nlWU58mrbtnerCQDwTFYXqsQmAGAKuv8AADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQ\nqgAAMphXqGqapmmaNxbIWiMAgItkDlVlWZZlud/vB5d3VFUVCzRNU5bldrvdbrfh1c4W9vt9WuDU\nW8TVpSsA4J6yhaoQekYKbDabw2tpMNput0VRHA6Htm0Ph0Px+iZ9+/2+ruvNZtO2bdu2m82mrut0\n9VA4XT1sEADgPvKEqqZp6rre7XYh0JxS9YTlIR4dDoewpKqq3W5XJH15IVHFp+FBXdenVg/V6Ldm\nAQBMJE+oCjlmPMSknX0dIR6lBcKmwv9DhOqsnqaufoHwOKYuAICpZev+G8lMqVMjzTebTX/h8Xgs\nToSq8DS8FIpdskEAgInc4+q/2FuXjjS/MIRF4+VFKADgse46pUIYSB5Gmh+Px4fkpMHrEC8xRWUA\ngKfx1R3eo6qqtm3TJU3TVFU12G0Xjb96s05NAACyeNjkn+mgqEuMl+wnsIkyGQDAoLnMqD4y2Hww\nfqWj1091FBpoBQDczT1CVZgXtJOK0mkUOrNSFa+nVBicHyFdPS0cDF4wCAAwnTxjqjoTJTRNE/NQ\nVVVhPvTtdrvb7dIYFLJU8TJh+na7DZN2hqlEN5tNTEW73a6u67CpIpntM7yapq6qqsLqhck/AYB7\nanOI8agj3ljmcDh0OuN2u11nI2mBuOKpt4gXEl64epBrf08pih/OLrmw2OCKl2zqLVsDgIymPu3O\nTdne92q4cN3fo7ZfltPub1l+bNsP40suLDa44iXv+JatAUBGU5925+beA9WnHudkHBUA8BBzufoP\nAGDRhCoAgAyEKgCADIQqAIAMhCoAgAyEKgCADIQqAIAMhCoAgAzy3PtvQcqyHHl1VRO/AgAZrS5U\niU0AwBR0/wEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBU\nAQBkIFQBAGSwuhsqE5Tlx8HlbfvhzjUBgOcgVK1XPz+dSloAwFm6/wAAMlhdS1VZliOvtm17t5oA\nAM9kdaFKbAIApqD7DwAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOh\nCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIIOvHl2BeyvLcuTVtm3vVhMA4JmsLlSJTQDAFHT/\nAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBksLp5qmYinYLUzFkA8ASEqvtJgtSH3vIP\nZSldAcCCCVX3MHhrnBCh0pfCY9EKAJbImKrJdXr62rYoio8xOcUlafnR+xMCAHOkpWpqX3r6zrY/\n9RquPpwsCgDMj5aqCcWEdHmP3kvD1avVAYD5E6qmckOiitJcJVoBwCIIVZNIktDHkWKjXo2yAgBm\nbnVjqsrRhNLmuPQubaN6Sx5KR1mZcAEAZm51oSpLbBrxll6/QW9MZgDAfawuVN3HDYmqLM90FGqs\nAoA5E6ryunEehLYdWzE2VslVADBbBqrnN0XuMc8CAMycUJXN1HFHGxUAzJlQldkdoo/GKgCYIaEq\nj/sEHZ2AADBbQlVeN0/1eSm5CgDmSajK4M75xuAqAJghoSqb+2cdjVUAMB/3DlVN0zRNc9urFxa4\nrWI3e0iy0VgFAHOTOVSVZVmW5X6/779UVVVZltvtdrvd9ss0TZO+2r9D336/H1k9vnVc/T7pKvtN\naa5343SjAEBe2UJVCD2nXq2q6ng87na7tm0Ph8Nms6nrOs092+22KIrD4RAKFK/vfLzf7+u63mw2\nbdu2bRtWT3NVKJyuHjZ4Hw9JVEasA8C8tDmEHLPb7eKDToGiKGIkikviu+92u+IlEg0uuXb1kWrc\nsHenFMWX/5IlP/TKdJecWviWOpzd/v0rBgB5T7vzl6elqqqqw+Ew2OtXFEVY3nl1s9nEx3Vdh40M\nrhIatNJXi6IIQSq81C8QHofNTm0Ow5s0VgHAw2Xr/uuEntRgKgpPYw9gmrGi4/F4yeqhWMfgBjOa\nTY6ZfGYsAOASD5tSYSSE3VB+6gh1yhyaqYLZhDwAWKkHz1M1co3eRDmpvNXrjUxRtRsZsQ4Ac/DV\nY99+pP1psFPv7dp8jUvzaaZqW4kKAB7sHi1Vg8np2nmkxsv3E9hEmayYd4PQnOsGAM/tfqGqk4o6\nw89HBpufXf1UR+GkA63m00wVzK0+ALA29wtVnSkVjsdjDD3p/AhBOqXC4PwI6SwM/SkbBi8YzGL+\nTUHzryEAPKU8Y6o6t+RrmibmoZBsdrtdXddVVYXlYbrzGIPChOnb7TZM2tk0TZg/Paaizurh/6Fw\n8Tp1VVUVVi96MS6jeTYLGVkFAI+UZQrR0NTUl06D3imTToAepL11nfnT37568Mb97U+h3itw/xnV\nf3j99PY51vNWDAByxYylKNt5trrcqmmakV6/snzT/oZ2oJENlOXHtv0wvuTUwlur1H/HLw/Set6/\nYgDwxtPu4jx4nqrsphhHtSxr+vQCwIw8W6iazuKGKy2uwgCwaELVdRbSDuSGgABwb0LVM9NYBQB3\nI1RdZHHpZCEtagDwPISqKywxqSwuDgLAQglVT2uJERAAlkuoOm/pjT1Lrz8ALIJQdaklNvwssc4A\nsFB57v23IOVou83zTvxqnnQAmNbqQtUlsaks03melh1H3GUZAO5jdaHqQvEWeC+J5OPS01VZ6g0E\ngAkZU/X8ZCkAuAOhasyTdZw92e4AwKwIVec9RUuPuwECwLSEqnXRWAUAExGqTnqy/PEU7W0AMF9C\n1RmyCABwCVMq3O71dFaLYW4FAJiCUDXswr6/OJ3VIpgIFACmo/tvzLO26IhWAJCdULUuzxoTAeDh\nhKrVWlLHJQDMn1A16JkDh8YqAJiCUHWS8AEAXE6oWi/D1QEgo9VNqVCORom2bdcQNcytAADZrS5U\ntZf16q2k789EoACQi+6/1VrkdPAAMFtC1drpBwSALISqV1aVMHT8AUBGQtUgXWMAwHWEKtbVPgcA\nExGqfrXCbKEHEAByEaq61pkzVhgoASAvoWrt1hkiASA7oYovNFYBwFsIVZnQnBMAABD+SURBVF+s\nOVJorAKAtxOqXll9vPjw6AoAwFKt7t5/z6csM8yq5RbLAPBGQlVRLLnvr20zty25xTIA3Eb3369W\nHiZWvvsA8EZCFV3LbbcDgAdaXagqB/z60kOrNgdueggAN1rdmKq218sVolTbFkWhA+wLI6sA4Fqr\na6linCwFALdZe6jS49cRJ2goy6IsP2aZrwEA1mDtoSrQPBO07Ye2/RCPRvb5GgDgiQlVnKQZDwAu\nt+pQJTScoukOAK616lAVCBCj9AACwEWEKobJmgBwldXNU8W1OnNWnboe0Kh2AFZuvaHKgKqz2nb4\nKPXzk5kXAGDt3X86uS4hgALAWWsPVZyjCQoALrLSUKXp5VqOGACMW2moCvT9XcJRAoBLrDpUcRWN\nVQAwYo1X/60zHLzlAr1TlwECANHqQlVZlkXR/vrwtfZJ+7pyTSIlWgHAKasLVW3bhmTwrPlpIhqr\nAGDc6kJVYeT1W5k5HQAGGKjOpYRRABghVHE1/YAA0CdUcQWNVQBwyhrHVN1br2GnLYqi/CZ5vryo\nUpZLrDUATEiomkxZFi+TN1xSsiiWkq4+GqsOAH26/yZQlt3WqbZN/yuLH359Or7ijC2npgBwD0JV\nVp1UFPPTiH66mn20WkaDGgDc1+xCVdM0TdPc9uolBSb0Ok5dHT2WFq0KjVUAkLhTqNrv9+WQtExV\nVWVZbrfb7XZbluV+v09fbZomfbV/i5nwFqdWn1YagG6IU6mFRCuNVQDQcdeWqt1ud3gtvlRV1fF4\n3O12bdseDofNZlPXddrmtN1ui6I4HA6hQPH61n37/b6u681m07Zt27Zh9by5qiw/dv6LL/xaKFfW\n6EeruZpx1QDgvtq72Gw24+9VFEWMRHFJXGW32xUviWpwyfjq6cJLalsUP1xSLBT98t81m7pi+xe8\nxaT6VU2XxHpdt0cArMPdYsZMzGJMVWhS6jQshRwW1HVdFEVVVYOrhAat9NWiKELqmnZ81X2mQogb\nn2tX4CwrBQD3du9Q1bzoLCx6qSg8jSXTjBUdj8cLV8+sM4jqTNkT/YaXm2tXoJFVABDdafLPkH46\no8sPh0MnCaWqqgoNVBca2VRm1wyiattM82SGNwpvPcfpzE0HCsDa3Xugeuh0DCPNt9vt2ZakkQKD\nbVdnDV6E2Lsm8ZuTFyrmusrvNmlX4DzMLt0BwIPcKVSFLBVHTVVVFXLV2Qv0RtqfQuvXbTUZd2J0\n+TzuJzO/XBXMrDoAcG8PG6iepqXB5HTtcKjJ5/ycQ6LqVGAeQebhxwMA5mAWV/8NDirvDD8fbJcK\nPYCXrP5W80lUwR1z1eWj7OeR8QDgQTJNzXC2Q607j1SY8iCOsuoXSJdcMk9VZ18G9+7C/R3o/nvc\nTFFjHjqF1euK/DCPigAwI3eLGTNx18k/N5tNmEg9RKL0WIclsUB4NU1RcUlcPQ1hndXD26Wrx41c\nUttuqJpNdhkwj7rFUPXoigAwI0LVVGKQCjrtUv0C/UiUXu53w+rtbaFq/mFhBjUMR2wGFQFgRtYW\nqsp2JoOE7qIsL9rfsvz4ZX6puQ2lOuXR9YxHLFRk5kcLgPu48LT7NGYxUH3u5v+BmNn1gPOoBQDc\nlVB12rKiwTxy1fzzJwBMRKg64dEdareYR66aTRUA4K6EqlELSlTBDHLV4o4ZAGQhVA1oi28eXYU3\nmEGumsf7A8BdCVWnLbfJ5dG5arlHDgBuJlT1PEcDy6Nz1QzeHADuSqh6bYnj0095aK56guMHAFcR\nqgaUxQ+PrkImM2iv0lgFwEoIVYmnPP8/LlfNINEBwP2sLlSVp30p8XwdVzPIVQDw9L56dAXu7eRN\niJ67OaVtv+xgWd4n6ZTlx5eH8Z6AL3dUBIBntLpQNeyJm6miu+eqGKFejq5EBcAzW13335gnTlTB\ng/oBn/64AkAhVBXFs3f8dTx09PiqjjQAa7P6ULWGjr+OR+zpeo4uAKu1+lAVrPOcr7EKAPJZd6ha\n7Rn+MZ2AHx/wngBwL+sOVcE6m6kekavWeaQBWIkVhyoNJo8btO7YA/B81hqqVjg+fdDdc5V71wDw\nrNYaqoKVJ6rgcbkKAJ7JKkOVRpKOBzUf+TsA8ExWGaoCDSap++YqnYAAPJ9Vhqq2lagGPChXAcBz\nWN0NlcvyY2dJvO8v97/p8n3fCgAmtLpQ1YlQ/Yy1dnfMVY+IcAAwldWFqrLXt1WW38THrXN78Zhc\nBQBLt7oxVe1rRfHD66cURfGAkeSiFQBLt7pQxaXuFTFdCQjAcxCqOGf6sKOJEIAnIFRxmk5AALiY\nUMWoe+UqnYAALJ1QxTlyFQBcQKjiAnfPVQCwOEIVl0lz1fRNSRqrAFgcoYqLpe1Ik6UenYAALJRQ\nxTXSe1HLVQCQEKq4nlwFAD1CFTeRqwDgNaGKW8lVAJAQqniDl9TTFt9M/A5yFQBzJ1TxNtqrAKAo\nCqGKDOQqAFhhqCpfK4pvXj/lFmXxw8sjuQqAlVpdqGpfK4ofXj/lVtqrAFi31YUqJpSknomGrstV\nAMyWUEVW09/KRq4CYJ6+enQFeBJl+fHl4Q9FnGShLIsJOlXb9kuimmbzAHALoYoM2vZDb9mHSYOP\nXAXA3Oj+YzJpR90EfXX6AQGYFaGKKSWNSG3xTfb4M3FsA4Ar6P5jYiH4xMgTHrykoWQkVrpGvzNx\nbPPptnUFAvAoQhV30bZl+fHXeRaS+NOJUIMx69y2v2yykKsAeBzdf9xR207XY2eIFQCPJVRxP2X5\nsSw//npPmzjzQiaGWAHwQLr/uJPX3XwfiiLpsSuKXJ12nSFW+TYMAGdoqeJxppl+Pe1jLLRaAXAv\nQhUPNdkoq360AoBJCVXMQGc6q6wbNtAKgPswpop56MyLUBTpePbiysmr+ts2lxUAU1tdqCp7jRVl\n+WvTSOt8+1hJ/PnSZNW2xU2TV/U3XEwyMh4Avlhd91/7WlH88PopjzblOHMD2AGYzupCFcvwOlpl\nvG+gawMBmMjquv9YksH7BhYZuu5Gb0gIALfQUrUk/QFhC3J75du2M2j9S/vSm4/GYKvVqa0u+uAX\n6v9Qi658of4PtejKr5CWKhZi8N5+OdquOq1WmbYKwOoIVSzNNOlqsswGwFoIVSxSWX4sXvoEX80X\nOmG6EqwAGCNUsVTJdKAfiqIoy4/D6aq4MWD105Xx7ACMEKoyTCzJXAy2MhVvbb46dQ3irdsD4DmV\nTzbjZdM0RVFUVTX4allet7/KZyz/lo0PBt/OjWtObn/w2pm2va0+4xfivL6WcEYHf23lZ1UZ5ZdV\nflaVWWH5pXuelqqmabbbbXy62WxCwOIJvOXGf6cGSbXpkou/850pGDoMbAdYsyeZpyomqsPh0Lbt\nbrc7Ho+n2qtYqTAt1anWrP5/F2+vv9UwsD3TdFoALMOTtFTt9/siuR1yeFrXddM0otWaneg3DK1U\nL43Sp6f7fLXWufc6NaBr5B06G7ikEACz9SSdnWVZ9vv7yrLc7XYhYMUls+o8XlX5O1Sm6Ey8XhTF\nwNCrj2HJme2/pX3pdWh785aGzeqPO7fys6qM8ssqP6vKrLD80j1JS1VxYnC6YVWrcuHQq9h8FR8M\nrNj/Fbg8HMVZra6c3qrslR19z7FR82v6EQOYi2cIVSPJ6Xg83rEiLEDMT2X5zUuT1ccLp9W46N9b\nb2iYaot8w696W7q20WzR5WdVGeWXVX5WlVl2+VX+2+4ZQtVVrr05pfIZy8+qMncof5U1/vwAz2ud\nt4JeV6haVc8uADzKOk+3zzClQhhNNdgJuNls7lwZAGCdniFUDRqfWh0AIK8nudaxqqrj8ZjuS38J\nAMB0nmRMVdM0ZVmWZXk4HMLT4/HY6ftbetvV0uu/dIs+/ouufN+id2fRlX8Ciz7+M6/82erNvP55\ntM8ixKlos9lc8tIi7Ha7zl8t3I1nWZZb8046X9YuLPHDE+q52+36L3X+FoNlHuvyys/zh+iSAxs+\nVIur//xPBMv98Jz9nZn/NzeX5wlVp8QvUrwt4Aw/kSMG618sLQ3Hb9T8z+gdoebhJyD+LR5dqUt1\nPu3zr3/609z/2e38LcLT+XyiLqz84XCIlZ/VqWW8/qn5n9f79Q8f/s1mk/6Qzuf4j1e+8+M5t8qf\nPUnN/Jub13x/XnMJf790SfiTL+Uv2j8LLqv+bfJztqxqty81T3+8DofDfH7LzlrWhyce7f5hD/on\n8vlkxPHKx69AunA+lW8vOPhR/AfSrELV2foPnghm8l2+5MPTWTirD8/Z35mZf/jzetqr/6L+4Kpw\nN8Cl3MEm/Ls2XbK4DuntdrvZbBZX7eLlo5LePrKqqvTp4sz5r1BV1eFwOHV4+3+LYk5zpoxXPrw6\n59+c8fpH+/3+eDx2fpHm4Gz9j8djp4tqv9/P5Lt84cGfrfGT1My/udk9f6gqFn5bwKqqOvUPn845\nnx1TI7OIzV9M5E3TzOcn+HLhLJIe/Jl/eK4d4jqrT9f4UR18dVanlks+FXVd73a7eX5+LvzwhO/y\nTD4z0Ujlw0t1XXeWz+fDM36Smv83N68nD1XPd1vAwQsbZ2u2/669SlVV2+22ruu6rsuyXFC02u/3\nm81mu92GBrYwz8jS/xypeZ7dLxFqvqDPUvFy15Fl1TkIJ4JwkXj4Lm+32wXdRCV8Z8OPT1VVZVlu\nNpvZJpJLTlLL/eae9eSh6sk0TbPdbovlBPw5/7v2QiF8x/7yoijqul7K8S9efryOx2Nd1wuK41dZ\n0J8jCOl2WV+NkKXaJc/8V9d1HOUTY8qjK3WRqqrCNzd8i4sZh5KrTlKL++ZeQqhajP1+HwYnLeV3\nbYn/Fh+UfvPDb/FSdqqqqvREEj452+32yX7LZnuCGRTbC5fyKSqKomma8A+kR1fkTdIUG2PKIpRl\nGeayjomwrusZfuyvPUnNcBfe7slD1UjH7YK+UUVR7Pf7uq7n3OTbETr+drtdkyiKIj5YqGX9CoSm\nqbTO4eAv6HQeLXpkZBTOjofDYVkfpND2UFVV+l0+Ho9LOf6DR3spw3rCtzXtta+qarfbzW0Ey6mT\n1HN8cy/35KFq0OImdQ0f1hBQHl2XS4WqhrELQfgJCEseXLmLLSt5P7fBU+Cyvsuhs6lt26VUuGOb\niEsWui8LMvIhn88ZYeQk9QTf3Ks8f6jabDadRD94hedsxYb3pVQ4aJqmM3tHOnPJo2t3qf7Pwcyv\nnuvrfPjn8yt8rcHe5AWNEgv1X9CHP9WfjKd4mXloEZ+owQvowtP5f5cHP/mzCiXjJ6mlf3Ov9ST3\n/hvRXHBbwDmL/yjsfCj7V7GSXfjn13a7DaMxwm9HsZxEHv5FES/9i2NI51n/TtdwuPS9SD7qu90u\nDCUJy2e1L+OVjz87/douov4PqtQVztb/cDjEy2CLl8M+k1Fi45UPv0Lhkx8DyqzOYmdPUnP+5ub3\nxslDF2H+t3waceoPt6y9aOc9l/e49MdroYc9Nds/wakzXHrMO2Xmsy/jlR85ec9kFy45+Km5fRcu\nqX/nRDCT6dTbCyofb+1y9u/yEKc+24v45mZXtstsjoY7a5pmEf9kP2Xp9QeYP6EKACCD5x+oDgBw\nB0IVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZC\nFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUA\nQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAGQhUAQAZCFQBABkIVAEAG\nQhUAQAZCFQBABv8fag3j5DOpoCIAAAAASUVORK5CYII=\n",
       "prompt_number": 7,
       "text": "<ROOT.TCanvas object (\"icanvas\") at 0x419f000>"
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "",
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}