docs/source/docs-factorybinned.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Build workspace using histogram templates\n",
    "\n",
    "Build a binned likelihood model version of the ex11 example\n",
    "\n",
    "* construct a histogram template SH(mgg) with a prediction for the binned signal shape\n",
    "* construct a histogram template BH(mgg) with a prediction for the binned background shape\n",
    "* construct a probability model model(mgg) = SH(mgg) + BH(mgg)\n",
    "\n",
    "This model can be 'seen' in two ways\n",
    "\n",
    "1. an extended probability model like ex11 that happen to have binned shapes, i.e.\n",
    "\n",
    "  * model(m<sub>γγ</sub>) = N<sub>sig</sub>/N<sub>sig</sub>+N<sub>bkg</sub> * pdf<sub>SH</sub>(m<sub>γγ</sub>) +  N<sub>bkg</sub>/N<sub>sig</sub>+N<sub>bkg</sub> * pdf<sub>BH</sub>(m<sub>γγ</sub>)\n",
    "  * $P(N) = N_\\mathrm{sig} + N_\\mathrm{bkg}$\n",
    "  * where pdf_SH,pdf_BH are probability density functions that follow shape of the unit normalized histograms\n",
    "\n",
    "2. A product of Poisson counting experiments for each bin\n",
    "\n",
    "  * model(vec_N) = Product(i=0..n-1) Poisson(N_i | S_i + B_i)\n",
    "\n",
    "  * where N_i i=0...N-1 == vec_N are the observed event counts in each bin and S_i and B_i are the predicted signal and \n",
    "background rate in each bin\n",
    "\n",
    "Both representations are mathematically equivalent, but expression (2) is in practice faster to calculate because it does not require a normalization calculation over pdf_SH and pdf_BH to happen. While for this very simple example it does not make a noticable difference because the normalization does not depend on any model parameters in scenarios where it does it will effectivelt double the calculation time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construct simulation workspace to generate template histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\u001b[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby\u001b[0m \n",
      "                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University\n",
      "                All rights reserved, please read http://roofit.sourceforge.net/license.txt\n",
      "\n"
     ]
    }
   ],
   "source": [
    "RooWorkspace wsim(\"wsim\") ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate two distributions, exponential distribution for the background, Gaussian distribution for the signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wsim.factory(\"Exponential::bkg(mgg[40,400],alpha[-0.01,-10,0])\") ; \n",
    "wsim.factory(\"Gaussian::sig(mgg,mean[125,80,400],width[3,1,10])\") ;\n",
    "\n",
    "RooDataHist* hist_sig = wsim.pdf(\"sig\")->generateBinned(*wsim.var(\"mgg\"),50) ;\n",
    "RooDataHist* hist_bkg = wsim.pdf(\"bkg\")->generateBinned(*wsim.var(\"mgg\"),10000) ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mock data distribution with mu=1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wsim.factory(\"expr::S('mu*Snom',mu[1.5],Snom[50])\") ;\n",
    "wsim.factory(\"SUM::model(S*sig,Bnom[10000]*bkg)\") ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that the sum is an extended model, no speficiation of the event count is needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RooDataHist* hist_data = wsim.pdf(\"model\")->generateBinned(*wsim.var(\"mgg\")) ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up binned likelihood model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RooWorkspace w(\"w\") ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First Import template and mock data histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing dataset genData\n",
      "[#1] INFO:ObjectHandling -- RooWorkSpace::import(w) changing name of dataset from  genData to template_sig\n",
      "[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::mgg\n",
      "[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing dataset genData\n",
      "[#1] INFO:ObjectHandling -- RooWorkSpace::import(w) changing name of dataset from  genData to template_bkg\n",
      "[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing dataset genData\n",
      "[#1] INFO:ObjectHandling -- RooWorkSpace::import(w) changing name of dataset from  genData to observed_data\n"
     ]
    }
   ],
   "source": [
    "w.import(*hist_sig,RooFit::Rename(\"template_sig\")) ;\n",
    "w.import(*hist_bkg,RooFit::Rename(\"template_bkg\")) ;\n",
    "w.import(*hist_data,RooFit::Rename(\"observed_data\")) ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now build signal and background models \n",
    "\n",
    "Note that we build _functions_ here and not pdfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing dataset template_sig\n",
      "[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing dataset template_bkg\n"
     ]
    }
   ],
   "source": [
    "w.factory(\"HistFunc::sig(mgg,template_sig)\") ;\n",
    "w.factory(\"HistFunc::bkg(mgg,template_bkg)\") ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also not that we don't need to declare mgg here, its definition was imported when we imported the histograms\n",
    "\n",
    "\n",
    "\n",
    "Now construct a 'amplitude sum' probability model, defined asa\n",
    "$$        \n",
    "pdf(x) =\\frac{c_\\mathrm{sig} * \\mathrm{hist}_\\mathrm{sig}(x) +  c_\\mathrm{bkg} * \\mathrm{hist}_\\mathrm{bkg}(x)}{c_\\mathrm{sig} * \\mathrm{SUM}(\\mathrm{hist}_\\mathrm{sig}(x)) +  c_\\mathrm{bkg} * \\mathrm{SUM}(\\mathrm{hist}_\\mathrm{bkg}(x))}\n",
    "$$\n",
    "here $c_\\mathrm{sig}$ and $c_\\mathrm{bkg}$ are coefficients scaling the predictions of the template histograms. If the template histograms encode the nominal event yield, one expects both coefficients to fit to 1 if the data matches the prediction\n",
    " \n",
    "**NOTE**: *If the bin width is not equal to one then the event count of a histogram is not identical to the integral over a histogram. In RooFit the __integral__ over the histogram is taken as the yield prediction, whereas one usually, interprets the histogram event count as the prediction. The simplest way to correct for this is to multiply $c_\\mathrm{i}$ with a constant which is 1/binwidth*\n",
    " \n",
    "In this case we choose $c_\\mathrm{sig} = \\mu$ (as usual) and introduce a Bscale as a nuisance parameter that can freely scale the background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[#1] INFO:ObjectHandling -- RooWorkSpace::import(w) Recycling existing object binw created with identical factory specification\n"
     ]
    }
   ],
   "source": [
    "w.factory(\"binw[0.277]\") ; // == 1/(400-30)\n",
    "w.factory(\"expr::S('mu*binw',mu[1,-1,6],binw[0.277])\") ;\n",
    "w.factory(\"expr::B('Bscale*binw',Bscale[0,6],binw)\") ;\n",
    "w.factory(\"ASUM::model(S*sig,B*bkg)\") ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the binned probability model to the binned data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[#1] INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.\n",
      "[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_model_FOR_OBS_mgg with 0 entries\n",
      "[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization\n",
      "[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (sig,bkg)\n",
      " **********\n",
      " **    1 **SET PRINT           1\n",
      " **********\n",
      " **********\n",
      " **    2 **SET NOGRAD\n",
      " **********\n",
      " PARAMETER DEFINITIONS:\n",
      "    NO.   NAME         VALUE      STEP SIZE      LIMITS\n",
      "     1 Bscale       3.00000e+00  6.00000e-01    0.00000e+00  6.00000e+00\n",
      "     2 mu           1.00000e+00  7.00000e-01   -1.00000e+00  6.00000e+00\n",
      " **********\n",
      " **    3 **SET ERR         0.5\n",
      " **********\n",
      " **********\n",
      " **    4 **SET PRINT           1\n",
      " **********\n",
      " **********\n",
      " **    5 **SET STR           1\n",
      " **********\n",
      " NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY\n",
      " **********\n",
      " **    6 **MIGRAD        1000           1\n",
      " **********\n",
      " FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.\n",
      " START MIGRAD MINIMIZATION.  STRATEGY  1.  CONVERGENCE WHEN EDM .LT. 1.00e-03\n",
      " FCN=-18731 FROM MIGRAD    STATUS=INITIATE        6 CALLS           7 TOTAL\n",
      "                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       \n",
      "  EXT PARAMETER               CURRENT GUESS       STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  Bscale       3.00000e+00   6.00000e-01   2.01358e-01   1.98597e+04\n",
      "   2  mu           1.00000e+00   7.00000e-01   2.24553e-01   9.86124e+01\n",
      "                               ERR DEF= 0.5\n",
      " MIGRAD MINIMIZATION HAS CONVERGED.\n",
      " MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.\n",
      " COVARIANCE MATRIX CALCULATED SUCCESSFULLY\n",
      " FCN=-27643.4 FROM MIGRAD    STATUS=CONVERGED      63 CALLS          64 TOTAL\n",
      "                     EDM=1.90326e-05    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  Bscale       1.00258e+00   1.02810e-02   5.15767e-04  -9.00358e-01\n",
      "   2  mu           1.54130e+00   4.86712e-01   1.62486e-02  -1.71092e-02\n",
      "                               ERR DEF= 0.5\n",
      " EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  2    ERR DEF=0.5\n",
      "  1.057e-04 -1.035e-03 \n",
      " -1.035e-03  2.386e-01 \n",
      " PARAMETER  CORRELATION COEFFICIENTS  \n",
      "       NO.  GLOBAL      1      2\n",
      "        1  0.20612   1.000 -0.206\n",
      "        2  0.20612  -0.206  1.000\n",
      " **********\n",
      " **    7 **SET ERR         0.5\n",
      " **********\n",
      " **********\n",
      " **    8 **SET PRINT           1\n",
      " **********\n",
      " **********\n",
      " **    9 **HESSE        1000\n",
      " **********\n",
      " COVARIANCE MATRIX CALCULATED SUCCESSFULLY\n",
      " FCN=-27643.4 FROM HESSE     STATUS=OK             10 CALLS          74 TOTAL\n",
      "                     EDM=1.90282e-05    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                INTERNAL      INTERNAL  \n",
      "  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   \n",
      "   1  Bscale       1.00258e+00   1.02822e-02   1.03153e-04  -7.28574e-01\n",
      "   2  mu           1.54130e+00   4.86745e-01   3.24972e-03  -2.77461e-01\n",
      "                               ERR DEF= 0.5\n",
      " EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  2    ERR DEF=0.5\n",
      "  1.057e-04 -1.038e-03 \n",
      " -1.038e-03  2.386e-01 \n",
      " PARAMETER  CORRELATION COEFFICIENTS  \n",
      "       NO.  GLOBAL      1      2\n",
      "        1  0.20664   1.000 -0.207\n",
      "        2  0.20664  -0.207  1.000\n",
      "[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization\n",
      "[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)\n",
      "[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAHYCAIAAAApvgy/AAAABmJLR0QAAAAAAAD5Q7t/AAAgAElE\nQVR4nO3dP4wjyXn4/WpDuJGDvSYDAQ4444N1yZ0BAac3sZywme2cTqkDATpnuzPDxNn8SdhMNHOJ\nsxkuJvxJjhQZknadsRm9AgzLFmBLhu3fYr2k4cAGhn2JdQL89hvUsLbYXdVs/u3u4veDwS6n2WxW\nkzPsZ6qeespLkkQAAACY/F7ZDQAAANVFoAAAAKwIFAAAgBWBAgAAsCJQwF4LgsDzPM/zoigq+JAo\nijyTIAiCIAjDcIvNFUIIYXv2MAxTZyGbuoMmLSRf55yWRFEkGx/N7KpphaSaV3JrgB1LgD2mfhF6\nvV7BhwyHw4W/VsPhcAdtNtJPRDa13W4v+xTD4XCDp6DaZmuJekmTJGm320u9HTtQ8eYB20aPAvaX\n/gduv99f9uGp3yV5cZV3dTqdTTWy+LPLa1i/31//T94wDDd+CkmS8Lc4UEcECthfMjjo9Xry2zUv\nY3LoIZn9Ab3ji2IQBFEUyVihCmMNinwdZMOKCIJge41ZR/FTABxDoIA9pS7kYRhu4/pqDBQKDsCv\nPE4vT2E0GhXZOecp5BGWyttYZ/BeBgfyXZC3jeFCFEXZPAxRIG+gYNtsu+lNymke4KzdjHAAVSM7\nEuSouT4IvdDCnY3D2NnMBuOAvereyNkt59lTbTPmKOS3JNsA22na9tdPvMgpqz3VA/UnlS+mGlhR\nZApF6vjZvIrso9QBc16Q7FusHiLvyn9NAMfwE489lbq02K40WQsDheyh1EPa7Xav19MvrvoD1Xa5\nm7rIpa6vOc+uHqI/r/7whS0ZDoepZuS8FHoL9QarR6U25hwt9XJln0IeIRsxyI3G11PfJxXQqKfI\neUH0o9maB+wDAgXso+zFXl5+ikwQsAUKejJj6opYcGPqGqY/l77RFiioi5w6ZjZQKNgS+WpYXoB0\n2/I3rjzzQm+JLQLIblQvlN5jlH2g2s3Y/WM8PrC3+E3APspeHoqPPiycHmkbdLAdR35rvLAZt+c/\ne7b/QG0p2JKkWKBgmyiY2r6RQCH1LDkbU11EqaOlAi/bC6JCrtXaDDiGZEbsI5msp6ekqdvrz1ZI\nTVCUt7Nj+amEOLlbNqFS7lYkP7Hdbg+Hw5z2F2xJQbJJtgZvdtKHsYUrNNv4kOx0hkpNGwFKR6CA\nvaOuYZ1ORy9uKDcWv0gYQ2951dGLEBS8ZGZjF8l2Ocw+dRRF+dfO3czYrMKMgIJzMuVuVWgwUGUE\nCtg7+aFAwbmFNraLcS2uRhuMJNZ8GQFUB4EC9o68hhlz1+UOG+l5Vhfdgl3x8s9fY5GA9Rujt2Sd\nHYpYtsLS9hSvJ7HlhgD1RqCA/aKCgJzh6m1cOQpGALbdNnLdlaecLVa92vnaqlRVoT/f9uypM7Xl\nf5CjAMzZbG4kUHH5i/oUmfuwQh2F7JYkMynRdljb0XJamGpqdsZEqiXZ1yRblch2cLGT6ZHGBudv\nNM53VfFWanIE0yOBHPQoYL/YcvUl9Zfo+n9T6n+8yokGnU5HViAOwzCbOxkEgbyMyeWY9d3a7fam\n/kBXLZFrQ8insL0m+ZWPVZP0BssszuzEit2TLR+NRnIB7jAMgyDIdh7IOKbf7+svyO5bC1Ra2ZEK\nsDtF5scbK/zoFvYoGJ+lSG3mgrsV/80tWMI5ezRb7cis7JiIsYzE7nsUEtOZ5veC5O8G7C0vWVS/\nBcCmqL/R5VKT+bvl77PtlqgdFnauFDypsqRORHYYZD/3Cu4G7CECBQAOsk02SUUAxt2iKOp0Ou12\nmwkRgGDWAwBXyQQFfYv8Vh9YGY1Go9Eo1WsiMy2Y+wBI9CgAcJNKS5TpCyqTUf/Qk50H2d3oTgAU\nehQAuCmZZTjKbgMxS7TU9wmCQC1dre9GlAAo9CgAAAArehQAAIAVgQIAALAiUAAAAFYECgAAwIpA\nAQAAWBEoAAAAKwIFAABgRaAAAACsCBQAAIAVgQIAALAiUAAAAFYECgAAwIpAAQAAWBEoAAAAKwIF\nAABgRaAAAACsCBQAAIAVgQIAALAiUAAAAFZfW+1hURTJf4MgCIJggw1ak+d5ZTcBALDXkiQpuwmb\n5C11PlEUhWE4Go2yd/V6vTAMN9auVXnecmcEAMAGuXcZKjr0EEWR53mdTicMwyRjOBzKHaoQKwAA\ngE0pGiiEYTgcDpMkMQ40BEEQRZFjMRQAAHCth8S9Ph8AQI24dxli1gMAALAiUAAAAFYECgAAwGrF\nOgpVZiul4NigEQAAO7Bcj4KssySFYSirLekbqyA7e1Mqu10AANTPEj0K8i91ecWVt9vtthCi0+lU\npNoSAADYrKKBgowDZJSg3xZCRFEkCzFtoXkAAKBMS1RmlP0HqdtCCFmCaf0BiGwpJzm6kQ1B5BoT\nhCYAAGxb0UAhCAK1xIN+W99hnXbIY+rRhud5/X5fCNHv9/WDh2HY6XTErKp01TIkAABwyRIFpGRe\nwnA4FEJ0Op12uy3/ppeX7XWyBeXghTy4jAnCMOz3++qYnuepNAjP89Ru2c4M90piAQBqxL3L0BKz\nHpIkabfbnU5HXtRHo5G83W6313xRZDqkviU1utFut2U0IGMF1cFgW8oSAABsxNLTI9VykVKSJGt2\n/gdBkJ00MRqN9OEG42CH3qp1GgAAAGxWLLi0ZkaCIrsEil/pjXmUKbaCSzkc6yYCAGBTNlDCOQzD\nFa7NQogoivr9vkx62KD8OkvUYgIAoLjNrPWg/5VfnBxukHMgZajR6XTy+ypSYxCyKyL/ISoUICYA\nAGBZG1jrYeV6BmEYqkGHIAj6/X673ZZXfZW9KKkRB7mbvn21pwYAAEVUaBaHPu9RTpiU3+q3RWaq\nZCqksM1LcW++CgCggty73Cwx9CDrIaqSiHK8QNp4s+RUiE6n43menDypxheGw2G/31fPS6cCAADb\ns1zBJdn/PxqN2u32aDSSeYhy5sKWAigZnRi3C1N2Aj0KAIASuXe5KXo+Mp9AXp5TYwFiftSgXAQK\nAIASuXe5WWJRKHVbL58MAAActsqiUKLaFYo8k7IbBQBALRUNFNQsg9RkSFVtqTodDJRUAgBgU5Zb\nFCq1dJMQQhY/4EoMAICTXMu5IJkRAFAi9y43mynhDAAAnESgAAAArBwPFCaTydXVVbPZFEI0m83L\ny8vJZFJ2owAAqA2XA4XJZPL8+fOHh4fpdCqEmE6ncRw/e/aMWAEAgIKWmx65qd124+7u7ujo6MWL\nF2rLYDB49erV7e1tia0CAKBGluhR8DxPruWYvSuKIlVQoXSqyNL19bUeJSiDwWD3rQIAoI6Wm8UR\nhmG/35e31QJR6tsqLOSo5qXEcdxoNGy7PTw85NwLAMBq3JseueL56EMMcu3pTTVoTZ7nCZE9o7mu\nDt/3ZdYCAACb5V6g8LXVHlapXATFPvSR6LHCycnJTpoDAEDtrRgo1NTZ2dnr16+73W7ZDQEAoB7c\nDBRUr4/ex+D7/pMnT+7v74+OjuZ3dqqPCACADXJqKMXzhBCefunXAoX0mbo3jAQAKJ17FxeXCy4B\nAIA1FQ0UgiDwPC8MwyrMgQQAALtRNFCIoihJkiAIZGElW+WlalAllzwh1BcAAFjackMPMj5IkiQM\nwwpHDMn87cRUWQEAACy2Yo6CHjFEUeR5XqXCBZlI4lg6CQAAu+dUcqY+60HmnTLrAQCwS+5dXJj1\nAAAArNwsuJRD62N419/gVvAHAMDG7FePgm0xiGqsjw0AQOWsGyhUKocRAABs1tKBgj7HIQiCTqcj\nCzFtvGVblSRCCI8RBwAA8i2dnOl5npjNPPQ8bzgcCiE6nU4Vkjxnsx5yPDay0WhOp9NGozGdPjze\nUX7zAQC1596sh1WSGWVwIHsRgiCQG6MoUrfLZpgeKd819e10OlX/AgAAm9VzFPr9fq/X22BTAABA\n1SzdQxIEwWg0arfbo9FIPlZuqUJPS07BpVSPgjZCkeg7AACwDveGHlZJZpRRguxOiKJI3V5NFEVh\nGGbXjAjm6ffKYY5lMyjjOM65l2EIAACyVhl6UKs8CCGCIFC3VxCGYafTkUFAp9PRMx5Go1HOQ8T8\n/IuUZrOp/lV8389pSaPRWL75AAA4bpVZD9mHyOkPKyQz6g+MokjNntBv5zxEPXB2l5BDD7aRBYYe\nAABbtddDD57nybmRXsZqz62KMchv9au+rY5TaqpFGIa2jgcAALC+JaZHykQE42QHmUaw7HPLYQt5\nW2YqiPlwQYUgvV4vZ3RjhZmZvu/Hcez7fm7eAgAA+26JQEFdqjdeh1HOmxCzCg1Kr9eTaYz9fl8+\nr0yl1B+obsdxLITMQtB7OB5vp3o9ZHwQx/FsB6e6iQAA2JSqDKVEUSQDAmOug5qBmUpKEPMpC8Vz\nFGzzJwEAWMde5yhIYRhmcxRWTlNQ1HRHY3aCnsegJyWkshwAAMBmLR0o9Pv9drudZKzw3DLmMN6V\nLaugtuvfsnYlAABbtZnpkas/veepREU5JVKOI+hLT6V2S91ut9vZ6ZG+35jlKj6WUWLoAQCwA+4N\nPawSKKxWMsEoDEOZqCjpsxv0zgY9GpDxhLpLb3/xEs4ECgCAbSBQeLxOy/kImwoXbKkGOSkIxrsI\nFAAA5SJQELasgiq8LusECvqkyiqcCwCgjtwLFJaooyA5dv5KkiTuvbsAAKxp6UChDh57PfTOD89b\nK8RJ9aMQTwAA9sQqq0eqUgphGMoVojfdqjU9TtdMkkSIx681L+1qCujKc0EBAKijpXsUZMmj4XAo\npyeooIHLJwAA7lm6R0FGCfqUB7Uw9EYbtl1qGcyyGwIAQKWtkqNQvbGGpWmlnMptCAAAlbZijoL+\nrfMLLnje45cQiXYbAAD3Ld2jMBwOO51OFEWj0ajdbsuUhV6vt43GVYEtJvA8ajQBANy3dI9CEASq\n314u5KgSG6uDFAQAADbCqdkKemVGbYsQpsqM2R2Mczeo+gwAKM69aYBL9yioCgpbaMymeIoQj1/0\nLgAAsIKlA4XhcNhut/v9vud5QRBUclZkohQvuNRsNuW/l5eXk8lkJ+0EAKDqVslRiKIoSZLhcCiE\n6HQ6le9gKGQ6ncp/4zg+PDzU+iQAANhf6w6lhGHY7/dFNZY/WCdHQVs90nwi5CgAABYiR+GRXOLB\n87x+v99ut2XvghOs764am1Bb6HUAADhvxbUehBAyPqh1naU4joXwi+x5cnL64sXj2IS+3bGwEQCA\nlFV6FIbDYZIkURTVOkoQQvi+NUrw/YYQXpKIy8srIbwXL17ssmEAAFTE0oGCMT6o5NyHtZycnMgb\ng8Gg3JYAAFCiooGC53kqGoiiSB+Yj6Ko0+lsvGXl6na7Qog4jlNjDQAA7JUVkxmrbZWCS6lhiMPD\nQ7mx0WhssaUAAFSbk4HCKgWXbD0HagwCAIA95GSgsEndbvf4+Pj09LTshgAAUAIChUeeJ4xFFFqt\n1v39ve/7cmwiZ6IEAADuIVBYrNVqXV9fqxrPZTcHAIDdKVppcmEyYBVKD61QwlnMVXFO7Zx+cVIl\nnB1bpBsAsD73SjgXrczY6/W22o4SaVGFlyTJZDK5urqT5ROazebJyUm32221WmU2EQCAkjgV+KzW\no6A93BuPx8+fPz86OlKlGE9PT9+8eXN/f99qtehRAADkc69HgRyFOXd3dy9fvtQLNg8Gg1evXt3e\n3hr3bzabl5eXk8lkVw0EAGCnnAp8Zj0KZgtXiPY8r9FoGNMVfd+fTqeyy+HwUA1DeGK+y2G95gMA\nao8ehVpYpeCSZJvUoGo5393dpe7K73IAAKDWnAp81s9RWNij0Gg0ptMH9YjUDuufAgCg1uhRqCvP\nM0yD1O59XAxiYSHnhV0OAAC4pPxAIYqiMAyDIMiuVS23h2GYfYhx+8rUUMV4PE4VbD47O3v69Klc\nTFIIYVsjiuWjAABOKjlQCMOw0+nIEKHT6QRBoO7yPK/f7wsh+v2+vl0+RMxWu86GFzpj909On1C2\nYPOTJ0/u7+/lYpLCvkZUars3L6eFAABUWclDKZ7nDYdDGQdEUdTpdGR7wjDs9/uqbZ7n9Xo92YWg\nP0Q9cLabSOUorNk2vQHaPeronhDi7Ozs9evXejBhPAIAYB+498lfZo+CvMCr3gL9qh9FUbvdVnu2\n2225XcYK6iFhGI5Gox00NUn0mRTvZLscAABwSZmBQhAE6ror0w7ELAgYjUb6cEMQBDkBQf7ow1ZN\np9ObmxuiBACAq4qu9bBVKg4YDof5e6Z6GvRgQrN0VoBj3UQAAGxKJQKFaKbT6aj8gzVw3QcAYDPK\nnx4pqemO+eMIqTGIVJYDAADYrDIDhTAMbYMEKntRUiMOqZigxOwEAAD2QcmBgvpXzHcPyOkMagbE\naDTS5zuoh/T7fT1lYecSWfCRQgkAAFeVPN1T1ktQ36piCam79O0ylUE9RG//9uoozG8375/d173Z\ntACAfO598lfifHJSDdS0ySIP2U2gICyxAoECAMC9T36nzmdTgUIqc8IeLnhJkuR0MLj34wIAyOfe\nJ39VZj1USjJvtYPIAKLZbF5eXk4mk022DwCAXXEyUPCMNv4E4rHvIe/I0+k0juNnz54RKwAA6sjJ\nQCEx2uYTCPV1eXmV2nkwGLx69er29naDDQAAYDecGkrZbDLjaprN5nT6oFqktvu+P51OS2kSAGBn\nyFFAnjiObdFAzl0AAFQWgcIm+b7faDSWvQsAgMoiUNiwk5OTpbYDAFBlTg2lVCFHYTKZHB62VIuE\nEGdnZ69fv76/vz88PCytWQCAnSBHAQu0Wi39W9/3nzx5QpQAAKgppwKfKvQoiLkCz+9e3oLVHgEA\nteZej8LXym7ANpirK5X7zslnd+8HCADgNicDBa7FAABsBjkKAADAikABAABYESgAAAArAgUAAGDl\nZDJjLTF/EgBQQQQKVcH8SQBABTH0sF2e55mrOgAAUAdO9iiUUHBp/imThc+YbSH9CACACnKyRyEx\n2sFTCktwMJlMrq6ums2mEKLZbBrjGPodAAAV5GSgUC2TyeT58+cPDw/T6VQIIf8FAKAWnEqdK31R\nKJmKaOkbMI5NyEcJ9S3JjABQa+59jNOjsDNO/dwAAPYEgcLmFQ8lGYYAAFScUz0kFRl6yGx8d3N2\nI5nf8m4kwr0+KwDYK+59jNOjUCbHfpgAAO4hUKgQNX/y8vJyMpmU3RwAANwMFDyjslslfN9X/xqp\n+ZNxHD979oxYAQBQOqeGUiqeo6BSEFJbMikLjy4uLq6vr7fYXADAppGjAINU4UXbwIHnCSFsVRYM\nBoPB5toIAMAqCBTWlS28uKmBgziOmT8JAChX+T0k0UwQBGEYqu1BEOi7hWGotkRRJL/V9xclDT1c\nXV0ZBwj0gQPLIlXyLuvQQ2Z/p/qyAMBJDD1sWBiGnU4niiIhRL/fVymHURSNRqOch8h9PM+Tjy2R\nbYBA354kj19CeNrtBc7Pz9VCU4792AEA6qLkwMfzvF6vpzoG1LdRFHU6HWPbPM8bDoeyd0H+q2KF\n3fcoxHHcaDRs9z48PKTuzUaaxh6Fs7Oz169f39/fHx4euhecAoDD3PvQLj9HQR9iaLfb8qpv6yeQ\nIYV6SBiGto6H3fB93xYo5NyVczT575MnT2SUsIEmAgCwhpIDhSRJ9EBhNBrpnQSqBEIqFyGl3NGH\nk5OTpbZL+kQJRaVD3tzcECUAAKrga2U34JEcaxCzPgOp1+sFQRBFUb/fl3dFUdRut9UOqYTHmaWr\nK63TTdTtdn/1q1998MEHKilBDhx0u13bQ+REiaOjo9SkBiosAQCqphKBQhAEo9FIjTuI+U4CFSvk\n9ytodjo81Gq17u/vb29vfd+P47jIwMHd3d3Lly+z229vb7fZUgAAllZ+zoX881/lJxqFYdjv95Mk\nUTfk9lTOYzUrM6q77I97fIjvN+I4zmY7lv4eAQAKcu9Du+QcBc/z2u12KlNBzHoRsvundit9bmRx\nyUxODaU4jnfZJAAAFioz8JH9ATIRQd8eBIH8+1vrKng3izJ1Wx+wqHKPgq7ZbM6HC4V6FFJ9Eo5F\nrADgBvd6FMo8HzmOkNqoLvz6dTGVviDTHiW9/XUJFC4vL29ubrQNjw85P7/44osv8oce3PsRBACX\nuPcpXenzkcGBMXfBeFddAoXJZPLs2TNtosTjQ96+HR8dHckjTCaTu7u7wWAwnU4bjcbJyUm32221\nWu79CAKAS9z7lC6/4FKOIAhsGY45d1WfnCjh+76ssKSoiRLbW2gKAIClOBX41KVHIfWQgotCnZ6e\n+r5/c3Pj0lsGAI5xr0fBqfOZBQpmOzjTNQMF9dBMtuMjWarBpbcMABxDoFBpde9RkA/NX2hKMN8B\nACrMvUCh0jkKNSLXpNBvrCx/oal1jgwAwLIIFDYjmbfm0YwLSp2dneUvNAUAwMYRKFRRt9s9Pj4+\nPT1VWxYuNAUAwDYQKFRRav6k7/tJknzzm9/81re+JYRoNpuXl5eeJ7JfAABsllM5F6UnM67AOOsh\ntcN4PJbLUr948WK22XqGNTp3AHCPe8mMlVhmej8Vz3m0LUsNAMC2ORX41LFHQYh3Qwa2HoVGo5Ep\nq5DuhMg/CABgN9zrUXAyR8EzKrtVq8tZmTrnLgAA1udkoJAYld2q1eXUX8ovzQQAwJqcDBRcs73y\nCc50ugAAtsSpoRRXcxTG4/H8stRiszkK7o2oAUBZ3PtEpUehBrJlFcpuEQBgXzgV+NS9RyF7j/6N\nPCsZq2b7D+hRAIAqcO8TlToKVebUjxoAoI4YeihfwdBzVqQ5IekQALAz9ChUQjZWIBoAAFSBk4GC\neaJfvQaNUo21nNFu2gIA2F9OBgr1CgkK0ZIWXUuTAQBUGTkK+2symVxdXTWbTTFbunoymZTdKABA\ntRAo7KnJZPL8+fOHhwe5WsR0Oo3j+NmzZ8QKAACdU/3YNa2jUFAq8yJ1ktqdnnGHlKurq+vr6+z2\ni4sL43YAQBHuDRA7dT5uBwr59EChyPk3m03jypO+77MiJQCszL1AgaGHGlt5Pac4jm3RQM5dAIA9\n5OSsh31hLOqcb2FI4fs+S1cDABR6FPZLkiQyvLi4uDDusL0lrQEAdeRkj4ILBZe2rdvt/upXv9KX\nrj47O3v9+nW32y23YQCASnGyRyExKrtV1ZJduvrJkyf39/eHh4dLHcebt53GAgBK41Ry5n7Oesjk\nKLx7T20zKlNJuevn6LqX5QsAq3Hv89DJoQc80rMdd/ak+fUeAAD1Uv7QQxRFYRgGQRCGYeou2/Yo\niozbIZVbj1kN9DDiAwAOKDlQCMOw0+lEUSSE6Pf7+h+jnuf1+325PQiC1EOEEFEUeZ4nH7vnUjEB\n9ZgBABtjTvzbFSFEr9fLftvr9fS26bsJIYbDobzdbrfb7ba2W1KBc9qR8Xh8eXkpax4cHBwIkcy+\n3rm4uEi0P+6V1Jb1X7KFTwEAe8K9T7/yhx703oJ2uy17CKIoarfb2e1yuEE9JAzD0Wi0q5ZWSGpJ\np6+++kq7U0UMiZr6qB7FcpEAgKWUHCgkSaIHCqPRSH6rbkhBEOQEBHs4+nB3d/fy5csXL17k7xbH\nU1WPeTfLRXre45cQiXYbAFBX5fcoSDLhQMz6DHJ203sa9GBC4y1ro6eyC6muAiGEWjQyRdVjzsYW\ng8Hggw8+uL293VSrbC9kDV9gAMCjSgQKQRB0Oh2ZbbCJ4y09ALOJJ90d+7pNnvaVZootxGAwuLm5\nEbPYqsizp2IsfSxjmZMAANRD+YGC53mj0Wg4HBYZQUiNQciHWPoVnLXUuk3y+t1oNHLWhHx4eCge\nMKk9x+Pxp59+qo9laPsIIby6BWAAAIOSAwXP82RHQupir7IXJTXikNptD7MTJNvSTQcHB0IIWZVZ\nmk4fhEjieJqdEyGtvFzk3d3d0dGRMU+CPgYAcEaZlSajKOp0Or1eL3X5D4JA3jUcDlO3hRCe5/V6\nPZnKIOMMFS7sTwnnyWTy7Nmz7JJO9/f3R0dHyYJVp9V9hpep4EvneZ6pl0I9OP0Uzr8jACC5V8K5\nzOmeslhCiqqLoN+r11oYDoe29u9bHYWLiwu1pNP5+fnbt2+T2RTeRqOhz5Oc/3p85WxfRdh+mjL9\nFoZjVucnEAA2zr2PtaoHPrJas3G7yIxE7E+Pgi67wtN0OjWNJjzu8/77/pdffmnsTnjcr8Crt36P\ngoNBNwC4+OHm1PkQKKhvm82m7SqeJHLdpvQlXI1WFAwULi4u4jien0xhDRQeHtKxi3u/SwAgXPxw\nK3/WA7bBlu24Qd1u982bN6enp0V2Xi1fEgBQOicDBUeqKq2j2+0eHx8XvIqvptVq3d/f+76v8iS2\n91wAgLI4GSiY0zHKbtVO5VzFZX3l1Q6bWi3i9va22+1m6ygAAJzhZKCwp1THibzRarWur683eBUv\nvloEfQwA4Iyvld0AbEx+r0mSGNZcWKqfRa4WoW+RmYzZ1SKm06mcfLFnAz4A4CB6FGos1YWwcP8k\nefyS9ZWXHY0xrhaRs92Ipa4BoF4IFGpshSSMZWMLxb4SlYjjuOBBdrPUNQBgg5ya7rmfdRQWKjip\nVwsb3t3SH2iqzSCEEL7vx3GcaHWjZamG1BYhxNXV1fX1dfYIFxcXxu0AUDvu1VFw6nwIFIyWDRSy\nl/nUjqnvz8/Pv/jiiyKBQk6owaQJAG5wL1Bg6AFm9nGJd78AZ2dnT58+7Xa7RQ6YP3hBoAAA1eRk\noEDBpUcrZyTkU7Mfnzx58jd/8zdHR0fyKRY+ylafceWlrgEA2+ZkoEDBpUcbegW82dcjlYp4c3Nj\nPL6KTlIHshWW3kHBaQDAapwaSiFHYR2pDIPUIlLGFIScI6S26E5Pz9SMyrOzs9evX9/f3x8eHm7q\nRACgROQoAOsaDO70wQuiBACoMqcCH3oU1rGNHgVh6VSwHQEA6s69DzdKOCNNlk3clIVBAwCgyhh6\nQFpqpiJlEwFgnxEoYIHsmk/bmHIJAKgmp4ZSyFFYR6aE8+OL6PuNbIXm/MvY8lcAACAASURBVCOk\ndtDiiXTew9qtBoBqce/DzckeBQouLS1nYSd51zqJC3tbxwIAHOBU4EOPwjpsPQpaqaV3/QH5RzDu\nYJxJsXpzAaCS3PtwY9YDsrb+I66XbnTsNwoAHEOggLWsNqRDcAAAdeFkjgJWYbl2GwIBNX+SxA8A\ncB6BAt5JkscvIbxGo2mMEoQQzWbD8/KiBP0uggkAqDUCBZgttaLjLLxIk1ECsQIA1BeBAsy63e7x\n8fHp6anacnbWXfloxWdXyr6K1BcAoCwECjBrtVr39/e+7+srPb59O1bDEzn5iNm7CpaFtsUExAoA\nUBanpnvO6iiYuXSm25aaB7zytOCrq6vr62ttw+NBLi4u57erJ7IeincPQC24V0fBqfOh4NKmbCpQ\naDab830J78pCp/oYZk802y+xbgGAKnMvUGDoAXNUuev1617HcTydToVIZl/ZuwAAVedU4EOPwpbY\nAuTUn/s5/QHyO3oUADiPHoVtCcMwtSWYF0WRuiuKoiAIsg9B1VxcXBq3LzX3EgBQokoEClEU9fv9\nVCgwGo2MO4dh2Ol05D6e5+mPQtV0u93j409PT89SSabd7uozLQEAu1RyD0kURWEYyphgOBwGQaC2\ndzodS3e3p/aU/6pYgaGHjUulKazwyk4mk9vb28FgEMfxCutPMvQAoF4Yeti8IAh6vV5qo62fQA43\nqHhCBRnYkmTeCkdotVrX19ekLgJATZUcKMhUg2y2gQwUvJn8dARGH0q0WvFEY/lFqioBQAVVepnp\nXq8n0xj7/b4QIgzDKIra7bbaQXUtzFt6Wp9j3UQ7k3qdGSYAAPdUNFDQOwlUrFB4mgPX/R1JkrlY\nYVMvO28fAFRH+TkKRaiegyAI9KQEGU9Y+hWwCwvXfSj4qCJLSAAAdq+igUKqcIK+Xf+W7IR6sc1l\nIDgAgMqqaKAwGo1ksQSp3+/LmREyUFBjEP1+X09ZQDXpZaHLbgsAYDlVme6pV0dQW9TtdrutOg9k\niQV11/zaRYI6ChVnrJSQ836RIAmgXtyro1Dp88lJQTDeRaDgHi1cfHeLNxhAZREoVBqBQk3lDkm8\nK+bo3q8fAPe490lV0emR2Cvylyr12yW/3VlWw/rFqgHASU4GCuY/UPnoRw5jsAIAcDJQ4KMeAIDN\ncDJQwH4pMjxB7AgAq6loHQUgpdlsyn8vLy8nk4naXjCJgQoOALAaAgXUg1yoejqdxnH87NkzPVYA\nAGyPU6lbTI+so8lkcnd3NxgMptNpo9H4/ve//3u/93s//vGP5bfT6cNsR1ufwLv5k0ZLlWwimRHA\nmtz7GKFHAWWaTCbPnz9/eHhQHQY/+9nP/vqv/1p9a3ug7/siXZfznRWacXV1ZRvdAIB9RjIjynR3\nd/fy5Ut9y9u3b4s8MI7j1JYkSVYL5GWwcnR0lBrduL+/b7Vayx4NABzjVA8JQw+102w2c7oNhBBq\nZCE19OD7fhzHekUmW+nGhUWgr66urq+vs098cXFh3A4AOdwbenDqfGaBgplLZ+qGOI4bjcaivcyB\nwvn5+RdffLFUoGDbwRas+L6/KIgBgDQChUqjR6F2VuhRODs7e/369f39/dHR0fqBQn6w8vDwUCCU\nAYB33AsUSGZEmU5OTgruKbMXfd9PkuSb3/zmt771LTErrrAO3/dtoUDOXQCwPwgUUKZut3t8fHx6\neqq2HB0dHR4eZveUHQ//+I//+O///u//+7//m50Woc9ZKDIDQu1g69IoHsQAgMMIFFCmVqt1f3/v\n+77qMPjud7/7ve99T32b2v/u7u7o6OjFixfZQ+lzFo6Pj8fjscgdhUqSRN47Ho9TwcrZ2dnTp0+7\n3e4GzhAAas6poRRyFGptyWWm1R3mbMeLi4ubmxv9CI1GUxZxOjk56Xa7cuqjfJbJZHJ7ezsYDOI4\n9n1f7mDs2ACAfO7lKDh1PgQKtbZkoCBmkUHR+ZNqh9PT0zdv3sgyCcYnXes0AOw39z5GnDofAoVa\ns12zF3UqmAMFyRgoSLJMAoECgM1y72PEqfMhUKi1nGv2fE5iKjJYukdB7TCdTgkUAGyWex8jTiYz\nekZltwqrS2ZkimIRv/3tb4UQOfMbVQUFFncAgBxOBgqJUdmtgpWK5FKBXTbCSy2+kJ0WIYT4wz/8\nw8PDw6+++kqYloRQ3nvvPXkvS1cDQA6nekgYetgH+lCCPrLg+404jg8ODj766KN/+Id/0B6Rl8Rg\nw88QgNW4N/Tg1PkQKOwD2yCSrNDcaDQyBZTSgcLBwYHsbzDtM3dAAFiWe4GCk0MP2F85K0e8//77\n8t8iUYLIm20BAHuEQAE1Y4zU1cb87EVRdMlKAMAjAgXUT5LIyMCTN/TQ4eTkRC/GbGRbxOH8/CJ1\ntK3yPMMXAFSNU0Mp5CjsidRUCL3cwng8fvbs2QcffDAYDNT9s93eFWye3+dxh7dvx7Jss74y9dZO\nwXoXP71ArZGjAJQvO/FVTaQ8PDx89eqVvspU9uGplaiU/MUdKMsBYD85FfjMehTMXDpTFJFaLUL1\nKKT20bscZhtFakvqmJtoW/opdtCNAWAH3OtR+FrZDdgGx94jbIYq4iT4EQGAwpwMFACD+TUd3m3O\neUhqhIHoAsAeIlDA3lk5wWDZB9qSLgGgRqqSzBiGoXFjEATZu6IoMm4HJH21iIIP0S/im7qgq1xL\nVhsBUF+VCBSiKOr3+1EU6Rs9z+v3+0KIfr8fBIHaHoZhp9ORj/I8L/UoQMxPi1Abs7UKVBmG7EVc\nvyt1L6tNAtgrJQcKsm9AXvh1srcgSZIoipIkGY1Gqv+g3+8Ph8MoiqIoarfb9Ctg21IxAatNAtgr\n5fcoBEHQ6/VSG2UQoL5tt9uy50CGBaqDIQzD0Wi0k2airmxd/sWHAu7u7lJbBoPBq1evbm9v9Y0U\nWgDgpJIDBZlqkO0VGI1G+nBDEAQ5AQGjD8iXGkdYtk6zVuQxbzsZCQCcVH6PwlJSPQ16MKHxlrWr\n5qNO1I/H/IqUifqK4zjnp2gymVxdXTWbTUFaA4A6c3J6JH/OYQNUcehGo2FZvDrRqz3OQgVVC7Il\nxA+F+KEQYjoVNzfi5kaMx5NWq6UfQosw+KkFUEU161FIjUHIQQdLvwKwGScnJ6enZ/n7FOyWOjxs\npbao0YqVmgYAW1fRQEFlL0pqxCEVE5CdgB3odrtv3rw5PT0Twpt9AcC+qGigIKczyDggiiI1PVIG\nCvpUST1lAdiG1GqT2RUp55MPPMvXo6llGAMAqqmigYKcM9npdDzP63Q6vV5P9SUMh8N+v6/Sx+hU\nwDakUhFvb2+73a68xutXernCVHZAIUez2SR/FkCNVH01TFmRybhdZEYiZstMV/ycUHWTyeT58+dH\nR0cvXryQW05PT9+8eXN/f394eKgvXW1iu+/dYtb62taNRnM6nTYajen0Qe0AoL7cW2baqfMhUMBG\nXF1dXV9fZ7dfXFzc3NwkSXJ1dXV9/UPTQ3MiCHOgoD3k3Q4A6otAodIIFLARzWbTmEng+34cx0mS\n2HZIOTg4+PrXvx7Hse/7cfy4P4EC4Db3AoWK5iish6pKWF0cx7YgII5jIUSmBJPVX/zFX2TTGgCg\nXpwMFBKjsluFevB9v9Fo2O6SP0u2HZSzs7OnT592u90VGkBUAaBSnAwUgLWcnJzkb7ftcHBwIITw\nff/Jkycy8zG7j5xJkUOPQugVA1A6p4ZSyFHARkwmk2fPnn3wwQdq2aezs7PXr1+ra3/ODkdHR6mf\nPzlgab/KL85RyBnyTEUP/OgDpSNHAXBftsJSqodg4Q6S6gZYqjPA8+a+8rFkJYBtcyrwoUcBG7fw\nj4PUDrb9m82mqpSgP1r+p0+LMD7JUm0AUCL3fh+dOh8CBWzcRgKFOI7z8x/n159c0CLTA536RQZq\nzb3fR4YegM3QBxpSYw35MynU7SR5/BLCs3/OMNYAYKcIFIDNyJ+Oa5wocXZ2ZptAMZlMPv30u/nP\nmFqQ4vLycn55KgDYACcDBQouYQNyeggK7qDrdrvHx8enp6dqi5woodda0C/8H3744X/+538a159U\nOz9//vzh4UHVdIrj+NmzZ8QKADbLqaEUchRQZZPJ5Pb2djAYyKLOJycnX3zxhb7Dp59+qq9ElTG3\nWsTl5eXDw4N959lj+GUAdsu9HAWnzodAAbVg/ByxrUSleXyI7zfiOLZ9GPm+P51O3fuoAurCvd8+\np86HQAG1YPwcKbDQVHYdKbOHh4dms8mvAVAK9wIFJ3MUgJrJWYlqWTkzLABgBQQKwO7Y8h83eHX/\n6KOPmAcBYIMIFIDdyZlCmb/Q1Pvvv68fZv7r0eeff/6Nb3zjww8/ZB4EgA1yaiiFHAXU18KFpnIn\nYHq+73/00Ue/+MUvsvddXFwsSpMEsDHkKADYioULTV1cXOY8fDqd/vM//7PxLhV5bA81SwCHORX4\nzHoUzFw6U7jNuH5EtstBnwcxnU5zshweHh52kOHo3h9SwArc+0VwskchMSq7VcDS9OTHw8PDVJeD\nvmf+chL6Xfz1D2ApTgU+5CjADQWXrNSu8nnX+/Pz85ubm+JPkYoeiv86ufeHFLAC934RnOxRAOpq\nqfUjFNlnNh6Ps8tJPH36VF9OovjRROEFKlmbCnAbgQJQIQvHy/RIInXXwnTIleWMVrA2FeA8p3pI\nGHrAXjH2OMgfflvn52Qyubu7GwwGMvnx5OSk2+22Wi3TwdNHWGqJCuZkYm8x9ACg0nLGK4r89V9k\nHEHvXbBFAzuYkwlgNwgUgLpa9o+Wu7u7ly9f6itTDwaDV69e3d7eym8LjiOoYZGc9Sk2uHoFgHI5\n1UPC0AP2mXEShP7LYFmg0vjb4uXvkCSP/au2RS/latcFWg24hqGHWvCMym4VsCPGdEjLn/i2j7Mk\nfwf1+2RbosK2HUDtOBX40KOAfaYu3raff9Nf/zm/Kt6iHd45PT3LLlGhz7bIr81QMIzn9xq1QI8C\ngKrzvNTXo5yxgCQRQniNRjPnqLMvg8HgLn9OZk5thuKdffq5FH0MgLU5FfjQo4B9lnP1lCkF4/HY\ntlpEo9GcTqcHBwdfffVbdTx9h+JZC7ktNOyw5EXfqY8sOIkehd2JoiiYp98bhmEQBGEYltM4oHpy\nPpo8TwiRHB62Xr16aVwtQnY2fPXVV8Weyst2MMinUN0YuoKlG2WvRpKI7BeAElU6UBiNRsa7PM/r\n9/tCiH6/nwoggH1W5Pp6c3OtZj/mHy217tRS9LGPw8NWdsrlyoMIVIkGds280mIFtNvtdrud3d7r\n9fRmCyF6vd7sdlLtcwJ2ajweC5Fkv8bjsXgcolNf7z4S1BbLDoZPkWJfcy4uLpLHWOZxB7nEZaPR\nuLi4GI/HlrMQQojT09Pj42N9H6A63LsIVbdHYTQaGXsLoihqt9vq23a7HUXRzloF1Mjd3Z1xmEBm\nGi7117zsXTg4OPjkk08yd3rn5xfLjhRkSzfaqjzd3d1lH6vXiQKwVdUNFIQQURTJzskgCFQ0kAog\ngiCwjVAAe85eR7lIP8FcJCGv4v/2b//2B3/wBwsXqNQiBmssEsex7EIwNluPA2xnQZVoYDcqHSiM\nRqPhcDgcDoUQnU6n8OPMBZdybO8UgLKsXUTZEEkUXKBSz160ee+99+I4ztlBxgE5Z0GVaGA3vlZ2\nA6wSrQdTdi2EYVhsmkNCmjTg+36j0Zi/lHrGXoTT07OFf5x72trWSZJcX1/bCjPIBSOOjo7yr+K/\n+93v8p9RxgGNRqPRaBiPJE9wQbsBrK3SPQo6chGAZZnqKHtnZ109pWA8nug9BOfnF9njJMlccpbq\nhzN2yGWXnko5ODgo0ngVB9iqQcdxTL8gsAMVDRRkEQV9i0pNSEUMqdxGAEq32z0+Pk6lFLx+/VpP\nKWi1WtfX7yZM3tzcLCxjkMqITt2bnzrw/vvvF6zWoOKDVAKEmCVGvH37NrFUezTKFKw07sOgJJBW\n0UBBpiiqgQZ5QwYKYRiORiMZK8haC5RdAowKphRsUG7eQCJE8uWX8cI8ylSCZKvV0u+VZ9Hv9weD\nwcIiToolLMg0cRZzFAw+gL2wxamX65H1EhRVLCF1l76dOgqATc4vxmY/EzJ5AwsrKzx+K4ckfN8/\nPz+XvQVaC+eqO4zH408//VQfklhYWcHWBt14PL68vLSVcwAKcu8iVPWS1LLnwFZQIbWdtR4Am53V\nn7+8vLy5uZnfZnveueUk3r4dHx0dqUbO9/y/W3IiSZKrq6vr6x8WaIvhCPLw2ZU2VQ6myq44PT1N\nDaPwyYIi3FvrwanzIVAAsvKXeN64yWSSWnoqtfC0ao7vN+I4Nq47pYcLSZKsnC2glqpKRQbZQOHq\n6ur6+jp7hIuLi5ubm5wXbccvL6rPvUChojkKADYl1Yu47acrnhiRymbYRiNVEoPakspSVHmL19fX\nxrEJ2TuSu5bVY4N38/IWRFYmNsipwGfWo2Dm0pkCtWD800q/TAshUoMC0mQyubu7GwwGspTCycnJ\nzU32z/30L7vv+9Pp1PM8Uy/Fu2fRd5DfTqfTRiNnBSxPCHF6evrmzZv7+/tUcmXOmZaumq1ynnsv\nu5M9CuZ0jLJbBeyR/FoLku0XU6YLpBacPD7+dDyeyOmaFxeXxj8JPvroo/xykKZ2CiGS3CjhkXGB\nCW2m5eMS2yxrCfc4GSgAKFmRMN1WzCBbsil1kc7Wh/j888+/8Y1vfPjhhwuLOheIJNJraOn09EZj\n4w8PW6lYYfejAHoJbVblxvoIFABUy8JVoLJpEP/yL//yX//1Xz/+8Y8zD0ovfzWLJMwX7B/84PP8\nthVZYCLV67CDJIZULJLtj9FX4wSW5dRQCrMegFowVTp6vJGzqqQQ4uHhIXWvHA9uNpvz1++cj4C8\nv+n/5E/+5De/+c1stSrDdAyZBjF76tQxDfuPx+NUskW32zVmORRnm2chXwrj9I3T01Pf943TOrBx\n5CgAwLpyqkTnLPWUuktPg8j8lW+LBhb0/P/mN7+Rh7q4MKx5IeYXmMg/VJIk4/F4I3/cpzoMUl0U\nqYGGv/zLv8weYTAYsCo3VkagAKBabKtApbbrORCm2MIzfS0g4wAhxHzNqPzKkmbNZvPP/uzPsskW\nH3zwgT42kV2BwvPmIgMxP3gxHo/1sKDb7f75n/+5HovYltJgVW6szKkeEoYeAAcsLNmUZSoHKb79\n7W//0z/901dffeX7/m9/+9sii1FZRhYWSg09GO81yhsikR9k+kzR999//xvf+Maf/umf/uhHP5I7\nffvb3/7lL39ZpIn6qWGrGHoAgO1aYS0r4zqZv/zlL2VwEMfxsktWivm6DjkW5j/ORwPL9Ul0u93D\nw9b19Q+n0we5pNb//b//9qMf/R+1T8Eo4ezszNZPAyzkZKDgGZXdKgBFtVqtm5sbmVQYx/EXX3yR\nv+KlMbbQV5Yaj8epSOLo6Eg/ZmrJSknPori8vDIOZ/zrv/6rfNICQxv5KZbpwZHpdHp3d2vZf4k/\nWLNriwNLcaqHhKEHYJ/lL7swmUxub28Hg0Ecx77vf//73/c876/+6q/kt3I+Qk44kplYUURqHoT5\ng2k8nrRaLc/zLi8viy12pR/T7ODg4Otf/3rBU8NmuTf04NT5ECgAyJHfs5jzwVFk0qYpkiiUtfCd\n73xHzrbwPC9J/r9FDzFMwlzImc/E/FiwItwLFJwcegAAA3N19wJV3otM2lw5CeDv/u7vZISRJInl\n8r8gJvjkk0/U7bOzs2UnaBSRHcktpeLk5eWlfLUbjcbFxQVVpHaDQAEAFls4adOUUGlMC0hfU3/3\nu99ldigyq/NdHPCd73xHJWcY0xqyMzCXkr1C/+3f/u2Or9nGFUCoOLkbBAoAsJhxYoWe/2hJqByr\ndMhGo7nUeIGidxgY3d3dqsTPgsdMVWvIr+VweHioX6H/4z/+47vf/e5kMtnlNXvhCiDYHqeGUshR\nALA9qXRIY5JgziC6sdhDvrOzs1//+tcff/yxnnT5xReGg6jnUc+f3ZLSaDRlYenp9CFnh4ODg48/\n/vjv//7vl2q52PQHsS2ZtIL1IdzLUXDqfAgUAGzVOsl02UJSBwcHOQUeVCxydHRke9KcsMDWtMwS\nFWKlUlEGBwcHv//7v7+9VS1ssiuAlItAodIIFABU2WQyWXaaYv7HWcFAYf6Km503YZvGuU6KoqHZ\nS30y6yUpG43G//zP/xiDKnoUdsDJHAUKLgGoolarJWdYDIdDYxkomfegl4pa/0lTqYjrH7AAc7NT\nH8M5H9HZ1MWPP/7YeEwqTu6AU4EPPQoA6iXbx/D27dvivQ4LexTkFffo6EhLA9xSj0LBj913x7y8\nvLQtwH11dfXw8KCnLkqffPKJypYwrgBShUIL9CgAADYm1ceQJMlqJRRt8x6zkwU0K1ZZ+PzzzwuX\navB+9rOfz1803z3k+vqHtnkTg8HA2OZf//rXOSuAFC+0MBqNbGe38lTSnGPWnVOBDz0KAPZKzmVM\nfgwuUy9SFO5RWLrzYI2iT0WLaWb7Tk5PT9+8eXN/f686KrL9N+PxWM+4tL2eOdcUY96JYxchehQA\noK5s1yO5PY5jU6Kf7dJrvSSrzAnZ7WFvzlydqD/+4z9Wd3z72/+P/VH55p7O931bMc27u7uXL3/+\n4sVA9VgMBnevXr1UhRZkJKHnNJyenqpujNFolF8HwthhYDym3L70iVZZfk3Tepn9SJXdDgCohoLZ\niz/4wQ/0y0J2NGE+v/LxS8lmZRa8/hT7euf8/Fw/u4JHkC4vL5dpkvWlGI/HqgG2Y15cXGz1bd0x\npy6qBAoAoLu4uDBeyeQVV02+WOn6PfdEmYOs4r333stvhjFeWfj18PCQ2GMmy5MavxWnp6fHx8cy\nVoiiKGcFkB2+yVu37zkK7qWnGnGajtmTM+U01yerPL169Urf+PTp01QaoBRFURAExdP3jK3+2c9+\n9r3vfW+VtmYOn/+klnYaEiOSRPz85z//7LPPlnnSvJXBv/OdP5ULfuYcqGploNbh1O8hgYINp+mY\nPTlTTnMjilSezjRp8WFzmmwrt7y8FbImDQ+Xtajtq3jnVKNamPhp3MerYBmodTj1ezgLFMyMZ8on\nkUv25DTF3pwpp7lZssNgB09UcFWLP/qjP/rv//7vL7/88v333//yyy9tu8k/zXNjlzVnZwjTVb9g\nfWvjva5NvnNy1oN5lKXsVgFAmXYTJQjLSpvy2VUVhPPz8yiK4jgeDodxHOcM9su77B/h5ihBPtEm\nFFny23FOBewMPdhwmo7ZkzPlNOvLWHHy6OgoSRJjx4atE+L8/Ny4PSe2EEL89Kc//eyzz5rNpmlh\nzKLLXvi+/9FHH/3iF7+Yf3h+B4PcgR4FAABy5VScNHZsyE6I1ManT592u13j8VVPg/Guzz77bFZD\nwjN2Cfz0pz+dv5KnZzcMh8PpdPqTn/wk2yrbky4qMlFjdY1kwzAUQgRBoP/MbbtHYdnAf3sHr05L\nODgHr2NLOHgFD55KuozjOH/ZC2MnxNnZ2ZMnT+R2W06lPHiSJAvrWopZ14h8yKyP4f/NPGlXPamT\nBYLr16MQRZHneVEURVHU6XRkxAAAqLVWq3V9fT2dTuVf80KI/KkZxk6I169fq04I48KSZ2dnant+\nXUvVKiGEatVPfvKT4+NPU70U+pO6Ka/IQiW12+12uy1v93o9/RRWKLi01N7LvlzbO3h1WsLBOXgd\nW8LB3Tj4eDy+uLhQeYvn5+d6RaacVbzXaYn+pDIrM1MGqo6X1jz1G3rwPG84HKoRB/1bhh44OAev\n8sGr0xIO7tjBZQdzdn9bDYmNtMSYlenk0EPNAoXsT4Pneb1eTw5AEChwcA5e5YNXpyUcfN8Onrqo\nb/OHXLgXKHyt7AZsQBRF8xu84iVIl917uUNv8+DVaQkH5+B1bAkH5+Dba4ljXAgUVJwox5LKbAoA\nAG6p36wHAACwMzULFGTnQWqsYWd1SQEA2Dc1CxSEEO12u9PpyNsyRCBQAABgS+oXKMjuBM/zPM8b\njUbD4bD4A4N5+r1hGAZBUOvyTcbG285Lvhp1PN9sm1Nvq97hVNPTjKLI9sa59IbaTtOxN1SdZibt\n2ql3U9jP1LE3VMn+jerYG/rOStUXyjccDlUJ8YJkdaa2Rt0lX4p2uy3/3WhLd0QGTKnXxHZe+kuR\nfVSVZU9TbtHfVnVvTU8z1Wz9l9SlN9R2mo69oak2F/nYqeNpJvYzdewNVbJtduwN1dU1UFhBKjhQ\nMuUdRa/X21mr1jccDtVHrf4jmHNe+p62l6VqbKcpP4aMD6njaSaZn0D1rWNvqO00HXtD9Tbrp+bY\nu5nYz9SxN1RSPdmpoEft4MAbqtujQMEWAaTeudq9kcPhsNfryR9T/QpqO6/UD3TOr3Gl2E4zdTq2\n7XU5zWT+YyXR3jjH3lDbabr0hmYbqc7asXcz50xdekMVeTXJiQDq/oam1C9HYR1yQSnP8/RxstFo\npA81BUEwGo1Kad5q5NBXdvRrqfPKjp5Wje009ZwVz/PyRwGrf5pCiCRJ9DdOvY+OvaG203TpDQ2C\nIJn1SKvKgE6+mzln6tIbKgVBoMoBK469oSkuFFwqTiU/hmFoLAy+D6IoUn34wpU5I71eTwZ//X5f\nCBGGoQOnKQuWC0uaqr5brc/UeJqOvaHqsrEw+brWpynsZ+rMGxqG4Wg0Kn6lr+lppuxRoKCHBbJr\nwfgXKupF/41Vn0QOvK3yA7fdbtfuj4+lZE/TyTc0mul0Ovqadu7JnqlLb6hsfPGpds7Yr6EHnfMf\nwTapPjH5IjjzyaVOpNanqab+FvkRre+ZFjlNN95QMRs7E4u6net+mmLRmdb6DZXnJedAytUfOp1O\nfpvreJpZ+xIoRJn1QNWQUipiSPUU1ZftvFKvQ92jpcA0PV3U+TQ9z5NpUKlTcOwNtZ2mS29oGIa2\nxYQcezdzztSxN1SOocg0BSFEu93Ov47U8TQNysqi3D2hzXrQE1b1ElGQ1QAAAlVJREFUqfnGagR1\nIUzzBo3nJean7tRrlocwzV3Wv9VPrXanKd+pXq83nJe49YbmnKZjb6jeZts7WPd3U7KdqWNvqE5k\nZoS69Ibq9ihQkMGBok+V1O+qVxEFXTbEsZ1Xaoxtx+1cU/Y09XPRfwnreJqpn9LUSTnzhuafpsNv\naJGPnTqeZpJ7pi69oTphmqrtzBuq85I9y/zPGSLKDk+4wXZeNR0tM8p/W2131RRvqO2uylrhY6eO\npynszXbsDc3h2Bsq7V2gAAAAituXZEYAALACAgUAAGBFoAAAAKwIFAAAgBWBAgBgT3meF80WC5S3\nZfEouXag2k1tlLX/9epS+mPlv7s/i20jUAAA7C+5LIWcANjpdKIoSpKk1+upxZ/kEg+yNEKv15Pr\nWkmywKgsFybXNnMSgQIAYH/JqszyhpgVPNCXqwjDsN1uq4UeVI1/tY8s6uzwYlEECgCA/bWwCJJa\nGCi1f+TEEtJFECgAAJDH4SCgCAIFAADypFaGlDeMS0g7iUABAAArPYExiiIVHKisBf1bJxEoAABg\nJRMY5RzITqeTWiWy3+9np1M6hkWhAABYQC3/GIZhFEWpwYggCKIo6nQ6Tl5S6VEAAMAqCALZYSD7\nDPr9vryhKizJbx0eeqBHAQCAPHopxna7rboTwjDU6y8Nh0MnByAIFAAAWEx1HhTc7gwCBQAAYEWO\nAgAAsCJQAAAAVgQKAADAikABAABYESgAAAArAgUAAGBFoAAAAKwIFAAAgBWBAgAAsCJQAAAAVgQK\nAADAikABAABYESgAAACr/x+qu8TWCRkODgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w.pdf(\"model\")->fitTo(*hist_data) ;\n",
    "TCanvas* c1 = new TCanvas();\n",
    "RooPlot* frame = w.var(\"mgg\")->frame() ;\n",
    "hist_data->plotOn(frame) ;\n",
    "w.pdf(\"model\")->plotOn(frame) ;\n",
    "w.pdf(\"model\")->plotOn(frame,RooFit::Components(\"bkg\"),RooFit::LineStyle(kDashed)) ;\n",
    "frame->Draw() ;\n",
    "c1->Draw();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now save the workspace with the data a modelconfig so that you can use RooStats to extract limits\n",
    "\n",
    "Create an empty ModelConfig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RooStats::ModelConfig mc(\"ModelConfig\",&w);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the pdf, the parameter of interest and the observables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mc.SetPdf(*w.pdf(\"model\"));\n",
    "mc.SetParametersOfInterest(*w.var(\"mu\"));\n",
    "//mc.SetNuisanceParameters(RooArgSet(*w.var(\"mean\"),*w.var(\"width\"),*w.var(\"alpha\")));\n",
    "mc.SetNuisanceParameters(*w.var(\"Bscale\"));\n",
    "mc.SetObservables(*w.var(\"mgg\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the current value $\\mu=1$ as an hypothesis "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== Using the following for ModelConfig ===\n",
      "Observables:             RooArgSet:: = (mgg)\n",
      "Parameters of Interest:  RooArgSet:: = (mu)\n",
      "Nuisance Parameters:     RooArgSet:: = (Bscale)\n",
      "PDF:                     RooRealSumPdf::model[ S * sig + B * bkg ] = 2.77715\n",
      "Snapshot:                \n",
      "  1) 0x7fe4dca0cd80 RooRealVar:: mu = 1 +/- 0.486745  L(-1 - 6)  \"mu\"\n",
      "\n"
     ]
    }
   ],
   "source": [
    "w.var(\"mu\")->setVal(1) ;\n",
    "mc.SetSnapshot(*w.var(\"mu\"));\n",
    "\n",
    "mc.Print();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import model in the workspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.import(mc);\n",
    "\n",
    "w.writeToFile(\"model.root\") ; "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ROOT C++",
   "language": "c++",
   "name": "root"
  },
  "language_info": {
   "codemirror_mode": "text/x-c++src",
   "file_extension": ".C",
   "mimetype": " text/x-c++src",
   "name": "c++"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}