Stats/hepstat_tutorial_hypo_p0.py.nbconvert.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ready-separation",
   "metadata": {},
   "source": [
    "# Hepstat_Tutorial_Hypo_P 0\n",
    "RooStats example: compute the p0 and significance (Hypothesis Test) \n",
    "The signal is a simple Gaussian and the background is a smoothly falling spectrum. To estimate the significance,\n",
    "we need to perform an hypothesis test. We want to disprove the null model, i.e the background only model against the alternate model,\n",
    "the background plus the signal. In RooStats, we do this by defining two two ModelConfig objects, one for the null model\n",
    "(the background only model in this case) and one for the alternate model (the signal plus background).\n",
    "\n",
    " Objectives of this tutorial are the following:\n",
    " * Compute the null hypo significance using the Asymptotic calculator\n",
    " * Compute the significance by hand using the asymptotic formula\n",
    " * Compute the significance using frequentist method\n",
    " * Plot the p0 scan as a function of the signal mass\n",
    "\n",
    "\n",
    "Based on the example [here](https://www.nikhef.nl/~vcroft/RooStats.html), and also [StandardFrequentistDiscovery.C](https://root.cern/doc/master/StandardFrequentistDiscovery_8C.html)\n",
    "\n",
    "\n",
    "**Author:** Lailin XU  \n",
    "<i><small>This notebook tutorial was automatically generated with <a href= \"https://github.com/root-project/root/blob/master/documentation/doxygen/converttonotebook.py\">ROOTBOOK-izer</a> from the macro found in the ROOT repository  on Wednesday, April 07, 2021 at 09:25 AM.</small></i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "available-bracket",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:45.503064Z",
     "iopub.status.busy": "2021-04-07T01:25:45.502300Z",
     "iopub.status.idle": "2021-04-07T01:25:45.504115Z",
     "shell.execute_reply": "2021-04-07T01:25:45.504549Z"
    }
   },
   "outputs": [],
   "source": [
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rubber-basis",
   "metadata": {},
   "source": [
    "Import the ROOT libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "iraqi-forth",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:45.508190Z",
     "iopub.status.busy": "2021-04-07T01:25:45.507622Z",
     "iopub.status.idle": "2021-04-07T01:25:46.000288Z",
     "shell.execute_reply": "2021-04-07T01:25:46.000498Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.22/07\n"
     ]
    }
   ],
   "source": [
    "import ROOT as R\n",
    "from math import pow, sqrt, fabs\n",
    "R.gROOT.SetStyle(\"ATLAS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cheap-conversion",
   "metadata": {},
   "source": [
    "Prepare the model\n",
    "=======================\n",
    "Signal mass point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "grateful-intent",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:46.002281Z",
     "iopub.status.busy": "2021-04-07T01:25:46.001981Z",
     "iopub.status.idle": "2021-04-07T01:25:46.109207Z",
     "shell.execute_reply": "2021-04-07T01:25:46.109592Z"
    }
   },
   "outputs": [],
   "source": [
    "mass = 125"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "genuine-guarantee",
   "metadata": {},
   "source": [
    "Worksapce"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "touched-quest",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:46.113520Z",
     "iopub.status.busy": "2021-04-07T01:25:46.112993Z",
     "iopub.status.idle": "2021-04-07T01:25:46.223670Z",
     "shell.execute_reply": "2021-04-07T01:25:46.224356Z"
    }
   },
   "outputs": [],
   "source": [
    "wsfile = \"test_hf_ws_{}.root\".format(mass)\n",
    "if not os.path.isfile(wsfile):\n",
    "  print(\"Error! No ws files found: {}\".format(wsfile)) \n",
    "  pyhist = \"hepstat_tutorial_histfactory_hists.py\"\n",
    "  if os.path.isfile(pyhist):\n",
    "    cmd = \"python3 {0} {1]\".format(pyhist, mass)\n",
    "    os.system(cmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "third-stanford",
   "metadata": {},
   "source": [
    "Open the workpace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "australian-operation",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:46.618092Z",
     "iopub.status.busy": "2021-04-07T01:25:46.617777Z",
     "iopub.status.idle": "2021-04-07T01:25:46.725182Z",
     "shell.execute_reply": "2021-04-07T01:25:46.725584Z"
    }
   },
   "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": [
    "tf = R.TFile.Open(wsfile, \"READ\")\n",
    "w = tf.Get(\"myws\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sustained-engine",
   "metadata": {},
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "neither-spelling",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:46.739220Z",
     "iopub.status.busy": "2021-04-07T01:25:46.738718Z",
     "iopub.status.idle": "2021-04-07T01:25:46.848218Z",
     "shell.execute_reply": "2021-04-07T01:25:46.848820Z"
    }
   },
   "outputs": [],
   "source": [
    "data = w.data(\"obsData\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "periodic-newark",
   "metadata": {},
   "source": [
    "The S+B model (Alternative hypo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "consecutive-seven",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:46.933332Z",
     "iopub.status.busy": "2021-04-07T01:25:46.932987Z",
     "iopub.status.idle": "2021-04-07T01:25:47.040058Z",
     "shell.execute_reply": "2021-04-07T01:25:47.040426Z"
    }
   },
   "outputs": [],
   "source": [
    "sbModel = w.obj(\"ModelConfig\")\n",
    "poi = sbModel.GetParametersOfInterest().first()\n",
    "poi.setVal(1.)\n",
    "sbModel.SetSnapshot(R.RooArgSet(poi))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "promising-competition",
   "metadata": {},
   "source": [
    "PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "lonely-craps",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:47.051575Z",
     "iopub.status.busy": "2021-04-07T01:25:47.051085Z",
     "iopub.status.idle": "2021-04-07T01:25:47.161729Z",
     "shell.execute_reply": "2021-04-07T01:25:47.162355Z"
    }
   },
   "outputs": [],
   "source": [
    "pdf = sbModel.GetPdf()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sensitive-northern",
   "metadata": {},
   "source": [
    "save snapshot before any fit has been done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "saving-bridal",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:47.182173Z",
     "iopub.status.busy": "2021-04-07T01:25:47.181558Z",
     "iopub.status.idle": "2021-04-07T01:25:47.292395Z",
     "shell.execute_reply": "2021-04-07T01:25:47.292975Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = pdf.getParameters(data)\n",
    "snapshotName_init = \"snapshot_paramsVals_initial\"\n",
    "w.saveSnapshot(snapshotName_init, params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eastern-completion",
   "metadata": {},
   "source": [
    "The B model (Null hypo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "floral-effort",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:47.316638Z",
     "iopub.status.busy": "2021-04-07T01:25:47.315935Z",
     "iopub.status.idle": "2021-04-07T01:25:47.429076Z",
     "shell.execute_reply": "2021-04-07T01:25:47.429683Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "RooWorkspace(myws) SR workspace contents\n",
      "\n",
      "variables\n",
      "---------\n",
      "(Lumi,alpha_background_shape,alpha_bkg_norm_uncertainty,alpha_signal_norm_uncertainty,binWidth_obs_x_SR_0,binWidth_obs_x_SR_1,mu,nom_alpha_background_shape,nom_alpha_bkg_norm_uncertainty,nom_alpha_signal_norm_uncertainty,nominalLumi,obs_x_SR,weightVar)\n",
      "\n",
      "p.d.f.s\n",
      "-------\n",
      "RooRealSumPdf::SR_model[ binWidth_obs_x_SR_0 * L_x_signal_SR_overallSyst_x_Exp + binWidth_obs_x_SR_1 * L_x_background_SR_overallSyst_x_HistSyst ] = 1.353\n",
      "RooGaussian::alpha_background_shapeConstraint[ x=alpha_background_shape mean=nom_alpha_background_shape sigma=1 ] = 1\n",
      "RooGaussian::alpha_bkg_norm_uncertaintyConstraint[ x=alpha_bkg_norm_uncertainty mean=nom_alpha_bkg_norm_uncertainty sigma=1 ] = 1\n",
      "RooGaussian::alpha_signal_norm_uncertaintyConstraint[ x=alpha_signal_norm_uncertainty mean=nom_alpha_signal_norm_uncertainty sigma=1 ] = 1\n",
      "RooGaussian::lumiConstraint[ x=Lumi mean=nominalLumi sigma=0.02 ] = 1\n",
      "RooProdPdf::model_SR[ lumiConstraint * alpha_signal_norm_uncertaintyConstraint * alpha_bkg_norm_uncertaintyConstraint * alpha_background_shapeConstraint * SR_model(obs_x_SR) ] = 1.353\n",
      "\n",
      "functions\n",
      "--------\n",
      "RooProduct::L_x_background_SR_overallSyst_x_HistSyst[ Lumi * background_SR_overallSyst_x_HistSyst ] = 2.706\n",
      "RooProduct::L_x_signal_SR_overallSyst_x_Exp[ Lumi * signal_SR_overallSyst_x_Exp ] = 0\n",
      "PiecewiseInterpolation::background_SR_Hist_alpha[ ] = 2.706\n",
      "RooHistFunc::background_SR_Hist_alpha_0high[ depList=(obs_x_SR) depList=(obs_x_SR) ] = 2.68\n",
      "RooHistFunc::background_SR_Hist_alpha_0low[ depList=(obs_x_SR) depList=(obs_x_SR) ] = 2.937\n",
      "RooHistFunc::background_SR_Hist_alphanominal[ depList=(obs_x_SR) depList=(obs_x_SR) ] = 2.706\n",
      "RooStats::HistFactory::FlexibleInterpVar::background_SR_epsilon[ paramList=(alpha_bkg_norm_uncertainty) ] = 1\n",
      "RooProduct::background_SR_overallSyst_x_HistSyst[ background_SR_Hist_alpha * background_SR_epsilon ] = 2.706\n",
      "RooStats::HistFactory::FlexibleInterpVar::signal_SR_epsilon[ paramList=(alpha_signal_norm_uncertainty) ] = 1\n",
      "RooHistFunc::signal_SR_nominal[ depList=(obs_x_SR) depList=(obs_x_SR) ] = 0\n",
      "RooProduct::signal_SR_overallNorm_x_sigma_epsilon[ mu * signal_SR_epsilon ] = 0\n",
      "RooProduct::signal_SR_overallSyst_x_Exp[ signal_SR_nominal * signal_SR_overallNorm_x_sigma_epsilon ] = 0\n",
      "\n",
      "datasets\n",
      "--------\n",
      "RooDataSet::asimovData(obs_x_SR)\n",
      "RooDataSet::obsData(obs_x_SR)\n",
      "\n",
      "embedded datasets (in pdfs and functions)\n",
      "-----------------------------------------\n",
      "RooDataHist::signal_SRnominalDHist(obs_x_SR)\n",
      "RooDataHist::background_SR_Hist_alphanominalDHist(obs_x_SR)\n",
      "RooDataHist::background_SR_Hist_alpha_0lowDHist(obs_x_SR)\n",
      "RooDataHist::background_SR_Hist_alpha_0highDHist(obs_x_SR)\n",
      "\n",
      "parameter snapshots\n",
      "-------------------\n",
      "NominalParamValues = (Lumi=1[C],nominalLumi=1[C],alpha_signal_norm_uncertainty=0,nom_alpha_signal_norm_uncertainty=0[C],mu=1,obs_x_SR=159,alpha_bkg_norm_uncertainty=0,nom_alpha_bkg_norm_uncertainty=0[C],alpha_background_shape=0,nom_alpha_background_shape=0[C],binWidth_obs_x_SR_0=0.5[C],binWidth_obs_x_SR_1=0.5[C],weightVar=0)\n",
      "ModelConfig__snapshot = (mu=1)\n",
      "snapshot_paramsVals_initial = (Lumi=1[C],nominalLumi=1[C],alpha_signal_norm_uncertainty=0,nom_alpha_signal_norm_uncertainty=0[C],mu=1,alpha_bkg_norm_uncertainty=0,nom_alpha_bkg_norm_uncertainty=0[C],alpha_background_shape=0,nom_alpha_background_shape=0[C],binWidth_obs_x_SR_0=0.5[C],binWidth_obs_x_SR_1=0.5[C])\n",
      "B_only_model__snapshot = (mu=0)\n",
      "\n",
      "named sets\n",
      "----------\n",
      "B_only_model__snapshot:(mu)\n",
      "ModelConfig_GlobalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_bkg_norm_uncertainty,nom_alpha_background_shape)\n",
      "ModelConfig_NuisParams:(alpha_background_shape,alpha_bkg_norm_uncertainty,alpha_signal_norm_uncertainty)\n",
      "ModelConfig_Observables:(obs_x_SR)\n",
      "ModelConfig_POI:(mu)\n",
      "ModelConfig__snapshot:(mu)\n",
      "coefList:(binWidth_obs_x_SR_0,binWidth_obs_x_SR_1)\n",
      "constraintTerms:(lumiConstraint,alpha_signal_norm_uncertaintyConstraint,alpha_bkg_norm_uncertaintyConstraint,alpha_background_shapeConstraint)\n",
      "globalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_bkg_norm_uncertainty,nom_alpha_background_shape)\n",
      "likelihoodTerms:(SR_model)\n",
      "obsAndWeight:(weightVar,obs_x_SR)\n",
      "observables:(obs_x_SR)\n",
      "observablesSet:(obs_x_SR)\n",
      "shapeList:(L_x_signal_SR_overallSyst_x_Exp,L_x_background_SR_overallSyst_x_HistSyst)\n",
      "\n",
      "generic objects\n",
      "---------------\n",
      "RooStats::ModelConfig::ModelConfig\n",
      "\n"
     ]
    }
   ],
   "source": [
    "bModel = sbModel.Clone()\n",
    "bModel.SetName(\"B_only_model\")\n",
    "poi.setVal(0)\n",
    "bModel.SetSnapshot(R.RooArgSet(poi))\n",
    "\n",
    "w.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cooperative-fighter",
   "metadata": {},
   "source": [
    "Asymptotic calculator\n",
    "=======================\n",
    "The Asymptotic calculator it is based on the Profile Likelihood test statistics,\n",
    "We will do:\n",
    " * create the AsymptoticCalculator class using the two models and the data set;\n",
    " * run the test of hypothesis using the GetHypoTest function.\n",
    " * Look at the result obtained as a HypoTestResult object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "coordinate-james",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:47.476864Z",
     "iopub.status.busy": "2021-04-07T01:25:47.476414Z",
     "iopub.status.idle": "2021-04-07T01:25:47.588748Z",
     "shell.execute_reply": "2021-04-07T01:25:47.589554Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[#0] PROGRESS:Eval -- AsymptoticCalculator::Initialize....\n",
      "[#0] PROGRESS:Eval -- AsymptoticCalculator::Initialize - Find  best unconditional NLL on observed data\n",
      "AsymptoticCalculator::EvaluateNLL  ........ using Minuit / Migrad with strategy  1 and tolerance 1\n",
      " **********\n",
      " **    1 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **    2 **SET NOGRAD\n",
      " **********\n",
      " PARAMETER DEFINITIONS:\n",
      "    NO.   NAME         VALUE      STEP SIZE      LIMITS\n",
      "     1 alpha_background_shape   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     2 alpha_bkg_norm_uncertainty   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     3 alpha_signal_norm_uncertainty   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     4 mu           0.00000e+00  3.00000e-01    0.00000e+00  3.00000e+00\n",
      " MINUIT WARNING IN PARAM DEF\n",
      " ============== STARTING VALUE IS AT LIMIT.\n",
      " MINUIT WARNING IN PARAMETR\n",
      " ============== VARIABLE4 IS AT ITS LOWER ALLOWED LIMIT.\n",
      " MINUIT WARNING IN PARAMETR\n",
      " ============== VARIABLE4 BROUGHT BACK INSIDE LIMITS.\n",
      " **********\n",
      " **    3 **SET ERR         0.5\n",
      " **********\n",
      " **********\n",
      " **    4 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **    5 **SET STR           1\n",
      " **********\n",
      " **********\n",
      " **    6 **MIGRAD        2000           1\n",
      " **********\n",
      " MINUIT WARNING IN MIGrad    \n",
      " ============== VARIABLE4 IS AT ITS LOWER ALLOWED LIMIT.\n",
      " MIGRAD MINIMIZATION HAS CONVERGED.\n",
      " MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.\n",
      " FCN=31.25 FROM MIGRAD    STATUS=CONVERGED      93 CALLS          94 TOTAL\n",
      "                     EDM=4.37767e-06    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  alpha_background_shape   1.09635e-01   8.77455e-01   6.86556e-04   1.19450e-03\n",
      "   2  alpha_bkg_norm_uncertainty  -1.51887e-01   7.06075e-01   5.51247e-04  -9.41809e-04\n",
      "   3  alpha_signal_norm_uncertainty   2.95562e-04   9.92555e-01   7.77531e-04   8.25712e-05\n",
      "   4  mu           2.99864e+00   2.04722e+00   6.36625e-02  -1.18746e-04\n",
      "                               ERR DEF= 0.5\n",
      "AsymptoticCalculator::EvaluateNLL -  value = 31.25\tfit time : Real time 0:00:00, CP time 0.040\n",
      "[#0] PROGRESS:Eval -- Best fitted POI value = 2.99864 +/- 2.04722\n",
      "[#0] PROGRESS:Eval -- AsymptoticCalculator: Building Asimov data Set\n",
      "[#1] INFO:InputArguments -- AsymptoticCalculator: Asimov data will be generated using fitted nuisance parameter values\n",
      "MakeAsimov: Setting poi mu to a constant value = 1\n",
      "MakeAsimov: doing a conditional fit for finding best nuisance values \n",
      " **********\n",
      " **    1 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **    2 **SET NOGRAD\n",
      " **********\n",
      " PARAMETER DEFINITIONS:\n",
      "    NO.   NAME         VALUE      STEP SIZE      LIMITS\n",
      "     1 alpha_background_shape   1.09635e-01  8.77455e-01   -5.00000e+00  5.00000e+00\n",
      "     2 alpha_bkg_norm_uncertainty  -1.51887e-01  7.06075e-01   -5.00000e+00  5.00000e+00\n",
      "     3 alpha_signal_norm_uncertainty   2.95562e-04  9.92555e-01   -5.00000e+00  5.00000e+00\n",
      " **********\n",
      " **    3 **SET ERR         0.5\n",
      " **********\n",
      " **********\n",
      " **    4 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **    5 **SET STR           1\n",
      " **********\n",
      " **********\n",
      " **    6 **MIGRAD        1500           1\n",
      " **********\n",
      " MIGRAD MINIMIZATION HAS CONVERGED.\n",
      " MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.\n",
      " FCN=31.4141 FROM MIGRAD    STATUS=CONVERGED      35 CALLS          36 TOTAL\n",
      "                     EDM=1.62946e-06    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  alpha_background_shape   1.16969e-01   8.66549e-01   6.79692e-04  -2.39926e-04\n",
      "   2  alpha_bkg_norm_uncertainty  -5.97097e-02   7.05154e-01   5.46861e-04  -8.26567e-03\n",
      "   3  alpha_signal_norm_uncertainty   9.11847e-03   9.93038e-01   7.79939e-04   2.52554e-03\n",
      "                               ERR DEF= 0.5\n",
      "fit time Real time 0:00:00, CP time 0.000\n",
      "RooDataSet::AsimovData[obs_x_SR,weight:binWeightAsimov] = 25 entries (100.401 weighted)\n",
      "Generated Asimov data for observables RooArgSet:: = (obs_x_SR)\n",
      "Generated Asimov data for global observables RooArgSet:: = (nom_alpha_signal_norm_uncertainty,nom_alpha_bkg_norm_uncertainty,nom_alpha_background_shape)\n",
      "[#0] PROGRESS:Eval -- AsymptoticCalculator::Initialize Find  best conditional NLL on ASIMOV data set for given alt POI ( mu ) = 1\n",
      "AsymptoticCalculator::EvaluateNLL  ........ using Minuit / Migrad with strategy  1 and tolerance 1\n",
      " **********\n",
      " **    7 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **    8 **SET NOGRAD\n",
      " **********\n",
      " PARAMETER DEFINITIONS:\n",
      "    NO.   NAME         VALUE      STEP SIZE      LIMITS\n",
      "     1 alpha_background_shape   1.16969e-01  8.66549e-01   -5.00000e+00  5.00000e+00\n",
      "     2 alpha_bkg_norm_uncertainty  -5.97097e-02  7.05154e-01   -5.00000e+00  5.00000e+00\n",
      "     3 alpha_signal_norm_uncertainty   9.11847e-03  9.93038e-01   -5.00000e+00  5.00000e+00\n",
      " **********\n",
      " **    9 **SET ERR         0.5\n",
      " **********\n",
      " **********\n",
      " **   10 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **   11 **SET STR           1\n",
      " **********\n",
      " **********\n",
      " **   12 **MIGRAD        1500           1\n",
      " **********\n",
      " MIGRAD MINIMIZATION HAS CONVERGED.\n",
      " FCN=31.7809 FROM MIGRAD    STATUS=CONVERGED      40 CALLS          41 TOTAL\n",
      "                     EDM=6.49628e-19    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  alpha_background_shape   1.16969e-01   9.57629e-01  -0.00000e+00   3.46438e-09\n",
      "   2  alpha_bkg_norm_uncertainty  -5.97097e-02   7.03011e-01   0.00000e+00  -6.41859e-09\n",
      "   3  alpha_signal_norm_uncertainty   9.11847e-03   9.93249e-01  -0.00000e+00   9.05622e-10\n",
      "                               ERR DEF= 0.5\n",
      "AsymptoticCalculator::EvaluateNLL -  value = 31.7809 for poi fixed at = 1\tfit time : Real time 0:00:00, CP time 0.000\n",
      "[#1] INFO:InputArguments -- AsymptotiCalculator: Minimum of POI is 0 corresponds to null  snapshot   - default configuration is  one-sided discovery formulae  \n"
     ]
    }
   ],
   "source": [
    "ac = R.RooStats.AsymptoticCalculator(data, sbModel, bModel)\n",
    "ac.SetOneSidedDiscovery(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "robust-pound",
   "metadata": {},
   "source": [
    "Get the hypo test result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "acoustic-gibson",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:47.610867Z",
     "iopub.status.busy": "2021-04-07T01:25:47.610336Z",
     "iopub.status.idle": "2021-04-07T01:25:47.722944Z",
     "shell.execute_reply": "2021-04-07T01:25:47.724018Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "[#1] INFO:Eval -- AsymptoticCalculator::GetHypoTest: - perform  an hypothesis test for  POI ( mu ) = 0\n",
      "[#0] PROGRESS:Eval -- AsymptoticCalculator::GetHypoTest -  Find  best conditional NLL on OBSERVED data set ..... \n",
      "AsymptoticCalculator::EvaluateNLL  ........ using Minuit / Migrad with strategy  1 and tolerance 1\n",
      " **********\n",
      " **   13 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **   14 **SET NOGRAD\n",
      " **********\n",
      " PARAMETER DEFINITIONS:\n",
      "    NO.   NAME         VALUE      STEP SIZE      LIMITS\n",
      "     1 alpha_background_shape   1.09635e-01  8.77455e-01   -5.00000e+00  5.00000e+00\n",
      "     2 alpha_bkg_norm_uncertainty  -1.51887e-01  7.06075e-01   -5.00000e+00  5.00000e+00\n",
      "     3 alpha_signal_norm_uncertainty   2.95562e-04  9.92555e-01   -5.00000e+00  5.00000e+00\n",
      " **********\n",
      " **   15 **SET ERR         0.5\n",
      " **********\n",
      " **********\n",
      " **   16 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **   17 **SET STR           1\n",
      " **********\n",
      " **********\n",
      " **   18 **MIGRAD        1500           1\n",
      " **********\n",
      " MIGRAD MINIMIZATION HAS CONVERGED.\n",
      " MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.\n",
      " FCN=31.6406 FROM MIGRAD    STATUS=CONVERGED      35 CALLS          36 TOTAL\n",
      "                     EDM=1.00451e-05    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  alpha_background_shape   1.21099e-01   8.58248e-01   6.75353e-04  -1.78006e-05\n",
      "   2  alpha_bkg_norm_uncertainty  -2.23946e-03   7.04058e-01   5.44613e-04  -2.24332e-02\n",
      "   3  alpha_signal_norm_uncertainty   5.93232e-11   9.93347e-01   7.82964e-04   2.94939e-10\n",
      "                               ERR DEF= 0.5\n",
      "AsymptoticCalculator::EvaluateNLL -  value = 31.6406 for poi fixed at = 0\tfit time : Real time 0:00:00, CP time 0.010\n",
      "[#0] PROGRESS:Eval -- \t OBSERVED DATA :  qmu   = 0.781164 condNLL = 31.6406 uncond 31.25\n",
      "[#0] PROGRESS:Eval -- AsymptoticCalculator::GetHypoTest -- Find  best conditional NLL on ASIMOV data set .... \n",
      "AsymptoticCalculator::EvaluateNLL  ........ using Minuit / Migrad with strategy  1 and tolerance 1\n",
      " **********\n",
      " **   19 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **   20 **SET NOGRAD\n",
      " **********\n",
      " PARAMETER DEFINITIONS:\n",
      "    NO.   NAME         VALUE      STEP SIZE      LIMITS\n",
      "     1 alpha_background_shape   1.21099e-01  8.58248e-01   -5.00000e+00  5.00000e+00\n",
      "     2 alpha_bkg_norm_uncertainty  -2.23946e-03  7.04058e-01   -5.00000e+00  5.00000e+00\n",
      "     3 alpha_signal_norm_uncertainty   5.93232e-11  9.93347e-01   -5.00000e+00  5.00000e+00\n",
      " **********\n",
      " **   21 **SET ERR         0.5\n",
      " **********\n",
      " **********\n",
      " **   22 **SET PRINT           0\n",
      " **********\n",
      " **********\n",
      " **   23 **SET STR           1\n",
      " **********\n",
      " **********\n",
      " **   24 **MIGRAD        1500           1\n",
      " **********\n",
      " MIGRAD MINIMIZATION HAS CONVERGED.\n",
      " MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.\n",
      " FCN=31.8234 FROM MIGRAD    STATUS=CONVERGED      39 CALLS          40 TOTAL\n",
      "                     EDM=7.70392e-11    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  alpha_background_shape   1.23078e-01   9.51316e-01   7.51876e-04  -2.95026e-07\n",
      "   2  alpha_bkg_norm_uncertainty  -9.72312e-03   7.02142e-01   5.53620e-04  -6.22950e-05\n",
      "   3  alpha_signal_norm_uncertainty   9.11847e-03   9.93347e-01   7.85188e-04  -2.05872e-08\n",
      "                               ERR DEF= 0.5\n",
      "AsymptoticCalculator::EvaluateNLL -  value = 31.8234 for poi fixed at = 0\tfit time : Real time 0:00:00, CP time 0.000\n",
      "[#0] PROGRESS:Eval -- \t ASIMOV data qmu_A = 0.084954 condNLL = 31.8234 uncond 31.7809\n",
      "[#0] PROGRESS:Eval -- poi = 0 qmu = 0.781164 qmu_A = 0.084954 sigma = 0  CLsplusb = 0.188393 CLb = 0.276803 CLs = 1.46929\n",
      "\n",
      "Results HypoTestAsymptotic_result: \n",
      " - Null p-value = 0.188393\n",
      " - Significance = 0.883835\n",
      " - CL_b: 0.188393\n",
      " - CL_s+b: 0.276803\n",
      " - CL_s: 1.46929\n"
     ]
    }
   ],
   "source": [
    "asResult = ac.GetHypoTest()\n",
    "asResult.Print()\n",
    "pvalue_as = asResult.NullPValue()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "historical-mauritius",
   "metadata": {},
   "source": [
    "By hand calculation\n",
    "======================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "conditional-frederick",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:47.732802Z",
     "iopub.status.busy": "2021-04-07T01:25:47.732108Z",
     "iopub.status.idle": "2021-04-07T01:25:47.842083Z",
     "shell.execute_reply": "2021-04-07T01:25:47.842658Z"
    }
   },
   "outputs": [],
   "source": [
    "w.loadSnapshot(snapshotName_init)\n",
    "sbModel = w.obj(\"ModelConfig\")\n",
    "pdf = sbModel.GetPdf()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "residential-cisco",
   "metadata": {},
   "source": [
    "Get the nuisance parameters and global observables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "operational-profile",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:47.859214Z",
     "iopub.status.busy": "2021-04-07T01:25:47.858332Z",
     "iopub.status.idle": "2021-04-07T01:25:47.968581Z",
     "shell.execute_reply": "2021-04-07T01:25:47.969172Z"
    }
   },
   "outputs": [],
   "source": [
    "constrainedParams = sbModel.GetNuisanceParameters()\n",
    "glbObs = sbModel.GetGlobalObservables()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "published-yeast",
   "metadata": {},
   "source": [
    "Create the neg-log-likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "alien-davis",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:48.021416Z",
     "iopub.status.busy": "2021-04-07T01:25:48.020948Z",
     "iopub.status.idle": "2021-04-07T01:25:48.130316Z",
     "shell.execute_reply": "2021-04-07T01:25:48.131046Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting NLL value: 31.426611090781595\n",
      "[#1] INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.\n",
      "[#1] INFO:Minization -- createNLL picked up cached consraints from workspace with 3 entries\n",
      "[#1] INFO:Minization --  Including the following constraint terms in minimization: (alpha_signal_norm_uncertaintyConstraint,alpha_bkg_norm_uncertaintyConstraint,alpha_background_shapeConstraint)\n",
      "[#1] INFO:Minization -- The following global observables have been defined: (nom_alpha_signal_norm_uncertainty,nom_alpha_bkg_norm_uncertainty,nom_alpha_background_shape)\n",
      "[#1] INFO:Eval -- RooAbsTestStatistic::initMPMode: started 2 remote server process.\n",
      "[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (signal_SR_nominal,background_SR_Hist_alphanominal,background_SR_Hist_alpha_0low,background_SR_Hist_alpha_0high)\n",
      "[#1] INFO:Minization --  The following expressions will be evaluated in cache-and-track mode: (background_SR_Hist_alpha)\n",
      "[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (signal_SR_nominal,background_SR_Hist_alphanominal,background_SR_Hist_alpha_0low,background_SR_Hist_alpha_0high)\n",
      "[#1] INFO:Minization --  The following expressions will be evaluated in cache-and-track mode: (background_SR_Hist_alpha)\n",
      "[#1] INFO:Minization -- RooNLLVar::evaluatePartition(nll_model_SR_obsData_GOF0) first = 0 last = 12 Likelihood offset now set to 226.555\n",
      "[#1] INFO:Minization -- RooNLLVar::evaluatePartition(nll_model_SR_obsData_GOF1) first = 12 last = 25 Likelihood offset now set to -197.886\n"
     ]
    }
   ],
   "source": [
    "nll_sb = pdf.createNLL(data, R.RooFit.Constrain(constrainedParams), R.RooFit.GlobalObservables(glbObs), R.RooFit.Offset(1),\n",
    "                                       R.RooFit.NumCPU(2), R.RooFit.Optimize(2))\n",
    "nllval = nll_sb.getVal()\n",
    "print(\"Starting NLL value:\", nllval)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "international-increase",
   "metadata": {},
   "source": [
    "Do the minimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "colored-catch",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:48.219438Z",
     "iopub.status.busy": "2021-04-07T01:25:48.166050Z",
     "iopub.status.idle": "2021-04-07T01:25:48.326919Z",
     "shell.execute_reply": "2021-04-07T01:25:48.327289Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " =========== Unconditinal fit =========\n",
      "\n",
      "[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_model_SR_obsData_with_constr) Summation contains a RooNLLVar, using its error level\n",
      "[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization\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 alpha_background_shape   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     2 alpha_bkg_norm_uncertainty   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     3 alpha_signal_norm_uncertainty   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     4 mu           1.00000e+00  3.00000e-01    0.00000e+00  3.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        2000           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=2.74779 FROM MIGRAD    STATUS=INITIATE       32 CALLS          33 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  alpha_background_shape   1.17167e-01   1.00000e+00   2.34356e-02   8.99701e-04\n",
      "   2  alpha_bkg_norm_uncertainty   0.00000e+00   1.00000e+00   0.00000e+00   5.83820e-01\n",
      "   3  alpha_signal_norm_uncertainty   0.00000e+00   1.00000e+00   0.00000e+00  -4.14982e-02\n",
      "   4  mu           1.00000e+00   3.00000e-01   0.00000e+00  -2.34664e-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=2.58022 FROM MIGRAD    STATUS=CONVERGED     103 CALLS         104 TOTAL\n",
      "                     EDM=3.70809e-05    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  alpha_background_shape   1.06975e-01   8.77485e-01   2.13729e-04  -1.59403e-02\n",
      "   2  alpha_bkg_norm_uncertainty  -1.54039e-01   7.10103e-01   1.71787e-04  -2.34441e-02\n",
      "   3  alpha_signal_norm_uncertainty   2.61021e-03   9.92763e-01   2.42185e-04   1.13689e-02\n",
      "   4  mu           2.99562e+00   2.89147e+00   1.83852e-02   2.57355e-04\n",
      "                               ERR DEF= 0.5\n",
      " EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  4    ERR DEF=0.5\n",
      "  7.780e-01  1.099e-03 -3.902e-04  1.221e-02 \n",
      "  1.099e-03  5.077e-01 -4.377e-03  1.370e-01 \n",
      " -3.902e-04 -4.377e-03  9.988e-01  3.628e-02 \n",
      "  1.221e-02  1.370e-01  3.628e-02  3.050e+00 \n",
      " PARAMETER  CORRELATION COEFFICIENTS  \n",
      "       NO.  GLOBAL      1      2      3      4\n",
      "        1  0.00800   1.000  0.002 -0.000  0.008\n",
      "        2  0.11041   0.002  1.000 -0.006  0.110\n",
      "        3  0.02246  -0.000 -0.006  1.000  0.021\n",
      "        4  0.11243   0.008  0.110  0.021  1.000\n"
     ]
    }
   ],
   "source": [
    "minim = R.RooMinimizer(nll_sb)\n",
    "strategy = R.Math.MinimizerOptions.DefaultStrategy()\n",
    "minim.setStrategy(strategy)\n",
    "minim.optimizeConst(2)\n",
    "minimizer = R.Math.MinimizerOptions.DefaultMinimizerType()\n",
    "algorithm = R.Math.MinimizerOptions.DefaultMinimizerAlgo()\n",
    "print(\"\\n =========== Unconditinal fit =========\\n\")\n",
    "status = minim.minimize(minimizer, algorithm)\n",
    "\n",
    "obs_nll_min = nll_sb.getVal()\n",
    "reverse = (poi.getVal() < 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "weighted-diameter",
   "metadata": {},
   "source": [
    "Fix POI to 0 (B-only model) and do the minimization again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "guided-jacksonville",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:48.355445Z",
     "iopub.status.busy": "2021-04-07T01:25:48.355036Z",
     "iopub.status.idle": "2021-04-07T01:25:48.462588Z",
     "shell.execute_reply": "2021-04-07T01:25:48.463314Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " =========== Conditinal fit =========\n",
      "\n",
      "[#1] INFO:Minization -- RooMinimizerFcn::synchronize: set of constant parameters changed, rerunning const optimizer\n",
      " **********\n",
      " **    7 **SET PRINT           1\n",
      " **********\n",
      " **********\n",
      " **    8 **SET NOGRAD\n",
      " **********\n",
      " PARAMETER DEFINITIONS:\n",
      "    NO.   NAME         VALUE      STEP SIZE      LIMITS\n",
      "     1 alpha_background_shape   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     2 alpha_bkg_norm_uncertainty   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     3 alpha_signal_norm_uncertainty   0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00\n",
      "     4 mu           0.00000e+00  1.00000e-01     no limits\n",
      " **********\n",
      " **    9 **FIX           4\n",
      " **********\n",
      " **********\n",
      " **   10 **SET ERR         0.5\n",
      " **********\n",
      " **********\n",
      " **   11 **SET PRINT           1\n",
      " **********\n",
      " **********\n",
      " **   12 **SET STR           1\n",
      " **********\n",
      " NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY\n",
      " **********\n",
      " **   13 **MIGRAD        2000           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=2.98063 FROM MIGRAD    STATUS=INITIATE       12 CALLS          13 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  alpha_background_shape   0.00000e+00   1.00000e+00   2.01358e-01  -8.13857e-01\n",
      "   2  alpha_bkg_norm_uncertainty   0.00000e+00   1.00000e+00   2.01358e-01  -2.72843e-07\n",
      "   3  alpha_signal_norm_uncertainty   0.00000e+00   1.00000e+00   2.01358e-01   0.00000e+00\n",
      "   4  mu           0.00000e+00     fixed    \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=2.97077 FROM MIGRAD    STATUS=CONVERGED      37 CALLS          38 TOTAL\n",
      "                     EDM=8.27911e-08    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  alpha_background_shape   1.21349e-01   8.58252e-01   2.22167e-04   1.66749e-03\n",
      "   2  alpha_bkg_norm_uncertainty   2.72444e-08   7.04240e-01   1.81801e-04   3.92180e-07\n",
      "   3  alpha_signal_norm_uncertainty   0.00000e+00   9.93347e-01   2.57293e-04   0.00000e+00\n",
      "   4  mu           0.00000e+00     fixed    \n",
      "                               ERR DEF= 0.5\n",
      " EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  3    ERR DEF=0.5\n",
      "  7.440e-01  5.245e-08 -2.313e-10 \n",
      "  5.245e-08  4.993e-01 -1.630e-17 \n",
      " -2.313e-10 -1.630e-17  1.000e+00 \n",
      " PARAMETER  CORRELATION COEFFICIENTS  \n",
      "       NO.  GLOBAL      1      2      3\n",
      "        1  0.00000   1.000  0.000 -0.000\n",
      "        2  0.00000   0.000  1.000 -0.000\n",
      "        3  0.00000  -0.000 -0.000  1.000\n"
     ]
    }
   ],
   "source": [
    "print(\"\\n =========== Conditinal fit =========\\n\")\n",
    "w.loadSnapshot(snapshotName_init)\n",
    "poi.setVal(0)\n",
    "poi.setConstant(1)\n",
    "\n",
    "status = minim.minimize(minimizer, algorithm)\n",
    "\n",
    "obs_nll_min_bkg = nll_sb.getVal()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "manufactured-smoke",
   "metadata": {},
   "source": [
    "The asymptotic statistic: q0 = nll(b-only, mu=0) - nll(s+b)\n",
    "Significance: Z = sqrt( 2*q0 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "assured-mining",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:48.467686Z",
     "iopub.status.busy": "2021-04-07T01:25:48.467023Z",
     "iopub.status.idle": "2021-04-07T01:25:48.577050Z",
     "shell.execute_reply": "2021-04-07T01:25:48.577621Z"
    }
   },
   "outputs": [],
   "source": [
    "obs_q0 = 2*(obs_nll_min_bkg - obs_nll_min)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "forbidden-specific",
   "metadata": {},
   "source": [
    "Check the sign: excess or deficit? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "funky-floating",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:48.584368Z",
     "iopub.status.busy": "2021-04-07T01:25:48.583617Z",
     "iopub.status.idle": "2021-04-07T01:25:48.693071Z",
     "shell.execute_reply": "2021-04-07T01:25:48.692427Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Unconditional NLL value: 31.250010639010046\n",
      "Conditional NLL value: 31.640569479530004\n",
      "==> Asymmptotic signficance:  0.8838086224064099\n"
     ]
    }
   ],
   "source": [
    "if reverse: obs_q0 = -obs_q0\n",
    "sign = 0\n",
    "if obs_q0!=0: sign = obs_q0 / fabs(obs_q0)\n",
    "obs_sig = sign*sqrt(fabs(obs_q0));\n",
    "print(\"\\nUnconditional NLL value:\", obs_nll_min)\n",
    "print(\"Conditional NLL value:\", obs_nll_min_bkg)\n",
    "print(\"==> Asymmptotic signficance: \", obs_sig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "future-event",
   "metadata": {},
   "source": [
    "The frequentist appproach\n",
    "======================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "deadly-northwest",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:48.711566Z",
     "iopub.status.busy": "2021-04-07T01:25:48.710953Z",
     "iopub.status.idle": "2021-04-07T01:25:48.820988Z",
     "shell.execute_reply": "2021-04-07T01:25:48.821723Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RooRealVar::mu = 1  L(0 - 3) \n"
     ]
    }
   ],
   "source": [
    "w.loadSnapshot(snapshotName_init)\n",
    "poi.Print()\n",
    "fc = R.RooStats.FrequentistCalculator(data, sbModel, bModel)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "violent-louisville",
   "metadata": {},
   "source": [
    "Test statistics: profile liekelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "deadly-fantasy",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:48.848835Z",
     "iopub.status.busy": "2021-04-07T01:25:48.848285Z",
     "iopub.status.idle": "2021-04-07T01:25:48.958756Z",
     "shell.execute_reply": "2021-04-07T01:25:48.959343Z"
    }
   },
   "outputs": [],
   "source": [
    "profll = R.RooStats.ProfileLikelihoodTestStat(sbModel.GetPdf())\n",
    "profll.SetOneSidedDiscovery(True)\n",
    "profll.SetVarName(\"q_{0}/2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "breathing-lawsuit",
   "metadata": {},
   "source": [
    "Need to throw toys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "destroyed-matter",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:48.987695Z",
     "iopub.status.busy": "2021-04-07T01:25:48.987140Z",
     "iopub.status.idle": "2021-04-07T01:25:49.097746Z",
     "shell.execute_reply": "2021-04-07T01:25:49.098421Z"
    }
   },
   "outputs": [],
   "source": [
    "toymcs = fc.GetTestStatSampler()\n",
    "if not sbModel.GetPdf().canBeExtended():\n",
    "    toymcs.SetNEventsPerToy(1)\n",
    "    print('\\nAdjusting for non-extended formalism\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "first-intention",
   "metadata": {},
   "source": [
    "Set the test stat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "occupational-cornell",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:49.109204Z",
     "iopub.status.busy": "2021-04-07T01:25:49.108463Z",
     "iopub.status.idle": "2021-04-07T01:25:49.219464Z",
     "shell.execute_reply": "2021-04-07T01:25:49.218715Z"
    }
   },
   "outputs": [],
   "source": [
    "toymcs.SetTestStatistic(profll)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "greek-gazette",
   "metadata": {},
   "source": [
    "Set number of toys: null toys, alt toys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "committed-jerusalem",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:49.228458Z",
     "iopub.status.busy": "2021-04-07T01:25:49.227788Z",
     "iopub.status.idle": "2021-04-07T01:25:49.337782Z",
     "shell.execute_reply": "2021-04-07T01:25:49.338384Z"
    }
   },
   "outputs": [],
   "source": [
    "fc.SetToys(2500,1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "external-engineering",
   "metadata": {},
   "source": [
    "Run the test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "authentic-sterling",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:55.870792Z",
     "iopub.status.busy": "2021-04-07T01:25:55.870444Z",
     "iopub.status.idle": "2021-04-07T01:25:55.978854Z",
     "shell.execute_reply": "2021-04-07T01:25:55.979241Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== Using the following for B_only_model ===\n",
      "Observables:             RooArgSet:: = (obs_x_SR)\n",
      "Parameters of Interest:  RooArgSet:: = (mu)\n",
      "Nuisance Parameters:     RooArgSet:: = (alpha_background_shape,alpha_bkg_norm_uncertainty,alpha_signal_norm_uncertainty)\n",
      "Global Observables:      RooArgSet:: = (nom_alpha_signal_norm_uncertainty,nom_alpha_bkg_norm_uncertainty,nom_alpha_background_shape)\n",
      "PDF:                     RooProdPdf::model_SR[ lumiConstraint * alpha_signal_norm_uncertaintyConstraint * alpha_bkg_norm_uncertaintyConstraint * alpha_background_shapeConstraint * SR_model(obs_x_SR) ] = 0.0237772/1\n",
      "Snapshot:                \n",
      "  1) 0x137dd2bf0 RooRealVar:: mu = 0  L(0 - 3)  \"mu\"\n",
      "\n",
      "\n",
      "=== Using the following for ModelConfig ===\n",
      "Observables:             RooArgSet:: = (obs_x_SR)\n",
      "Parameters of Interest:  RooArgSet:: = (mu)\n",
      "Nuisance Parameters:     RooArgSet:: = (alpha_background_shape,alpha_bkg_norm_uncertainty,alpha_signal_norm_uncertainty)\n",
      "Global Observables:      RooArgSet:: = (nom_alpha_signal_norm_uncertainty,nom_alpha_bkg_norm_uncertainty,nom_alpha_background_shape)\n",
      "PDF:                     RooProdPdf::model_SR[ lumiConstraint * alpha_signal_norm_uncertaintyConstraint * alpha_bkg_norm_uncertaintyConstraint * alpha_background_shapeConstraint * SR_model(obs_x_SR) ] = 0.0237772/1\n",
      "Snapshot:                \n",
      "  1) 0x137dd2bf0 RooRealVar:: mu = 1  L(0 - 3)  \"mu\"\n",
      "\n",
      "[#0] PROGRESS:Generation -- Test Statistic on data: 0.390558\n",
      "[#1] INFO:InputArguments -- Profiling conditional MLEs for Null.\n",
      "[#1] INFO:InputArguments -- Using a ToyMCSampler. Now configuring for Null.\n",
      "[#0] PROGRESS:Generation -- generated toys: 500 / 2500\n",
      "[#0] PROGRESS:Generation -- generated toys: 1000 / 2500\n",
      "[#0] PROGRESS:Generation -- generated toys: 1500 / 2500\n",
      "[#0] PROGRESS:Generation -- generated toys: 2000 / 2500\n",
      "[#1] INFO:InputArguments -- Profiling conditional MLEs for Alt.\n",
      "[#1] INFO:InputArguments -- Using a ToyMCSampler. Now configuring for Alt.\n",
      "[#0] PROGRESS:Generation -- generated toys: 500 / 1000\n",
      "\n",
      "Results HypoTestCalculator_result: \n",
      " - Null p-value = 0.1792 +/- 0.00767039\n",
      " - Significance = 0.918418 +/- 0.0293136 sigma\n",
      " - Number of Alt toys: 1000\n",
      " - Number of Null toys: 2500\n",
      " - Test statistic evaluated on data: 0.390558\n",
      " - CL_b: 0.1792 +/- 0.00767039\n",
      " - CL_s+b: 0.286 +/- 0.01429\n",
      " - CL_s: 1.59598 +/- 0.105004\n",
      "\n",
      "Results HypoTestCalculator_result: \n",
      " - Null p-value = 0.1792 +/- 0.00767039\n",
      " - Significance = 0.918418 +/- 0.0293136 sigma\n",
      " - Number of Alt toys: 1000\n",
      " - Number of Null toys: 2500\n",
      " - Test statistic evaluated on data: 0.390558\n",
      " - CL_b: 0.1792 +/- 0.00767039\n",
      " - CL_s+b: 0.286 +/- 0.01429\n",
      " - CL_s: 1.59598 +/- 0.105004\n"
     ]
    }
   ],
   "source": [
    "fqResult = fc.GetHypoTest()\n",
    "fqResult.Print()\n",
    "fqResult.GetNullDistribution().SetTitle(\"b only\")\n",
    "fqResult.GetAltDistribution().SetTitle(\"s+b\")\n",
    "fqResult.Print()\n",
    "pvalue_fq = fqResult.NullPValue()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "multiple-journey",
   "metadata": {},
   "source": [
    "Plot the distributions of the test statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "dated-pension",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:56.008138Z",
     "iopub.status.busy": "2021-04-07T01:25:56.007735Z",
     "iopub.status.idle": "2021-04-07T01:25:56.119302Z",
     "shell.execute_reply": "2021-04-07T01:25:56.118751Z"
    }
   },
   "outputs": [],
   "source": [
    "c = R.TCanvas()\n",
    "plot = R.RooStats.HypoTestPlot(fqResult)\n",
    "plot.SetLogYaxis(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eight-fundamental",
   "metadata": {},
   "source": [
    "add chi2 to plot, to check the asymptotic behavior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "strong-karen",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:56.241048Z",
     "iopub.status.busy": "2021-04-07T01:25:56.240723Z",
     "iopub.status.idle": "2021-04-07T01:25:56.393329Z",
     "shell.execute_reply": "2021-04-07T01:25:56.393072Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Info in <TCanvas::Print>: png file test_p0_1.png has been created\n",
      "Info in <TCanvas::SaveAs>: ROOT file test_p0_1.root has been created\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAHYCAIAAAApvgy/AAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dYZarqBaAUXir55XckSU1siQj4/0gEgQOiiKifnvd1V1FCKKVyBEQtTFGAQAApPxv7woAAIB+ESgAAAARgQIAABARKAAAABGBAgAAEBEoAAAAEYECAAAQESgAAAARgQIAABARKAAAABGBAgAAEBEoAAAAEYECAAAQESgAAAARgQIAABARKAAAABGBAgAAEBEoAAAAEYECAAAQESgAAAARgQIAABARKAAAANF/e1fgVLTWe1cBAHAVxpgGWyFQqMz92bTWC/6EC97VbEO8i3dVfJfN33MNeddeG+JdM9/V7NKUoQcAACAiUKhMD/yfGZIAMIc9XTyfz70rgl7ErUn7ZoWhh8pWDj0AAOBkRugYegAAAPsjUNiV1uE/AAf3fr+VUs/n0/6wppz1hQDr0T1eU9lwQyos0K1ud+kfYze+8x2NNXvU89GIe4Mfj8fMOQfv9/vfv39Kqdvt9vl8XPrtdsuECz0fjcaudiia7S89Cnsz5vcPwJHd73f38+12sz/8/f3FzbydjCY1/36UYH/1SwYaI1AAgDpsA3+73Ywx7/fbXe0tGD54vV7GGFdCEDoALREoAEBNtgPADjfYxr70dsfb7ea6EF6vl/2BeyaxFwIFAD2KZ/r6U34zr273b9Lj8bA/fD6fv7+/YFGE9/t9H9iU5/Npfw2CAH+ggUEH7I51FPrCPAXnUpOSJnE0xjo9Grbhfz6f/kjB39+fGm6CiCcfrN8onw2HQ7ERehS64F+vcJskoMbTfDv5N8f9frezE16vlz+f0b50G7j89teg28Cf08DtkdjdtW4m2dqS2yONcWGBUVoppZVRinsgcHLnu5PN3dzo3w/p7pYMdtamv14vP0RwJfiFSCUA3B55OZwEgENzTb6bneDaeDd3YT5XyOISgFoIFACgDneHgm/+gkuOPzZhf+WWB+zobL1/+1oz9GDMN4WhB1zB+YYeHDurwN3jkMmTv6Ph/X7bVZwr1w9n0exLdNrv6i4IFICZThwoAG0wRwEAAOyPQAEAAIgIFAAAgIhAAQAAiFjCGRtirRgAODp6FAAAgIhAAQAAiAgUAKAOu+hy+wdDP59Pf71noC4CBQCoqcrDo4vwhElsikABAACIuOtBKaXsaurSmup2xXU/JwDk2fNG5nEPyfxKqfgt/oMh5hTrzlfBM6zjRGAWA2OUUrfbLflS8Bg3pdTr9cqUU7RV871v0Hzfp5RR6vfr8fExg+SUnwr3gQ/OG5mThjEmfoS0fzpyT6SMs7k8bnNBNZKFVN1j7KnZl4ihh1x8/Xw+7XDj6/V6vV72i/rv379mdQNwRME0hcxJ4/1+//39xW+Pz0txNuncFV/e0CeKVdrEIx16PB7+1ynZo2Bf8q8G7FukqLzseNKjgAs75aciPp+4k0zmpBGcZ9xbbIrrDPCPWJAS9CjE/QeZsxyOq9mX6Lo9Cu/3Oz852Ybet9vND9ttYhzaA4DjruDjH5Iej4c7z0hviYceJHFRQTpQ5NKBgo2VpK+fP4HIkb6BACrTurt/KyQvS/ITDNecZGwfg92oG25g3AHLXDdQmGS/Y8TgAA7HxQTuZop44gIwE4FCMft9o0cB2JYx3f1bIdlOu+sQ/1q/yn2MfslBvwJQikChMj2fUiy4CpySa6fjH5I+n098Y8LK7kx/9GF9adhaQdvxa0MaYcGlhaRvnZl/2cHC7MBJfT6f4DwuXdC/Xi9782RwC2UwjXqB5/PpymTcoX8FbcegWaxAj4JI+mq1X8gdwLEEZ49MG3C/3/0bIN3b1w9uxrdrAcvQozDBrpYap9OPByDgBwTzl3C+3+/GGCm/fTWzISVPmbrdbkzKxnoECiLbcRd8A5nDCGCO+U95WJZ/klsqhnEHrESgILJfWjvJKJicnFn5RBo0WjD+BADL+Ccixh0OreWkRQmBQs7j8fj7+/v3758dRHQReuaLR0AAoB/+mo84okyb0iyGIFDIeT6fNjjwJyTHM48AoCtcsaAizedpDtuFMDmIqHXJ8bTBoDHD/78pWpnvr8fnAl4+ZgiUfVkARJp9ifiu1kSgECBQgIRAAVip2ZeIdRQAAICIOQqVcdcDAKAW7no4IQICAEAtPdz1wNADAGDC+gcR3e/36o8ysgVy/+fWCBQAoI73+22bropLuK5vCxe00PYtjTfaAIHFMgw9AEDvGj+Lzm7O32hmOdqZ7LI0KwvBLggUAKAO1xDGD5Ozz3xSwrqu7tXFT3zIlCAFGe4tUq18yVqVVrtKIZNP25IKzIQp64//yRnUU3Y8lTLfiSrm+z6ljFK/X4+Pjxkk5/tUxM9esunxWq6Px8N/Y3xafr1eyTfa9FhcgtvEzPTMS3aj8Xc53t/8Rl1+V0Ly+dqZIxzkd50c7l3JZXOlgzn5lv6pVl+is31X9xV/4OJvl5+bQAGXdb5PRdA5b1svvx3yW9ZMg+1/Zfy33G63ZAMW5Alau2QJ0vMkM2/JV6xoo67myTpIsUJm7Xz3FimDfXtQqznHvwdSDVvWs6PDcQJlfzYCBVzYKT8VrjFzLXqmdfR/ja+J40hC2micIShBqkPRRpPVzpQZbzRIcXFVpp4+96o7tq4OQUzmCnSbCN4S73jm79WzZvVkjgIAbMXODwj6Ev7+/tR4HsPn87Gj48/n0yxdi2V+CX6G9/u95jnUi6sdT494Pp+Zmrgj6Q7a6/XyH9d3v9+D/bLHOV959xZbn8aTRg+jTTxyEWXHkx4FXNjkp2KnM2LO5E7FV6iZ0jKjAP4khsmtJ+9H8EtIXiXntxtvNEhZsNEgxf6cn5TgizcRFxKMLwSHOrlR6S0za7WvZvWkRwEANidNC7Dz7Z/Pp38tay+FZ17oP59Pe00flODf1BDwlzewTf7kxff6jTbgOhhcx0N+v97vd/CW9/tNp0JCm3jkIsqOJz0KuLDJT8VOZ8ScyZ2SehRmXjf7V7fxbLuZJfjX+jYxuJSPx/IXzFEo3WickrwJ4na7SccqPpKu2jYxnpEwOUchnvfAHIUkVmYE0KM2Z8Ai8yvvrqptw+NfpPpLFj6fT/uz7TywV7RFR8mtBelKWDDhoPQtbqP2qn3ZRt1UA/fef//+fT4f6YLeHUk/f34Tpd0kijkKkr2/eqdSdpzpUcCFnfJTEdwMGae4nydvj0zO7U/eHimVEN9NoKJFESbfIt0eKZUQX6yraDZGvubBLITksU1uNJNBulEis9xkP0spSDUMDua2dWizmYso+7MRKODCzvqpiE/i+QWXks2bn8FvzJKtV7KEZBe9KzloIP0M0kbj73LpRhssuBRMBAnGJoIM8V7cbre4zJ6pVl8ibbocCzworUuOp51PZMzw/2+KVub76/G5OVN8zBAo+7IcSnKN4fzCw3bEQcowuWjxZAlqmH7oXpq5EPL6jaqpoY05G4o3urjacYbSCvSj2ZfotN/VXRAoBAgUIDlxoAC00exLxGRGAAAgIlAAAAAiAgUAACAiUAAAACICBQBoxy1+vHdFgLmYeFyTv4J6IHGcuesBF3bNux6CU8QFjwBKZZoV1eojRI9CZdKCFXvXC8DObC+CXevQrvxzuBv30V5mHaRmdSBQAIB27MI+8dCDfe6DHzrYQQqbWHGo4vl8BqW5Z08ASVfs/dsOCy4FGHqA5JpDD879fv98Po/Hw7XZ9svyer1srJBsudcfsefzaR+V5BcVVwaHwMqMh0SgECBQgOTKgYJtmJX3vbApt9vNrk/smvPb7WYfKWnzuwylXLHu6Yj+wX+/3/ZJjJf9ixwUgcIhESgECBQguWygYL8UQZMfdCfEkYT/VfLf6D+wIEiJN+oLDr7NQKfCsbCEMwAcjJ1V4Dfk7/fbb3pde+zncRlcG+/6D1ye4LmI/wbujfZXqb/hNpBqbp+aaLsxgFDJoyYxoex48phpXNj5PhXxw5rN8BWwPwcPd1bDs4xt+x08ezp4onTwjGZX1OPxSG5XEj/rOahq8knW6FOzLxE9CgBQwfP5tGfV2+1mOwlsD4FryKXLfdt/4A8ZxI+oDvoYXCfE39+f+9mPGBZbNgcC50agAAAVuAb7+Xy+32/XuvvTCIILtTmt8vP5tEMMwbQGFxbYrTwej5WrMtgohEABsf/2rgAAJOQXMdzr1QzXTtupiDZumLzKz09CdLMa4/sd7vf77XZzdzEwCRHboUehMi3Yu14A2rFNuL25cXEhWmvXW5C80HdRgqIn4LykNqVls0KPQmXzrz8AnJVtwuc03v7AhB9VuGbA3TMZCLoQ/v37x8nnlDJ/1maxAj0KAHoUDOd38uoc7pkOC97rl2DZGywdl+huZXR3QPiv2ivOovGIeE4l8DX39gjMUHY8uT0SF3biT4X9zNtbH+ewIYWff/KM7X61dzMGv7qJEf4tl/62kgffLwGH0OxLRI8CAFRjL+Ifj4c/gSDPXsTPz+/6CdwcCBcZuCWYSsWLPgHORVdR3QhLOAdYwhmSsy7hbD/zxhj7yEfb7gbzD6R31V1BWWs9v0CWcD4ilnAGgIOx0YCdNGCXVbbpk61v9RWU7TzKmd0DbtIlUQKSzhnU74UehQA9CpCcskfBdSf4v6p5V+rBc6HWsE+DpDvh9OhRAIAjeb/ft9vNX2HJ9hO4FZ3z7BurLIdwv9+NMUWt/sxK4ppOGNTviB6FAD0KkJyyRwFoiR4FAACwPwIFAAAgYglnAPvgGSjAIRAoANgBExSAlZqF2gQKlUl/OU6LAIBSPXS8EShURkAAAKgl06bw9EgAALA/AgUAACAiUAAAACICBQAAICJQAAAAIgIFAAAgIlAAAAAiAgUAACAiUAAAACICBQAAICJQ6JBWSvewvjcAAAQKAABARKAAAABEBAoAAEBEoAAAAET/7V0BxIxSSn4EOQAA7dCjAAAARAQKAABAxNBDZdL6B4axBABAoR7W1CFQqIyAAABQS6ZNaRZDMPQAAABEBAoAAEDE0EOnXJ8SYxkAgB0RKHTHKK2U2n/6CgAADD30r4cprwCAy6JHoSfGKKWGwID4AACwP3oUOsXUBABADwgUAACAiKGHfn27FOhaAADshx4FAAAgIlAAAAAiAgUAACAiUAAAACICBQAAICJQAAAAIm6P7JpW35UaWX8JALALehQAAICIQKFTwaOgeDQUAGAXBAoAAEDEHIXu/GYjaGV4iCQAYFf0KAAAABGBwmEwTQEA0B6BAgAAEBEo9I4VFAAAO2Iy4yEY5U9yBACgFXoUlng+n8/nc+9aAACwOU3P9gJa69vt9n6/4/SC42knJxoz/L88Q/fcBEw+ZgBQV1mLswJDD8Xu93vT7Wn9/SAMNz1oRiIAAK0QKMz1fD7f7/fn89m7IgAAtEOgMNcOUULcacBSCgCAtpjMONf7/TbGGGMej8eO1dDf/2itlf8PAIAtECgAAAARQw+VlS60vGJK4ne+K90JAHB0PS/ST6BQWfHtkQCAy1two2Oz2IKhh4PhpkgAQEsECgfWc1cVAOAcCBQAAICIOQqVSVf5FRfaNMbQlwAAV9DD2Z5AobI2K2/z6AQAuILM2b5ZDEGgcEDhh4OgAQCwFQKF04ljTLofAABL8Zjpmio/ZnpeCaOHSSZ7ovb7E/OYaQDYCI+Zxgruo9PBLBgAwKERKFTW4K6HeJtKKa25ageAs+GuhxOitQYA1NLDXQ8suHQePQSeAICTIVAAAAAiAoUTYLADALAV5igcnlGMOAAAtkKPAgAAENGjUFnT2yONUYm1ErRN4v4LADi6HmapEyhU1k/zrLXupSoAgEV6uD2SQOHwgk+R/8nRTHQEAKzDHIXzITYAAFRDoHBKxAoAgDoYejin77BWpREsnlwNAJdFoHB+K5/13MGUWwDAbggUKtvj6ZFzrXl4OU+uBoD2uD3yhHoICDK+90xWrWTyY9z3YQCAY+jh9kgmM57fpk12B8EuAGBD9Cic2hCKmo1jTz/kJXQAgDMhUDin1H0Kv/sguIsBADATQw+Xo7//4cIfADCNHoWzifsG5LEAbYxhpAAAkEGgcGlaH+BxEAyUAMCOGHq4nGM1snR4AMC+6FGorOcFlzzmWHMUWO4JwDWx4NIJdRYQiFw99XAfxFFqDgDX0cOCSwQK+HLzFcy4s8EopQ82XjHCwpEAsAaBAnz6ZA1oB512AHBsBApXZ4yJ+69+AxE9tbSL68LCkQCwGIHCVdgGMtlhkIwV5pS2piZRHSYyaK2MoZkHgNYIFKBUNLfR/n+L6Y2TLX0mSmiwdQBAgEDh/EZNrI5SZN9nUmdKq1GlOR0MtaIE+iQAoBSBAvbgNdcm+H+UsWKUAAAoxcqMGDHbL91oZq/1RJQAALujRwFJv7se9PC/ylMWplZbrNi6EyUAwGIECpUdYwnnoJLz6qa1Tu+FDhdoGv1/qiLJfPQBAIBiCedT6isgqC0RK6z4EH+jgdkrJ3bwfQGApljCGc3NbubdugthU65TN06UL2lU2mcw/56Frb87k+WfOlYEcDkECpjmxwrfH7RW67pPFkQJM3MufnV9+QBwPgQKSPBaZTP637izK7nQgnst+avxU+bOjSi+RpdGLmqto8DICIDr4PZI7GReo32UOySZfQngrOhRQJnvXZPDr/FIhL2PIfFMKfuqP8NhKlYgSgCA3REoYEmnefRsiO/P0qwFreyDp6JX5ViBKAEAesDQAxqxMwTCRCFEIUoAgE7Qo3Btq1u5YCRiajvGeE3r1usoECUAwHoECqhgPOKgVaZRt2MNo4BhltJ7FogSAKAKhh7QnDFK643WUcjnr3UHI1ECgOugRwGtfZ8vtWIdhbi9zywUGZQzPORqOn9GPgpZWTgAdIVAAdUl5i3o4f7JX6ufjxW8xtbec7myvW0zEsGaSwDOh0ChsmM8PXJLbt6iio5GeP+k0LTGN0cEhfut/synP7j85U+lKCi/3gOzAEApnh55StcJCOrIX4a7uMBTcb5ClT6AredDALgynh6Jy/E/2d+RCClramyibpSw3tZRCADsjkABLeRWcsy0puPG9lhRAgCcA4ECFml2sTzECs2ihMlbKlaWDwDHwjoKaMoY4w+5aa21UnrGYFvRPICKUULd8gHgcOhRQKHVLaG/jkI4EqHkKQvpctJNuzGJNl/PnmmamkOZyJCsDwCcDD0KaCpcR2FG/olyZr9tl7UgAeDoCBTQRHKewfcKf7Sywnckwo5KrOzhTy5uEFZqRflL8wPAgTD0gCa+z3eImlNjhx+kQQH9CyOWrdUo91tMjFzMxjoKAM6NHgW0oHW6zdZa2fZ62+2m0rlDEgDmoEcBLSSf7+Cu6cePnR49KmJ8N8TQ8ZAKOGJha93q+RF1o4SZt2sCwEboUUBb0roI4/4GMxi/uazjocFToOanL94KAOyLHgU0J89XsEMF80sSW+WgffUfITV57+OK7W404lD0BCwAqIseBbQmzVdQyo4t6DDBCx209b0tongVpqIHMPS2jgJRAoBdECigqcl1FIK5jTZ/chwiORIxaq2Tz4mYvX5DV3dIMjsSwF4IFNDEjHkJv7xKnscwtZFk+VI5R1lHgSgBwI6Yo4AmptZRmJnfrv3st++j51ar1OyE5HZ/L6bT52uwjgJRAoAdESigBf/5DkG6+o41hOlSfmOMsY1w1H7GMyGlckorn3kpM1+h/ayCBVtcE4WsD7MA9I+hB7SQnB/wuwVBhY1qPn+G9v8J8xJqtd+9rbbUOC5hciVwEfQoVCY9LnnugwtP7/fcyNR8AumeyTh/mCPJ62IIypnXyGU6IxLpQy3cCIixMVBbQa1ajowQOgDVSW1KSwQKlREQTCtdR0GeZ+BafaO+7VXuS5Wdr5C0IEroSoORET7vwKYyp6xmMQRDD2itaB2FTH6hAyK17oJbfqFwHYXikQUzPK9yg/az9JzQ28gIgIMiUEBTy9ZRSObPtIK5dRfWraMw3cpuM4mx1h2bRAkAShEooIla6yhk2+C56y4sXUdh7nh/7ViBKAHAjggU0MR3fkAiPf2MaTl/pnlPjkSMnx/5G4VQ04FEUHxJetUG+RBRQpeTNABUQKCAFjLzDIKxhsn8ueZdyG+Gx1GOXiqs/6o+hob2ihLoqwDOirse0MJonoHX828T7G2L/p2EmfwTm0nml25zHNrzYQGndOnidrVOrAXp3Ysx3siomjN9QyjtflVaGb9AP8Yywf+D+iwyp/KniRLW/KWAE6NHAW2Z9PyDZL9CJv/vjcL8gOT8huCeiKioqadM5TY8lX+BqX6J9BGbV585fR4z85yjNWXoBJDQo4DmxusZGL8fQV5zKb+Owqz8yVihaHaknx7V3r4QbNf+5spZ0Bpp19OhdbhBN9VCWBihyjoKcvEqrM/xLf4zASdGjwJam5h/UJIuXisXlBP0MWhvymP5vISSetZUHiXMd5o+g5mutr/AJAIFNCWti1CariZn7c0uRw3BQtQJUTASMb1dZec0aKO00t6/oBThpelr3GCqZr17Hy7VahIlADECBTQhzxtYkJ7ZSGn5U7Ue3VCZaEWm6ulHCZMby7xyiDskT4DjA8QIFNCEya2LUJqe3kSV8sORCL/4wvKThdjCjbzGs31p+GfnBhTNRtwrSmBcHzgrAgW0UHFeQq6JXV3+UFJ6Deigg2GifEH8TO28ZP5Mq79XlMC1OHBWBApoocq8hLA1isfyy8sf3WEYzSHI3Ew5WtsxOy8hMRch2y+SSMwfhylECQDWIFBAW0vnJcxtjQrLEddv8N+aGZBYU89sOVJ+aZ6EJHePRonSkQ4Ap8E6CmjOTK5zMJXufk42U4Xlf2MFY9TQHI76A7z8w3IGbqFE46/RZILyR/GIW0ghrFSCEdaRjEqdzJ+frzA/VsjkP2WUwHwLwEeggNa0/rapqTv6dkhXKt3WyvmNMYnGRNsooWi7XjFm9Gu4pnWG1t/8cfqc+Qpz7sbYLhqIt146qrLm7VKZRVEUcHoMPaCp5fMStkn/qVGOneUYzHac2K4gGBOZbs5VHJ2ksvW0jsLKxniLtpyRFCBGoIAmqsxL2CB9ZJjbOL8ck1imye1xtLaj8Ext490saby7I23+7EDE6C3epmdFCf4bjUn8U01GFuL6l75d1asnUQIQI1BAE6bmOgrbpYvrFkyUk76d0r+h0i8/yiaTtltvHYVafQ979dXTBwBsjUABLei66yhslq7kdQvmlCOECzbnL1ZIl5N+l7hdqXXcK0rYq7UmSgC2RqCAFrqal+Cnp6/v185XkNZrmip/lGnGduc5cZQAoAEChSXe7/feVTisDuYljIbzpbGAGuUH4YIxSrmnU47fOqt8W2FhHYUGIxFF5Zfaa+QCwCQChWLv9/vfv3971+LI+puXkG6mKpVvhDmP37UdVXp9yYKRhSC/Hj2CMhkGVexjqBUldNUnQdQC+AgUCrzf7+fzSZSwUofzEvabD/H7N6ueKQvybzoSUarDKKGr+gC7Sy1gB8H9fv98Pvbn5HHTyQUBJfbMPqyUc8q/g+tid4clHI/vIF25gz+scWTMVtsdVjtIN+rGmFw5KvqgBB+deLvRZ2tNlBBvLV+dSZMFzn+7C5MWf4+CEk78rcRplLU4azZ04kDh+Xy6/8be77ebaiDlkYr9+/sjUJjjFyio754mxtd3Tf9WL2rLt9uu9oLM8VyFEePaK1fOnJY52O44w8q+hG4DBVvPld+jlZUB2msWKJx5Cee/v7/b7ZZ8ye8bsDlfr9f9fm9Uswsyq5/vsE160FJ/5zZutt1h7qRLFzsYvqcAqfxvH0Mkn3/LEYe9xvUZKQC2dto5CplW//l82ijh9Xq9Xq/H46GUYubBpnqbl5CdN9C6/GGyYxRVDEs1xR0PmdZdqs/WUcJerTVRArC1s/UoPJ/P9/vt9xbE/v7+lFKuC+F+v9u3PAfx3Y/P55P+hjW+Z/Px4P/vx/3S089e2my7+fLts6OSV+Zaa78UoxLdENP1kfIHqXH52daYa3rg3M4WKExGCXY6wu128xt+ey/D39+fDRS2reLFDW1Y2LrslB6NBWy7XRsF2LbYaK3i+RDDhEatVRwLjNruX6q83eHHWCZKSJJeKYoSMsMTDUYukptoGeLEFSDAwiGcbejh/X7bVW7sgEIyg4oGJtyvrKTUgrnWOgqZ7WbWejJGhcs7jouwwYb6zWNIlSOYaN3N+MFQ2fy1ooStm8yiQGSLqGWvORzAemfrUZhk+xu2G0fIzGNPuuAVhdbfPvaw5dspXan4hsittjvqzNBD/0HiCdG/cvw7S5Ny2xXzC+npLdRsxRMBU8ORi/D4p/h9P1tUYNPycVylbUdLZ+tRWMzeHzGnR+H5fGbuSDHzXTdKUEqFc+72Sv/poz759KCPwZvsGDY+ZfMSgvSGJ6xaUUJplauMpCzAfA5ICtqOXxvSCIECmrBD7sI4+o7pI+Y7FtBhPcdDEMZ4PQ2WHyusihKUfHxqqxglFJVTZSRlGaIEHBGBwgi3NmwlM47eU7ptbPupTz49rPzvXsrkwgwl0YO0XdmCa/reogQASZcLFKQlmPL3SmAl3dl6CVK6EuY29lbP7/yGqF/By6CD/KV9DEUD6Tu21kQJwNYuFyhY0lyE9T0KWrCy2KPrcLzfpmfuO+iqnpn0zICl7VtY3scwr1E90DU9UQIOR2pTWjYrlwsU7DIJQaBQ8a7I3Wed9K6z8f78PYr91HNe+viTFnYqzOpjUPK6C0m1ooTLx9JAGpMZd2D7DD6fjx8c2OhBWnoBlXU23r/XOgobpf9OIlFW/1qkVh9A6cmqaARkLw1mcQIHcrlAQQ0Bwb9//+wDJN2jH1iTsY3exvuV6qs+FdONEm/BTT4/QrJpK95hlHCp/QUmXW7BJeU9D8J/ENTr9dqxStcRjrsPp8y90n++be3O9amcru0rLjEzGGFc/pZNWa1NlV6j1xpJKUWUgCNq9DTrPtkuhPv9XuuuyLKngytlXGkAABwNSURBVOvvWXz4f5Uq9MU1RUZ99zTcTa33Tf9Wz/zSdfz8hQ7quTB9/NnS36UW0o2kcd0SQ8yhtTLqV0LwQY0/t/lPcjL//LdnChz2a+7Wk/kn926lrcvHBZW1OGs2dOVAobrMNNTEcb5UoCBdx6s5k+83TI9P30b1WM+F6ckPpNDBMM4ybn5TJSRbPjU7UIirvyxQUF7zP2frbrsECjiE/N0NbVrwK85R2NTu01P7pPsbv5daweTcxt7quWC/YnYhhuRL8ZrQeVLQUit/xl7bBdrgrgdcxWgcffA7a++XfoJ1FHLpxihjtBo9CjLIP70Mg/o+qVKLUcVhWmuiBGABAgW0NbRhYWuxU7q9cO6/ntXS3Y5H6fkrFKn/s1aUMLvzAkBrBApoznSx3oCTHGvosJ7V0rOtuxm6OlMv6m8Hg3e4Sq/RS+uzi62jFqIiHAuBAlrrcfy+p/o0SJ95TZ8ZB624fGyHUcKm9eltf4FJBApoqqPx+0Af9Wmwv9PX9GE50iQG27mwas35Wq1maRWk/EQJQIzbI2vi9sjAUdZRcOzdcyr+c+xdzyrp7iOnUi9+f1K/mwiDcrQ4SGPflDieqWyp7Xobz+SPq+1VdlasYIYbIzPhYlB4vMWZ1ZNqmyl/U2sqP6e06gWe8pS4QA+3R15xZcZNEXilfcfLo4Ozd3rwHfzObeyvnlXSw0SllHSNK5WjlLGzOpRS40mgWgcpE1+Ein0JReUs3m69wZYd1K189UNx6GO7tUyb0uwBkgw9oIUOx+mHtlD5/5TqtJ510iNSq5ksx+U0ShulzTgWiPspbIoxLr82SiutldbhdrVWWvvZ5rflv1qN/47BX9b+87ebzD9nc8e9FjDp+2SXl6aqHpC61UNFBApoIdkG7zVOf930sen0qJyYDRfSyzDoYRJD/K4gSqhq+f6uK/8oah3vox8HzEeggLaGtic8y5DeJn2QOcvnyokvw0d55bskEkUHGepcThIl5FWs/6GPA4oQKKA509m6AldLL20tsuWnstsuhmgownYvqOT6VnVanVpRQuaau1Y9d3H0KAd7IVBAa52O318pvah1nCg/lX/obihZhmF141kxSti0Nd2xtT5ElLBXFIUM7nqoLPcMX8Tj38NhIb1NuhJai28WHbx5RvnJcjz2N61c6S6nHuWJyi9ClHAOevadrtfR7NaGDAKFyggI0nTq/v6hbSC9TboSrC1fKmf8VhsPJM962o8VvE3NdJQoQR3kmr5UrYaMKCqJ2yNxGaazcfpLpkepSkln59LyZ5/ljcrdJWFnMBSd/A4UJZxSxePG8e8WgQJaqDi+Tvry9Eim1SwtXywnRWuVmcGgSi6VdokSOugM7gLR1UUQKKAF449zD8Lxb9K3Th+bTp9XvprTWgv5bbAgtDXf+yRK+1cbRAmHbh0ZKUApAgW0NbQx6fFv0rdOH2TO8kXlSxblr/CwSqKEvIr1P/RxQBEmM6I5Y8e5o9MM6W3Ss9f6Cdnyk+9It0bZWOHXxzD8HgQHibskMtuNNmS8sifqmT0kVVpHW75ROjEdY8vm9+hRDvZCjwJaazQeT/qc+QQzrtQnyk/lF1vfTDOVqIk4KJGsdCZKkDZYpY+hlC3flE3ZrOMQUUKtkRFURKCAptqNx5M+fz6BMVqNH500v/yxBem/ao3LH5LNnEkMifKN+e3X6nrWkhiRMTwK6Yc+jz4RKFSmBXvXa2+bjruTPi9dsrZ8qZzxW+dvN5k/f5dEspehVjRTUVj+Kc4MtXaCKCFJalNaNivMUaiMBZfSTGfj9JdMDxO9F9eWX2tWYKae+jvdUcoQTmIorOc+rZTdnSOfNGz1q7RZRz4MG2LBJVxFxfF10penS2qUL7W+0gaL8iuV3G52EkNhPUujnGpqNbN7OHiQg7kIFNCC6Wmc/rrpY78s68pXc0YWSvMnpdpUM4hK078pDNny6/SFrNE2VmCkAKUIFNDWcE5cOy5O+rL0Qa3yJWvzS+RME0s92o1EI7v7RwlWqya3Yv2JEq6DQAHNGTvOTfpO6dn5AaXpidIXzD+o1OpoO41BGpD4ZdNe/nQ5klr13AV9AFiGQAGtVR53J73ifIIF5afKyZUvWd14Btv1hiTSW42GJtLlVLdja32IKGGvKAoZBApoqv64O+l15xOUli+VMzv9V61x+UUmt2uU8gcmgvIn+xhq4Zo+j+PTJwIFNLHpuDvp89Ila8uXyhm/df52pfxz6y+n//oYovxD18K2t6eH2z3F5XOtnSBK6BaBApownY3TXzI9Sv29uLb8WrMCM0UI9Re3O1XPKs+gWmsqhutfxdadKKFbBApooeL4OunL0yU1yp8YWViXXymxTRV3a149XR9DXJfkDIbpepY6cqxAH8BFsDJjZdK1yMVXbDT+OPdwKMLxb9K3Tne/B9aVr5RY/mhSgJQ/yFIeK0zkl45DarvJ0rXWNpTQ3yUfbf7UG5Z9zdvGCowUHEsPTwCgR6EyI9i7Xt0Yzolrx8VJX5YuWVp+SDqpTbWFBfXMCjfibTcbJQx5lLjco40L7DyG+u16q1NExdads1obUpvSslkhUEBzpsfx+0ukG+9BhSZ6buGi8qPUaFv+FuX8VU566VZQLjoR/fwOiBQwKDX0I2i18MGPe10i0geAZQgU0FrlcXfSO5z3UJpfVWg8MyMLxfmNMd/nRYiXbq6DIfWI6yX1bOAQUUIHHe0IESigqXCcm/Sjp0vm5f9lifIXyZS/LF2rb6zgGPefxLvmVpdr+jyOT58IFNBElfF10vtMH4TN5Yz8yfJV4bVv9Sjh+6tK9ItIPQxa/W6TyFQ1fOcpLp9r7QRRQrcIFNCE6WacnvQt0jPzA0rzS4TmSGxdytd7kOqfzuz+k9jE7JZTLv8oKrbuRAnd4vZItKC1MvacOD4ZkH6idBVbkl9I0coYbzqAa17F1qVku9LL2k5UEGMFv76jPEGskJugnqqnt/WcUf9HKnO+6Y3fMllgXIK0iXzhcwQlBG8vOjiThU/mr2vOse0qbKJHAS2YPsfXSa+RroST2vL8gnjewASp/tF2k/n97Sbz+/UxUveCX4wSVnBKbbeoo6F65upRQl0ry9+3E+eIXUj0KKCt4fopPMuQfsx0JVibfwgmop6DiXkMYZ+EmZHfeN0BXj1H+fX4xeFX5Q7NN8UMxeZaA7uCU7i3qbe4OgTHZ7L/QNr+dLQUVWpZVKfGlS8SbFd6e7C/pZtbXL0q4g6SfeuTR48CmjOdja+TviI9SlVKCfMDMvnL+xgK8pdEJ66eZeUn9ssolVsYZ+bkx9x2s/mLylmQv6g+pRbs76b12Vr/9SdQQGtay32tpB8uPak0v0o0VrVaiyWt5tT1fSL/NHF0Yqr7oVqruWE0lq1PEaKEDhEooKnft2J8Lib9qOmSefl/WUrzR+lFrdp0erb/Y7IcWbqPYVgaOgwatm41r5a/N0epP4ECmtBaSePWpB89fRA2rDPyJ8tXhdeyW6V79U8qjWaUt19i98I4VqgSLUlq5e+tPkdxoPozmbEynh6ZZuy4b3QQSD9HujvrpWKFXP5U+XFm+4aCchbkT7bKrv5TkweD8md0RtgS3FzRoCNBe6WZ8Uti+XtFCVXqI+WfE3XNzJ9Umr+W+VEXT488od0f89Un3ef4Oun10sWzXml+SUk5S/KnNi/VP5O/+Ls+xCJzOhhqXYPWjRK2c+I+ifnHk6dH4ipMn+PrpNdIV5NnvdL8Eqk+6/L/smT2d0b+KONsQzlmaBVSlfxOY5h5o0QGUcKy/FvrrT4+AgW0NZwTw28F6cdMl6zN76WPDJkK8g9Js1o7aX+n8q813pgxRp7D8K3IsqCBKMHpqlXuOUpQBArYgbF9raSfIT1KVUqp4vzlfQwF+eUTsJi/wjoKObOzf3uY8/njWEH+s8yLrmbkP3qU0JX+60+ggNZ0f+PrpC9PTyrNryqM91dsBaW2P5d/dT3z+acGp38dDHWiq0X5k0o7XIgSOkSggKZ+34rxuZj0o6a73+0/J5/fS0/mV1obFZXp53dbHP5Nt2pBfqXjakv7K+3m1q2anD8zNlHQMifKHw6OdPx32t86+XtzlPoTKKAJ3eP4Oul10iVT+aXyJa1bHa/+ReVndiLRKmfly8/Ofp+exJCOEmSl+0uUkHeg+rOOApowdtw3+lqQfvR093PytCeVI50lbdssNEffcoJX5Wbqt13/Lbny0/VJTy3M9mFMBQCj8ouaC6F8W0J6q1rruP65zQ6zIoxXYD7aS+5vaSu4Nuqayp9Umr+WBVHmjuhRQAu6z/F10vdLF8+SC8pP6mm/cgrCiony53cwTNSz5O+16TXxifskdjmeixEooAXT5/g66Xukq8mzZGn5kj72a9q8WKGwFUn1gnzLSfeF+PWZ35ewnRNHCUm91cdHoIC2hnNi+K0g/RrpkrXle+WM7L1fc001EctaEamPQY8lNqd+scIpowS1OKrbRs9RgiJQwA6MHfcl/YrpUervxbXlS2fbXfcrqTS0WNBqRsQOBvWdxKCD9+lc/sNHCV3pv/4ECmhN9zdeTnq7dEmN8sU+hh33K1WfTVvBTP45zwiYjGFqRC2Vy++8lc07RP0JFNDU71sxPueSfpX0sV+WdeWryR7ytvsl6SS/FzFkJz9GcUMn9V+cvzdHqT+BAprQPY6Xk75P+qBW+ZId90sS5s/WP5E/+6bFrWz+0RLjZ163qM9G+XtzoPoTKKAJ0+N4Oem7pUtnyUXlJ0pXqq/9kkzFOkXlr25lxXjBf3blRvWR8hdEXVP5k0rz17IgytwRgQJa0H2Ol5O+X7p4llxQflJP+5VTEivUugadKEdrlb27snp9MhU5a59EhSizIQIFtGD6HC8nfY90NXmWLC1f0sd+TZsXKzSKEpQditDqFyyEub37Kje8/j1xlJDUW318BApoazgnht8K0q+RLllbvlfOyN77NddUE9EuSrDZ1DdWUEqZ4dmVQoGbxAoNWv2uWuWeowTFsx6qk742+VuSrsXY8d3ogJB+gfQw0XtxbfnubJuKFfbarySt5amDUv7CVrOonGR+P1aw3L6nznJ6XM6q093R+wZK5eu/UShWhB6Fyoxg73p1RPc3Xk56u3RJjfLFVnDH/UrVZ9NWMJNfOj7l5U+e0/RogKIEUUJAalNaNisECmjq960Yn3NJv0r62C/LuvLVZCvYdr8kZ8o/u7kqiBSIEvpEoIAmdI/j5aTvkz6oVb5kx/2ShPmz9U/kz75pr6hCa6WU8S994yyug6FBfY7iQPUnUEATxo7vkk76V/osuaj8ROlK9bVfkqlYp6j8vaKEVLqRpilkYoW1UddUfqE+ZflrWRBl7ohAAS3oPsfLSd8vXTxLLig/qaf9yimJFWpdg06UU/L3ypQz59mVq+rZPH9FFaLMhggU0ILpc7yc9D3S1eRZsrR8SR/7NW1erNAoSlBKDeso5POX1GdJB8OJo4Sk3urjI1BAW8M5MfxWkH6NdMna8r1yRvber7mmmoh2UYLNpn6xwuooQQ3lZCY/ar+jYUH5i+pTln9TPUcJikABOzCdjZeT3jA9Sv29uLZ86Wy7634llYYWC1rNonKS+XXuMVGrWrXJeyXccyVm3lvZeSs7qf/6EyigNd3feDnp7dIlNcoXW8Ed9ytVn02vlTP5q/QN1IpaZt5dWX0eQ28OUX8CBTT1+1aMz7mkXyV97JdlXflqshVsu18S8guMjRyG/4bFSB0Mh2hlM45SfwIFNKF7HC8nfZ/0Qa3yJTvulyTMn61/In/2TXu1+tXrk+lpaPnsyq0dqP4ECmjC9DheTvpu6dJZclH5idKV6mu/JFOxTlH5e0UJVeoT5x/ChWQp2nUzxPUpUpq/lgVR5o4IFNCC7nO8nPT90sWz5ILyk3rar5ySWKHWNehEOSV/r02viYfypXDB5mk3ElTRLsdzMQIFtGD6HC8nfY90NXmWLC1f0sd+TZsXKzSKEpRS9ddRWCIoPzP50buvsuAhVL21yr3Vx0eggLaGc2L4rSD9GumSteV75YzsvV9zTTUR7aIEm039YoUeooRIro9BzXu0RFetcs9RgiJQwA5MZ+PlpDdMj1J/L64tXzrb7rpfSaWhRWkrIpU/N7qyiZutozBpZvlzbrAsfsr1HjqPEhSBAtrT/Y2Xk94uXVKjfLEV3HG/UvUpbfVr5a/SN1AraqlYvhmLMsx9tMQu+o8SFIECGvt9K8bnXNKvkj72y7KufDXZCrbdLwn5m+QXByZ6ixUOESUoAgU0onscLyd9n/RBrfIlO+6XJMw/1W4Vlb9XK95bfayiR0vs4ihRglJKT66gifm0Hh/P5EfQZbCvGjP8P1XCwbkvobHn1lnNxe+96UNRWM450hNHo8t6TqdL52X/ezFVzuhoFDTLJfUsTc9/2X9f9+h3V85UCek8xmit7QV0XJ3oZDNR/e8ERndSUnUaUR1VL8wQVS94SWX3zkvXri9h9KoXDhiV26sGp9/kHqVfTf25x5kbtRf0KGymsz6ufekF474qcQyXlEN6f+mSheUn9bRfOfPLqXdKmbiW7ebcteCaO5FfazWMRkwW1qyDYTr47OavoAgUNmfM79+FmTbjxCdO93RRn8Xpxmg1/lIMP6w9PkOZ3/J336/5XL/CnNPF6pPJdOvr1lHwj2fw94rS/V+18v4K29Uzyp/j1X+IGPqa/Jje3z6aDwIFtDWcW8NvBen5dKX6qg/p+fRS41gkFryop6+NxXLm1FOrb6wwc39d/loWRAmz8nvH2csv9jU0ixXWfn42RqCA5oxR/d3fT/qR0pM9t73VU5ZufeRYoda1tVROMr9WZftr8xfUUrZVlGAJWees/LhR3NB5lKAIFLbTrNNq2We35buiQjYcJ/5+lwvL+V4iFW5Xh5uetV39+2HD43CU/Vr29/L+1+t+KVGmVZD+XrVazVEfgLdfYvnj/dJDX4IY0ySPz8x6/7by/XMV5V9ntK38Ik46+mCUban07/s95s3akzQCBTT1+1aE56A66UptW/4onf1iv+T9Sppu1aTy56mfv3R/5x2H5fVZl3++yQ4G/4fFSo//XggU0MRwLTL6VgzfgYrpIxuUH36r2S/2a3K/xlni+iQyDenJVkQqf6tWVqqP/Hc06c6RSvVZmn+xyT6GZQMTBcd/bwQKaMKcdBya/XLp7FcmXeK1wfF20+TySxuUOvUU5isYVTZvQ6pPQdSVzy+8ML9xr/toCSlK6KD7IIFAAS3oVuPQB0hnv46VXmO/crbOnzJRgLRfyYxCfarMbVzQZyDmF/5eCy7XXbiw5u5KKUrooPsg4VTrAO5Oa+2OplbCYmnBW8briC1YaWvZ4lxt3jVamXFIcsdk9K2I0r/bmp0/ruGc/PYX+675+RNHg/1iv4L9UtN+X3+l9XfNwF9K+j3DmcIrZMm9kpl3ffd7+Na6/frulHAcdFSeMd9TQN3qLXuXrd7wh9Qz3zV/W1JQEJRcuq3hqKY/Tto/tW6JQKGm3eemAgCug0ABAIDjSXSmHhmBAgAAEDGZEQAAiAgUAACA6L+9K4Cv9/v9fr+VUvf7/X6/71yb7ZXu7/v9fj6fcfr9fk+mn5Ld04vs7/ydvexnw32J1DU+FaX7e9kPRn0GHbja32XB/j4eDz7DSqnb7bZ3LRqZv7MX/Gy8Xq94Zx+Px9712sqy/b3gB2Mj9Cjsz86Pvd1u7hLq8/ksW+fgEJbtr72SeDweV+huSbrUjhft7AU/G//+/VPel+j9fv/9/f39/amTdi0s298LfjC2snekcnUuUvYTbcrr9dqpUhtavL+32+2sxyTv8XjYfbfO3aOwbGev9tmw+xscnOQ36xwW7+/VPhjbOeGn6ljsRznoQ7M9ZqdsEhbv71lPgpP8hvOsnwpn2c5e7bMhBdZnvcBYvL9X+2Bsh7sedvb5fFTU0Wo70+xLJ1Nlf92Epit4v9/2u5oZcD2N9Tt7hc+GDaeu051eZX+v8MHYDnMUunCd77y1eH+DRbJfr9fVDh0k1/lsJNs8l3i+vV65v9f5YGyHHgUcgH+muN1ubiT7379/p5y6hfn4bCil3u+3ne53hW4nNW9/+WDUtN+oB4yRR9Fs+lmHG6V0aX+TkxhOPHsr6cQzV2Lzd5bPhmssL/XZmNxfPhgVcbx2lm8429dna3X3Nzk18qwIFIpc4bPhX0+f76IiVmV/r/DBqI6hBxyYHWtkmhJip/9s3O93u5CAbfNOP+5ea39P/8HYApMZu/B+v0//PfddbX+Buu73u71LyJx0WbbA1fa3N/Qo7Mz2gwXhrZ1rE9xTfg7L9vd+v2ut44sA97SIqnXEkVzws2EXM7XjMnvXpYVl+3vBD8aG9hrzgMXKjCa1v6+B/VUasb7UZ5g5ChafjRPvWtLM/eWDsR2O1/7cAqX2U+4urPeu11bm7K9N8b/kLsW+y01rus6kJAIFi89G/trvfBcYM/eXD8Z2TtsaHUvc6753jbY1ub/xd/5qj8uLEShYF/9sJHc22XCew/z9vfgHY1OnfUTh4bhHrd/v9yuMny3b36sdJczHZwNJfDDWI1AAAAAi7noAAAAiAgUAACAiUAAAACICBQAAICJQAAAAIgIFAAAgIlAAAAAiAgUAACAiUAAAACICBQAAICJQAAAAIgIFAAAgIlAAAAAiAgUAACAiUAAAACICBQAAICJQAAAAov/2rgAAAFBKqff7/e/fP2OMn/J+v+3P9/v9fr+3rxWBAgAAXXAxgRqCBv/Vv7+/2+3m52mDoQcAALrw9/f3eDzszzZKuN1ur9fLGGPTP59P+04F7XdxAACAvWitX6/X/X5/Pp9/f39KqWAYwkYPjRtuehQAAKjj/X7f73ettW3slVLP59P+MMlmsx0GdnzB9S5Yri+h8egDPQoAAFRwv98/n0/ypTlNrdb68Xi4cOHz+djehSCPUipO3xQ9CgAArPV8Pm2UYKcUuFkFRfw+A2OMFA00nqZAoAAAwFp2SoF/rf98Pm+328y3++MOSe/323YnzC+zFgIFAABW8Zc68NNnzk6wJWQigOfz6W6CaH97JOsoAABQQdzSzx8jsDMS4nR/NQU3g6ExAgUAAFZZeZUvjTu42ZF7hQgWQw8AAKyycnZhctzBRQmv12vHKEHRowAAwEo2UIjvjZzZ0xCPO7h7KHpYwoAeBQAA6ggu/YOHNSTZYCLok3D3UNSr2nIsuAQAwFpu1uHtdrPhgh8lZJpafylGx94JKWHBJQAADuZ+v9sOgM/n8+/fPxslzFlzaZfnPBWhRwEAgGre77d94oNt/m3fQKapnf8wiL0QKAAAsJUgUHBjCrssnbQMQw8AALRg+xher9fj8fh8Pp13JDj0KAAAsBXXo2BnO7p5iHaZhEM0wfQoAACwlWAlJTdvsfMJjD4CBQAAtmIfGK2iGyCTd0X2iUABAACICBQAAGgtuSBjnwgUAADYXDDWcIhBB4tAAQCAzdlAwd4S+X6/P5/PnHUbe8DtkQAAtOCeB6EOteASgQIAAO0caHaCRaAAAABEzFEAAAAiAgUAACAiUAAAACICBQAAICJQAAAAIgIFAAAgIlAAAACi/wPDq6Fd7T4kggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nPOI = 1\n",
    "fchi = R.TF1(\"f\", \"1*ROOT::Math::chisquared_pdf(2*x,{0},0)\".format(nPOI), 0, 20)\n",
    "fchi.SetLineColor(R.kBlack)\n",
    "fchi.SetLineStyle(7)\n",
    "plot.AddTF1(fchi, \"#chi^{{2}}(2x,{0})\".format(nPOI))\n",
    "plot.Draw()\n",
    "c.Draw()\n",
    "\n",
    "c.SaveAs(\"test_p0_1.png\")\n",
    "c.SaveAs(\"test_p0_1.root\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acceptable-gossip",
   "metadata": {},
   "source": [
    "Draw all canvases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "passing-boating",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-07T01:25:56.398175Z",
     "iopub.status.busy": "2021-04-07T01:25:56.397827Z",
     "iopub.status.idle": "2021-04-07T01:25:56.525011Z",
     "shell.execute_reply": "2021-04-07T01:25:56.525257Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAHYCAIAAAApvgy/AAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dYZarqBaAUXir55XckSU1siQj4/0gEgQOiiKifnvd1V1FCKKVyBEQtTFGAQAApPxv7woAAIB+ESgAAAARgQIAABARKAAAABGBAgAAEBEoAAAAEYECAAAQESgAAAARgQIAABARKAAAABGBAgAAEBEoAAAAEYECAAAQESgAAAARgQIAABARKAAAABGBAgAAEBEoAAAAEYECAAAQESgAAAARgQIAABARKAAAANF/e1fgVLTWe1cBAHAVxpgGWyFQqMz92bTWC/6EC97VbEO8i3dVfJfN33MNeddeG+JdM9/V7NKUoQcAACAiUKhMD/yfGZIAMIc9XTyfz70rgl7ErUn7ZoWhh8pWDj0AAOBkRugYegAAAPsjUNiV1uE/AAf3fr+VUs/n0/6wppz1hQDr0T1eU9lwQyos0K1ud+kfYze+8x2NNXvU89GIe4Mfj8fMOQfv9/vfv39Kqdvt9vl8XPrtdsuECz0fjcaudiia7S89Cnsz5vcPwJHd73f38+12sz/8/f3FzbydjCY1/36UYH/1SwYaI1AAgDpsA3+73Ywx7/fbXe0tGD54vV7GGFdCEDoALREoAEBNtgPADjfYxr70dsfb7ea6EF6vl/2BeyaxFwIFAD2KZ/r6U34zr273b9Lj8bA/fD6fv7+/YFGE9/t9H9iU5/Npfw2CAH+ggUEH7I51FPrCPAXnUpOSJnE0xjo9Grbhfz6f/kjB39+fGm6CiCcfrN8onw2HQ7ERehS64F+vcJskoMbTfDv5N8f9frezE16vlz+f0b50G7j89teg28Cf08DtkdjdtW4m2dqS2yONcWGBUVoppZVRinsgcHLnu5PN3dzo3w/p7pYMdtamv14vP0RwJfiFSCUA3B55OZwEgENzTb6bneDaeDd3YT5XyOISgFoIFACgDneHgm/+gkuOPzZhf+WWB+zobL1/+1oz9GDMN4WhB1zB+YYeHDurwN3jkMmTv6Ph/X7bVZwr1w9n0exLdNrv6i4IFICZThwoAG0wRwEAAOyPQAEAAIgIFAAAgIhAAQAAiFjCGRtirRgAODp6FAAAgIhAAQAAiAgUAKAOu+hy+wdDP59Pf71noC4CBQCoqcrDo4vwhElsikABAACIuOtBKaXsaurSmup2xXU/JwDk2fNG5nEPyfxKqfgt/oMh5hTrzlfBM6zjRGAWA2OUUrfbLflS8Bg3pdTr9cqUU7RV871v0Hzfp5RR6vfr8fExg+SUnwr3gQ/OG5mThjEmfoS0fzpyT6SMs7k8bnNBNZKFVN1j7KnZl4ihh1x8/Xw+7XDj6/V6vV72i/rv379mdQNwRME0hcxJ4/1+//39xW+Pz0txNuncFV/e0CeKVdrEIx16PB7+1ynZo2Bf8q8G7FukqLzseNKjgAs75aciPp+4k0zmpBGcZ9xbbIrrDPCPWJAS9CjE/QeZsxyOq9mX6Lo9Cu/3Oz852Ybet9vND9ttYhzaA4DjruDjH5Iej4c7z0hviYceJHFRQTpQ5NKBgo2VpK+fP4HIkb6BACrTurt/KyQvS/ITDNecZGwfg92oG25g3AHLXDdQmGS/Y8TgAA7HxQTuZop44gIwE4FCMft9o0cB2JYx3f1bIdlOu+sQ/1q/yn2MfslBvwJQikChMj2fUiy4CpySa6fjH5I+n098Y8LK7kx/9GF9adhaQdvxa0MaYcGlhaRvnZl/2cHC7MBJfT6f4DwuXdC/Xi9782RwC2UwjXqB5/PpymTcoX8FbcegWaxAj4JI+mq1X8gdwLEEZ49MG3C/3/0bIN3b1w9uxrdrAcvQozDBrpYap9OPByDgBwTzl3C+3+/GGCm/fTWzISVPmbrdbkzKxnoECiLbcRd8A5nDCGCO+U95WJZ/klsqhnEHrESgILJfWjvJKJicnFn5RBo0WjD+BADL+Ccixh0OreWkRQmBQs7j8fj7+/v3758dRHQReuaLR0AAoB/+mo84okyb0iyGIFDIeT6fNjjwJyTHM48AoCtcsaAizedpDtuFMDmIqHXJ8bTBoDHD/78pWpnvr8fnAl4+ZgiUfVkARJp9ifiu1kSgECBQgIRAAVip2ZeIdRQAAICIOQqVcdcDAKAW7no4IQICAEAtPdz1wNADAGDC+gcR3e/36o8ysgVy/+fWCBQAoI73+22bropLuK5vCxe00PYtjTfaAIHFMgw9AEDvGj+Lzm7O32hmOdqZ7LI0KwvBLggUAKAO1xDGD5Ozz3xSwrqu7tXFT3zIlCAFGe4tUq18yVqVVrtKIZNP25IKzIQp64//yRnUU3Y8lTLfiSrm+z6ljFK/X4+Pjxkk5/tUxM9esunxWq6Px8N/Y3xafr1eyTfa9FhcgtvEzPTMS3aj8Xc53t/8Rl1+V0Ly+dqZIxzkd50c7l3JZXOlgzn5lv6pVl+is31X9xV/4OJvl5+bQAGXdb5PRdA5b1svvx3yW9ZMg+1/Zfy33G63ZAMW5Alau2QJ0vMkM2/JV6xoo67myTpIsUJm7Xz3FimDfXtQqznHvwdSDVvWs6PDcQJlfzYCBVzYKT8VrjFzLXqmdfR/ja+J40hC2micIShBqkPRRpPVzpQZbzRIcXFVpp4+96o7tq4OQUzmCnSbCN4S73jm79WzZvVkjgIAbMXODwj6Ev7+/tR4HsPn87Gj48/n0yxdi2V+CX6G9/u95jnUi6sdT494Pp+Zmrgj6Q7a6/XyH9d3v9+D/bLHOV959xZbn8aTRg+jTTxyEWXHkx4FXNjkp2KnM2LO5E7FV6iZ0jKjAP4khsmtJ+9H8EtIXiXntxtvNEhZsNEgxf6cn5TgizcRFxKMLwSHOrlR6S0za7WvZvWkRwEANidNC7Dz7Z/Pp38tay+FZ17oP59Pe00flODf1BDwlzewTf7kxff6jTbgOhhcx0N+v97vd/CW9/tNp0JCm3jkIsqOJz0KuLDJT8VOZ8ScyZ2SehRmXjf7V7fxbLuZJfjX+jYxuJSPx/IXzFEo3WickrwJ4na7SccqPpKu2jYxnpEwOUchnvfAHIUkVmYE0KM2Z8Ai8yvvrqptw+NfpPpLFj6fT/uz7TywV7RFR8mtBelKWDDhoPQtbqP2qn3ZRt1UA/fef//+fT4f6YLeHUk/f34Tpd0kijkKkr2/eqdSdpzpUcCFnfJTEdwMGae4nydvj0zO7U/eHimVEN9NoKJFESbfIt0eKZUQX6yraDZGvubBLITksU1uNJNBulEis9xkP0spSDUMDua2dWizmYso+7MRKODCzvqpiE/i+QWXks2bn8FvzJKtV7KEZBe9KzloIP0M0kbj73LpRhssuBRMBAnGJoIM8V7cbre4zJ6pVl8ibbocCzworUuOp51PZMzw/2+KVub76/G5OVN8zBAo+7IcSnKN4fzCw3bEQcowuWjxZAlqmH7oXpq5EPL6jaqpoY05G4o3urjacYbSCvSj2ZfotN/VXRAoBAgUIDlxoAC00exLxGRGAAAgIlAAAAAiAgUAACAiUAAAACICBQBoxy1+vHdFgLmYeFyTv4J6IHGcuesBF3bNux6CU8QFjwBKZZoV1eojRI9CZdKCFXvXC8DObC+CXevQrvxzuBv30V5mHaRmdSBQAIB27MI+8dCDfe6DHzrYQQqbWHGo4vl8BqW5Z08ASVfs/dsOCy4FGHqA5JpDD879fv98Po/Hw7XZ9svyer1srJBsudcfsefzaR+V5BcVVwaHwMqMh0SgECBQgOTKgYJtmJX3vbApt9vNrk/smvPb7WYfKWnzuwylXLHu6Yj+wX+/3/ZJjJf9ixwUgcIhESgECBQguWygYL8UQZMfdCfEkYT/VfLf6D+wIEiJN+oLDr7NQKfCsbCEMwAcjJ1V4Dfk7/fbb3pde+zncRlcG+/6D1ye4LmI/wbujfZXqb/hNpBqbp+aaLsxgFDJoyYxoex48phpXNj5PhXxw5rN8BWwPwcPd1bDs4xt+x08ezp4onTwjGZX1OPxSG5XEj/rOahq8knW6FOzLxE9CgBQwfP5tGfV2+1mOwlsD4FryKXLfdt/4A8ZxI+oDvoYXCfE39+f+9mPGBZbNgcC50agAAAVuAb7+Xy+32/XuvvTCIILtTmt8vP5tEMMwbQGFxbYrTwej5WrMtgohEABsf/2rgAAJOQXMdzr1QzXTtupiDZumLzKz09CdLMa4/sd7vf77XZzdzEwCRHboUehMi3Yu14A2rFNuL25cXEhWmvXW5C80HdRgqIn4LykNqVls0KPQmXzrz8AnJVtwuc03v7AhB9VuGbA3TMZCLoQ/v37x8nnlDJ/1maxAj0KAHoUDOd38uoc7pkOC97rl2DZGywdl+huZXR3QPiv2ivOovGIeE4l8DX39gjMUHY8uT0SF3biT4X9zNtbH+ewIYWff/KM7X61dzMGv7qJEf4tl/62kgffLwGH0OxLRI8CAFRjL+Ifj4c/gSDPXsTPz+/6CdwcCBcZuCWYSsWLPgHORVdR3QhLOAdYwhmSsy7hbD/zxhj7yEfb7gbzD6R31V1BWWs9v0CWcD4ilnAGgIOx0YCdNGCXVbbpk61v9RWU7TzKmd0DbtIlUQKSzhnU74UehQA9CpCcskfBdSf4v6p5V+rBc6HWsE+DpDvh9OhRAIAjeb/ft9vNX2HJ9hO4FZ3z7BurLIdwv9+NMUWt/sxK4ppOGNTviB6FAD0KkJyyRwFoiR4FAACwPwIFAAAgYglnAPvgGSjAIRAoANgBExSAlZqF2gQKlUl/OU6LAIBSPXS8EShURkAAAKgl06bw9EgAALA/AgUAACAiUAAAACICBQAAICJQAAAAIgIFAAAgIlAAAAAiAgUAACAiUAAAACICBQAAICJQ6JBWSvewvjcAAAQKAABARKAAAABEBAoAAEBEoAAAAET/7V0BxIxSSn4EOQAA7dCjAAAARAQKAABAxNBDZdL6B4axBABAoR7W1CFQqIyAAABQS6ZNaRZDMPQAAABEBAoAAEDE0EOnXJ8SYxkAgB0RKHTHKK2U2n/6CgAADD30r4cprwCAy6JHoSfGKKWGwID4AACwP3oUOsXUBABADwgUAACAiKGHfn27FOhaAADshx4FAAAgIlAAAAAiAgUAACAiUAAAACICBQAAICJQAAAAIm6P7JpW35UaWX8JALALehQAAICIQKFTwaOgeDQUAGAXBAoAAEDEHIXu/GYjaGV4iCQAYFf0KAAAABGBwmEwTQEA0B6BAgAAEBEo9I4VFAAAO2Iy4yEY5U9yBACgFXoUlng+n8/nc+9aAACwOU3P9gJa69vt9n6/4/SC42knJxoz/L88Q/fcBEw+ZgBQV1mLswJDD8Xu93vT7Wn9/SAMNz1oRiIAAK0QKMz1fD7f7/fn89m7IgAAtEOgMNcOUULcacBSCgCAtpjMONf7/TbGGGMej8eO1dDf/2itlf8PAIAtECgAAAARQw+VlS60vGJK4ne+K90JAHB0PS/ST6BQWfHtkQCAy1two2Oz2IKhh4PhpkgAQEsECgfWc1cVAOAcCBQAAICIOQqVSVf5FRfaNMbQlwAAV9DD2Z5AobI2K2/z6AQAuILM2b5ZDEGgcEDhh4OgAQCwFQKF04ljTLofAABL8Zjpmio/ZnpeCaOHSSZ7ovb7E/OYaQDYCI+Zxgruo9PBLBgAwKERKFTW4K6HeJtKKa25ageAs+GuhxOitQYA1NLDXQ8suHQePQSeAICTIVAAAAAiAoUTYLADALAV5igcnlGMOAAAtkKPAgAAENGjUFnT2yONUYm1ErRN4v4LADi6HmapEyhU1k/zrLXupSoAgEV6uD2SQOHwgk+R/8nRTHQEAKzDHIXzITYAAFRDoHBKxAoAgDoYejin77BWpREsnlwNAJdFoHB+K5/13MGUWwDAbggUKtvj6ZFzrXl4OU+uBoD2uD3yhHoICDK+90xWrWTyY9z3YQCAY+jh9kgmM57fpk12B8EuAGBD9Cic2hCKmo1jTz/kJXQAgDMhUDin1H0Kv/sguIsBADATQw+Xo7//4cIfADCNHoWzifsG5LEAbYxhpAAAkEGgcGlaH+BxEAyUAMCOGHq4nGM1snR4AMC+6FGorOcFlzzmWHMUWO4JwDWx4NIJdRYQiFw99XAfxFFqDgDX0cOCSwQK+HLzFcy4s8EopQ82XjHCwpEAsAaBAnz6ZA1oB512AHBsBApXZ4yJ+69+AxE9tbSL68LCkQCwGIHCVdgGMtlhkIwV5pS2piZRHSYyaK2MoZkHgNYIFKBUNLfR/n+L6Y2TLX0mSmiwdQBAgEDh/EZNrI5SZN9nUmdKq1GlOR0MtaIE+iQAoBSBAvbgNdcm+H+UsWKUAAAoxcqMGDHbL91oZq/1RJQAALujRwFJv7se9PC/ylMWplZbrNi6EyUAwGIECpUdYwnnoJLz6qa1Tu+FDhdoGv1/qiLJfPQBAIBiCedT6isgqC0RK6z4EH+jgdkrJ3bwfQGApljCGc3NbubdugthU65TN06UL2lU2mcw/56Frb87k+WfOlYEcDkECpjmxwrfH7RW67pPFkQJM3MufnV9+QBwPgQKSPBaZTP637izK7nQgnst+avxU+bOjSi+RpdGLmqto8DICIDr4PZI7GReo32UOySZfQngrOhRQJnvXZPDr/FIhL2PIfFMKfuqP8NhKlYgSgCA3REoYEmnefRsiO/P0qwFreyDp6JX5ViBKAEAesDQAxqxMwTCRCFEIUoAgE7Qo3Btq1u5YCRiajvGeE3r1usoECUAwHoECqhgPOKgVaZRt2MNo4BhltJ7FogSAKAKhh7QnDFK643WUcjnr3UHI1ECgOugRwGtfZ8vtWIdhbi9zywUGZQzPORqOn9GPgpZWTgAdIVAAdUl5i3o4f7JX6ufjxW8xtbec7myvW0zEsGaSwDOh0ChsmM8PXJLbt6iio5GeP+k0LTGN0cEhfut/synP7j85U+lKCi/3gOzAEApnh55StcJCOrIX4a7uMBTcb5ClT6AredDALgynh6Jy/E/2d+RCClramyibpSw3tZRCADsjkABLeRWcsy0puPG9lhRAgCcA4ECFml2sTzECs2ihMlbKlaWDwDHwjoKaMoY4w+5aa21UnrGYFvRPICKUULd8gHgcOhRQKHVLaG/jkI4EqHkKQvpctJNuzGJNl/PnmmamkOZyJCsDwCcDD0KaCpcR2FG/olyZr9tl7UgAeDoCBTQRHKewfcKf7Sywnckwo5KrOzhTy5uEFZqRflL8wPAgTD0gCa+z3eImlNjhx+kQQH9CyOWrdUo91tMjFzMxjoKAM6NHgW0oHW6zdZa2fZ62+2m0rlDEgDmoEcBLSSf7+Cu6cePnR49KmJ8N8TQ8ZAKOGJha93q+RF1o4SZt2sCwEboUUBb0roI4/4GMxi/uazjocFToOanL94KAOyLHgU0J89XsEMF80sSW+WgffUfITV57+OK7W404lD0BCwAqIseBbQmzVdQyo4t6DDBCx209b0tongVpqIHMPS2jgJRAoBdECigqcl1FIK5jTZ/chwiORIxaq2Tz4mYvX5DV3dIMjsSwF4IFNDEjHkJv7xKnscwtZFk+VI5R1lHgSgBwI6Yo4AmptZRmJnfrv3st++j51ar1OyE5HZ/L6bT52uwjgJRAoAdESigBf/5DkG6+o41hOlSfmOMsY1w1H7GMyGlckorn3kpM1+h/ayCBVtcE4WsD7MA9I+hB7SQnB/wuwVBhY1qPn+G9v8J8xJqtd+9rbbUOC5hciVwEfQoVCY9LnnugwtP7/fcyNR8AumeyTh/mCPJ62IIypnXyGU6IxLpQy3cCIixMVBbQa1ajowQOgDVSW1KSwQKlREQTCtdR0GeZ+BafaO+7VXuS5Wdr5C0IEroSoORET7vwKYyp6xmMQRDD2itaB2FTH6hAyK17oJbfqFwHYXikQUzPK9yg/az9JzQ28gIgIMiUEBTy9ZRSObPtIK5dRfWraMw3cpuM4mx1h2bRAkAShEooIla6yhk2+C56y4sXUdh7nh/7ViBKAHAjggU0MR3fkAiPf2MaTl/pnlPjkSMnx/5G4VQ04FEUHxJetUG+RBRQpeTNABUQKCAFjLzDIKxhsn8ueZdyG+Gx1GOXiqs/6o+hob2ihLoqwDOirse0MJonoHX828T7G2L/p2EmfwTm0nml25zHNrzYQGndOnidrVOrAXp3Ysx3siomjN9QyjtflVaGb9AP8Yywf+D+iwyp/KniRLW/KWAE6NHAW2Z9PyDZL9CJv/vjcL8gOT8huCeiKioqadM5TY8lX+BqX6J9BGbV585fR4z85yjNWXoBJDQo4DmxusZGL8fQV5zKb+Owqz8yVihaHaknx7V3r4QbNf+5spZ0Bpp19OhdbhBN9VCWBihyjoKcvEqrM/xLf4zASdGjwJam5h/UJIuXisXlBP0MWhvymP5vISSetZUHiXMd5o+g5mutr/AJAIFNCWti1CariZn7c0uRw3BQtQJUTASMb1dZec0aKO00t6/oBThpelr3GCqZr17Hy7VahIlADECBTQhzxtYkJ7ZSGn5U7Ue3VCZaEWm6ulHCZMby7xyiDskT4DjA8QIFNCEya2LUJqe3kSV8sORCL/4wvKThdjCjbzGs31p+GfnBhTNRtwrSmBcHzgrAgW0UHFeQq6JXV3+UFJ6Deigg2GifEH8TO28ZP5Mq79XlMC1OHBWBApoocq8hLA1isfyy8sf3WEYzSHI3Ew5WtsxOy8hMRch2y+SSMwfhylECQDWIFBAW0vnJcxtjQrLEddv8N+aGZBYU89sOVJ+aZ6EJHePRonSkQ4Ap8E6CmjOTK5zMJXufk42U4Xlf2MFY9TQHI76A7z8w3IGbqFE46/RZILyR/GIW0ghrFSCEdaRjEqdzJ+frzA/VsjkP2WUwHwLwEeggNa0/rapqTv6dkhXKt3WyvmNMYnGRNsooWi7XjFm9Gu4pnWG1t/8cfqc+Qpz7sbYLhqIt146qrLm7VKZRVEUcHoMPaCp5fMStkn/qVGOneUYzHac2K4gGBOZbs5VHJ2ksvW0jsLKxniLtpyRFCBGoIAmqsxL2CB9ZJjbOL8ck1imye1xtLaj8Ext490saby7I23+7EDE6C3epmdFCf4bjUn8U01GFuL6l75d1asnUQIQI1BAE6bmOgrbpYvrFkyUk76d0r+h0i8/yiaTtltvHYVafQ979dXTBwBsjUABLei66yhslq7kdQvmlCOECzbnL1ZIl5N+l7hdqXXcK0rYq7UmSgC2RqCAFrqal+Cnp6/v185XkNZrmip/lGnGduc5cZQAoAEChSXe7/feVTisDuYljIbzpbGAGuUH4YIxSrmnU47fOqt8W2FhHYUGIxFF5Zfaa+QCwCQChWLv9/vfv3971+LI+puXkG6mKpVvhDmP37UdVXp9yYKRhSC/Hj2CMhkGVexjqBUldNUnQdQC+AgUCrzf7+fzSZSwUofzEvabD/H7N6ueKQvybzoSUarDKKGr+gC7Sy1gB8H9fv98Pvbn5HHTyQUBJfbMPqyUc8q/g+tid4clHI/vIF25gz+scWTMVtsdVjtIN+rGmFw5KvqgBB+deLvRZ2tNlBBvLV+dSZMFzn+7C5MWf4+CEk78rcRplLU4azZ04kDh+Xy6/8be77ebaiDlkYr9+/sjUJjjFyio754mxtd3Tf9WL2rLt9uu9oLM8VyFEePaK1fOnJY52O44w8q+hG4DBVvPld+jlZUB2msWKJx5Cee/v7/b7ZZ8ye8bsDlfr9f9fm9Uswsyq5/vsE160FJ/5zZutt1h7qRLFzsYvqcAqfxvH0Mkn3/LEYe9xvUZKQC2dto5CplW//l82ijh9Xq9Xq/H46GUYubBpnqbl5CdN9C6/GGyYxRVDEs1xR0PmdZdqs/WUcJerTVRArC1s/UoPJ/P9/vt9xbE/v7+lFKuC+F+v9u3PAfx3Y/P55P+hjW+Z/Px4P/vx/3S089e2my7+fLts6OSV+Zaa78UoxLdENP1kfIHqXH52daYa3rg3M4WKExGCXY6wu128xt+ey/D39+fDRS2reLFDW1Y2LrslB6NBWy7XRsF2LbYaK3i+RDDhEatVRwLjNruX6q83eHHWCZKSJJeKYoSMsMTDUYukptoGeLEFSDAwiGcbejh/X7bVW7sgEIyg4oGJtyvrKTUgrnWOgqZ7WbWejJGhcs7jouwwYb6zWNIlSOYaN3N+MFQ2fy1ooStm8yiQGSLqGWvORzAemfrUZhk+xu2G0fIzGNPuuAVhdbfPvaw5dspXan4hsittjvqzNBD/0HiCdG/cvw7S5Ny2xXzC+npLdRsxRMBU8ORi/D4p/h9P1tUYNPycVylbUdLZ+tRWMzeHzGnR+H5fGbuSDHzXTdKUEqFc+72Sv/poz759KCPwZvsGDY+ZfMSgvSGJ6xaUUJplauMpCzAfA5ICtqOXxvSCIECmrBD7sI4+o7pI+Y7FtBhPcdDEMZ4PQ2WHyusihKUfHxqqxglFJVTZSRlGaIEHBGBwgi3NmwlM47eU7ptbPupTz49rPzvXsrkwgwl0YO0XdmCa/reogQASZcLFKQlmPL3SmAl3dl6CVK6EuY29lbP7/yGqF/By6CD/KV9DEUD6Tu21kQJwNYuFyhY0lyE9T0KWrCy2KPrcLzfpmfuO+iqnpn0zICl7VtY3scwr1E90DU9UQIOR2pTWjYrlwsU7DIJQaBQ8a7I3Wed9K6z8f78PYr91HNe+viTFnYqzOpjUPK6C0m1ooTLx9JAGpMZd2D7DD6fjx8c2OhBWnoBlXU23r/XOgobpf9OIlFW/1qkVh9A6cmqaARkLw1mcQIHcrlAQQ0Bwb9//+wDJN2jH1iTsY3exvuV6qs+FdONEm/BTT4/QrJpK95hlHCp/QUmXW7BJeU9D8J/ENTr9dqxStcRjrsPp8y90n++be3O9amcru0rLjEzGGFc/pZNWa1NlV6j1xpJKUWUgCNq9DTrPtkuhPv9XuuuyLKngytlXGkAABwNSURBVOvvWXz4f5Uq9MU1RUZ99zTcTa33Tf9Wz/zSdfz8hQ7quTB9/NnS36UW0o2kcd0SQ8yhtTLqV0LwQY0/t/lPcjL//LdnChz2a+7Wk/kn926lrcvHBZW1OGs2dOVAobrMNNTEcb5UoCBdx6s5k+83TI9P30b1WM+F6ckPpNDBMM4ybn5TJSRbPjU7UIirvyxQUF7zP2frbrsECjiE/N0NbVrwK85R2NTu01P7pPsbv5daweTcxt7quWC/YnYhhuRL8ZrQeVLQUit/xl7bBdrgrgdcxWgcffA7a++XfoJ1FHLpxihjtBo9CjLIP70Mg/o+qVKLUcVhWmuiBGABAgW0NbRhYWuxU7q9cO6/ntXS3Y5H6fkrFKn/s1aUMLvzAkBrBApoznSx3oCTHGvosJ7V0rOtuxm6OlMv6m8Hg3e4Sq/RS+uzi62jFqIiHAuBAlrrcfy+p/o0SJ95TZ8ZB624fGyHUcKm9eltf4FJBApoqqPx+0Af9Wmwv9PX9GE50iQG27mwas35Wq1maRWk/EQJQIzbI2vi9sjAUdZRcOzdcyr+c+xdzyrp7iOnUi9+f1K/mwiDcrQ4SGPflDieqWyp7Xobz+SPq+1VdlasYIYbIzPhYlB4vMWZ1ZNqmyl/U2sqP6e06gWe8pS4QA+3R15xZcZNEXilfcfLo4Ozd3rwHfzObeyvnlXSw0SllHSNK5WjlLGzOpRS40mgWgcpE1+Ein0JReUs3m69wZYd1K189UNx6GO7tUyb0uwBkgw9oIUOx+mHtlD5/5TqtJ510iNSq5ksx+U0ShulzTgWiPspbIoxLr82SiutldbhdrVWWvvZ5rflv1qN/47BX9b+87ebzD9nc8e9FjDp+2SXl6aqHpC61UNFBApoIdkG7zVOf930sen0qJyYDRfSyzDoYRJD/K4gSqhq+f6uK/8oah3vox8HzEeggLaGtic8y5DeJn2QOcvnyokvw0d55bskEkUHGepcThIl5FWs/6GPA4oQKKA509m6AldLL20tsuWnstsuhmgownYvqOT6VnVanVpRQuaau1Y9d3H0KAd7IVBAa52O318pvah1nCg/lX/obihZhmF141kxSti0Nd2xtT5ElLBXFIUM7nqoLPcMX8Tj38NhIb1NuhJai28WHbx5RvnJcjz2N61c6S6nHuWJyi9ClHAOevadrtfR7NaGDAKFyggI0nTq/v6hbSC9TboSrC1fKmf8VhsPJM962o8VvE3NdJQoQR3kmr5UrYaMKCqJ2yNxGaazcfpLpkepSkln59LyZ5/ljcrdJWFnMBSd/A4UJZxSxePG8e8WgQJaqDi+Tvry9Eim1SwtXywnRWuVmcGgSi6VdokSOugM7gLR1UUQKKAF449zD8Lxb9K3Th+bTp9XvprTWgv5bbAgtDXf+yRK+1cbRAmHbh0ZKUApAgW0NbQx6fFv0rdOH2TO8kXlSxblr/CwSqKEvIr1P/RxQBEmM6I5Y8e5o9MM6W3Ss9f6Cdnyk+9It0bZWOHXxzD8HgQHibskMtuNNmS8sifqmT0kVVpHW75ROjEdY8vm9+hRDvZCjwJaazQeT/qc+QQzrtQnyk/lF1vfTDOVqIk4KJGsdCZKkDZYpY+hlC3flE3ZrOMQUUKtkRFURKCAptqNx5M+fz6BMVqNH500v/yxBem/ao3LH5LNnEkMifKN+e3X6nrWkhiRMTwK6Yc+jz4RKFSmBXvXa2+bjruTPi9dsrZ8qZzxW+dvN5k/f5dEspehVjRTUVj+Kc4MtXaCKCFJalNaNivMUaiMBZfSTGfj9JdMDxO9F9eWX2tWYKae+jvdUcoQTmIorOc+rZTdnSOfNGz1q7RZRz4MG2LBJVxFxfF10penS2qUL7W+0gaL8iuV3G52EkNhPUujnGpqNbN7OHiQg7kIFNCC6Wmc/rrpY78s68pXc0YWSvMnpdpUM4hK078pDNny6/SFrNE2VmCkAKUIFNDWcE5cOy5O+rL0Qa3yJWvzS+RME0s92o1EI7v7RwlWqya3Yv2JEq6DQAHNGTvOTfpO6dn5AaXpidIXzD+o1OpoO41BGpD4ZdNe/nQ5klr13AV9AFiGQAGtVR53J73ifIIF5afKyZUvWd14Btv1hiTSW42GJtLlVLdja32IKGGvKAoZBApoqv64O+l15xOUli+VMzv9V61x+UUmt2uU8gcmgvIn+xhq4Zo+j+PTJwIFNLHpuDvp89Ila8uXyhm/df52pfxz6y+n//oYovxD18K2t6eH2z3F5XOtnSBK6BaBApownY3TXzI9Sv29uLb8WrMCM0UI9Re3O1XPKs+gWmsqhutfxdadKKFbBApooeL4OunL0yU1yp8YWViXXymxTRV3a149XR9DXJfkDIbpepY6cqxAH8BFsDJjZdK1yMVXbDT+OPdwKMLxb9K3Tne/B9aVr5RY/mhSgJQ/yFIeK0zkl45DarvJ0rXWNpTQ3yUfbf7UG5Z9zdvGCowUHEsPTwCgR6EyI9i7Xt0Yzolrx8VJX5YuWVp+SDqpTbWFBfXMCjfibTcbJQx5lLjco40L7DyG+u16q1NExdads1obUpvSslkhUEBzpsfx+0ukG+9BhSZ6buGi8qPUaFv+FuX8VU566VZQLjoR/fwOiBQwKDX0I2i18MGPe10i0geAZQgU0FrlcXfSO5z3UJpfVWg8MyMLxfmNMd/nRYiXbq6DIfWI6yX1bOAQUUIHHe0IESigqXCcm/Sjp0vm5f9lifIXyZS/LF2rb6zgGPefxLvmVpdr+jyOT58IFNBElfF10vtMH4TN5Yz8yfJV4bVv9Sjh+6tK9ItIPQxa/W6TyFQ1fOcpLp9r7QRRQrcIFNCE6WacnvQt0jPzA0rzS4TmSGxdytd7kOqfzuz+k9jE7JZTLv8oKrbuRAnd4vZItKC1MvacOD4ZkH6idBVbkl9I0coYbzqAa17F1qVku9LL2k5UEGMFv76jPEGskJugnqqnt/WcUf9HKnO+6Y3fMllgXIK0iXzhcwQlBG8vOjiThU/mr2vOse0qbKJHAS2YPsfXSa+RroST2vL8gnjewASp/tF2k/n97Sbz+/UxUveCX4wSVnBKbbeoo6F65upRQl0ry9+3E+eIXUj0KKCt4fopPMuQfsx0JVibfwgmop6DiXkMYZ+EmZHfeN0BXj1H+fX4xeFX5Q7NN8UMxeZaA7uCU7i3qbe4OgTHZ7L/QNr+dLQUVWpZVKfGlS8SbFd6e7C/pZtbXL0q4g6SfeuTR48CmjOdja+TviI9SlVKCfMDMvnL+xgK8pdEJ66eZeUn9ssolVsYZ+bkx9x2s/mLylmQv6g+pRbs76b12Vr/9SdQQGtay32tpB8uPak0v0o0VrVaiyWt5tT1fSL/NHF0Yqr7oVqruWE0lq1PEaKEDhEooKnft2J8Lib9qOmSefl/WUrzR+lFrdp0erb/Y7IcWbqPYVgaOgwatm41r5a/N0epP4ECmtBaSePWpB89fRA2rDPyJ8tXhdeyW6V79U8qjWaUt19i98I4VqgSLUlq5e+tPkdxoPozmbEynh6ZZuy4b3QQSD9HujvrpWKFXP5U+XFm+4aCchbkT7bKrv5TkweD8md0RtgS3FzRoCNBe6WZ8Uti+XtFCVXqI+WfE3XNzJ9Umr+W+VEXT488od0f89Un3ef4Oun10sWzXml+SUk5S/KnNi/VP5O/+Ls+xCJzOhhqXYPWjRK2c+I+ifnHk6dH4ipMn+PrpNdIV5NnvdL8Eqk+6/L/smT2d0b+KONsQzlmaBVSlfxOY5h5o0QGUcKy/FvrrT4+AgW0NZwTw28F6cdMl6zN76WPDJkK8g9Js1o7aX+n8q813pgxRp7D8K3IsqCBKMHpqlXuOUpQBArYgbF9raSfIT1KVUqp4vzlfQwF+eUTsJi/wjoKObOzf3uY8/njWEH+s8yLrmbkP3qU0JX+60+ggNZ0f+PrpC9PTyrNryqM91dsBaW2P5d/dT3z+acGp38dDHWiq0X5k0o7XIgSOkSggKZ+34rxuZj0o6a73+0/J5/fS0/mV1obFZXp53dbHP5Nt2pBfqXjakv7K+3m1q2anD8zNlHQMifKHw6OdPx32t86+XtzlPoTKKAJ3eP4Oul10iVT+aXyJa1bHa/+ReVndiLRKmfly8/Ofp+exJCOEmSl+0uUkHeg+rOOApowdtw3+lqQfvR093PytCeVI50lbdssNEffcoJX5Wbqt13/Lbny0/VJTy3M9mFMBQCj8ouaC6F8W0J6q1rruP65zQ6zIoxXYD7aS+5vaSu4Nuqayp9Umr+WBVHmjuhRQAu6z/F10vdLF8+SC8pP6mm/cgrCiony53cwTNSz5O+16TXxifskdjmeixEooAXT5/g66Xukq8mzZGn5kj72a9q8WKGwFUn1gnzLSfeF+PWZ35ewnRNHCUm91cdHoIC2hnNi+K0g/RrpkrXle+WM7L1fc001EctaEamPQY8lNqd+scIpowS1OKrbRs9RgiJQwA6MHfcl/YrpUervxbXlS2fbXfcrqTS0WNBqRsQOBvWdxKCD9+lc/sNHCV3pv/4ECmhN9zdeTnq7dEmN8sU+hh33K1WfTVvBTP45zwiYjGFqRC2Vy++8lc07RP0JFNDU71sxPueSfpX0sV+WdeWryR7ytvsl6SS/FzFkJz9GcUMn9V+cvzdHqT+BAprQPY6Xk75P+qBW+ZId90sS5s/WP5E/+6bFrWz+0RLjZ163qM9G+XtzoPoTKKAJ0+N4Oem7pUtnyUXlJ0pXqq/9kkzFOkXlr25lxXjBf3blRvWR8hdEXVP5k0rz17IgytwRgQJa0H2Ol5O+X7p4llxQflJP+5VTEivUugadKEdrlb27snp9MhU5a59EhSizIQIFtGD6HC8nfY90NXmWLC1f0sd+TZsXKzSKEpQditDqFyyEub37Kje8/j1xlJDUW318BApoazgnht8K0q+RLllbvlfOyN77NddUE9EuSrDZ1DdWUEqZ4dmVQoGbxAoNWv2uWuWeowTFsx6qk742+VuSrsXY8d3ogJB+gfQw0XtxbfnubJuKFfbarySt5amDUv7CVrOonGR+P1aw3L6nznJ6XM6q093R+wZK5eu/UShWhB6Fyoxg73p1RPc3Xk56u3RJjfLFVnDH/UrVZ9NWMJNfOj7l5U+e0/RogKIEUUJAalNaNisECmjq960Yn3NJv0r62C/LuvLVZCvYdr8kZ8o/u7kqiBSIEvpEoIAmdI/j5aTvkz6oVb5kx/2ShPmz9U/kz75pr6hCa6WU8S994yyug6FBfY7iQPUnUEATxo7vkk76V/osuaj8ROlK9bVfkqlYp6j8vaKEVLqRpilkYoW1UddUfqE+ZflrWRBl7ohAAS3oPsfLSd8vXTxLLig/qaf9yimJFWpdg06UU/L3ypQz59mVq+rZPH9FFaLMhggU0ILpc7yc9D3S1eRZsrR8SR/7NW1erNAoSlBKDeso5POX1GdJB8OJo4Sk3urjI1BAW8M5MfxWkH6NdMna8r1yRvber7mmmoh2UYLNpn6xwuooQQ3lZCY/ar+jYUH5i+pTln9TPUcJikABOzCdjZeT3jA9Sv29uLZ86Wy7634llYYWC1rNonKS+XXuMVGrWrXJeyXccyVm3lvZeSs7qf/6EyigNd3feDnp7dIlNcoXW8Ed9ytVn02vlTP5q/QN1IpaZt5dWX0eQ28OUX8CBTT1+1aMz7mkXyV97JdlXflqshVsu18S8guMjRyG/4bFSB0Mh2hlM45SfwIFNKF7HC8nfZ/0Qa3yJTvulyTMn61/In/2TXu1+tXrk+lpaPnsyq0dqP4ECmjC9DheTvpu6dJZclH5idKV6mu/JFOxTlH5e0UJVeoT5x/ChWQp2nUzxPUpUpq/lgVR5o4IFNCC7nO8nPT90sWz5ILyk3rar5ySWKHWNehEOSV/r02viYfypXDB5mk3ElTRLsdzMQIFtGD6HC8nfY90NXmWLC1f0sd+TZsXKzSKEpRS9ddRWCIoPzP50buvsuAhVL21yr3Vx0eggLaGc2L4rSD9GumSteV75YzsvV9zTTUR7aIEm039YoUeooRIro9BzXu0RFetcs9RgiJQwA5MZ+PlpDdMj1J/L64tXzrb7rpfSaWhRWkrIpU/N7qyiZutozBpZvlzbrAsfsr1HjqPEhSBAtrT/Y2Xk94uXVKjfLEV3HG/UvUpbfVr5a/SN1AraqlYvhmLMsx9tMQu+o8SFIECGvt9K8bnXNKvkj72y7KufDXZCrbdLwn5m+QXByZ6ixUOESUoAgU0onscLyd9n/RBrfIlO+6XJMw/1W4Vlb9XK95bfayiR0vs4ihRglJKT66gifm0Hh/P5EfQZbCvGjP8P1XCwbkvobHn1lnNxe+96UNRWM450hNHo8t6TqdL52X/ezFVzuhoFDTLJfUsTc9/2X9f9+h3V85UCek8xmit7QV0XJ3oZDNR/e8ERndSUnUaUR1VL8wQVS94SWX3zkvXri9h9KoXDhiV26sGp9/kHqVfTf25x5kbtRf0KGymsz6ufekF474qcQyXlEN6f+mSheUn9bRfOfPLqXdKmbiW7ebcteCaO5FfazWMRkwW1qyDYTr47OavoAgUNmfM79+FmTbjxCdO93RRn8Xpxmg1/lIMP6w9PkOZ3/J336/5XL/CnNPF6pPJdOvr1lHwj2fw94rS/V+18v4K29Uzyp/j1X+IGPqa/Jje3z6aDwIFtDWcW8NvBen5dKX6qg/p+fRS41gkFryop6+NxXLm1FOrb6wwc39d/loWRAmz8nvH2csv9jU0ixXWfn42RqCA5oxR/d3fT/qR0pM9t73VU5ZufeRYoda1tVROMr9WZftr8xfUUrZVlGAJWees/LhR3NB5lKAIFLbTrNNq2We35buiQjYcJ/5+lwvL+V4iFW5Xh5uetV39+2HD43CU/Vr29/L+1+t+KVGmVZD+XrVazVEfgLdfYvnj/dJDX4IY0ySPz8x6/7by/XMV5V9ntK38Ik46+mCUban07/s95s3akzQCBTT1+1aE56A66UptW/4onf1iv+T9Sppu1aTy56mfv3R/5x2H5fVZl3++yQ4G/4fFSo//XggU0MRwLTL6VgzfgYrpIxuUH36r2S/2a3K/xlni+iQyDenJVkQqf6tWVqqP/Hc06c6RSvVZmn+xyT6GZQMTBcd/bwQKaMKcdBya/XLp7FcmXeK1wfF20+TySxuUOvUU5isYVTZvQ6pPQdSVzy+8ML9xr/toCSlK6KD7IIFAAS3oVuPQB0hnv46VXmO/crbOnzJRgLRfyYxCfarMbVzQZyDmF/5eCy7XXbiw5u5KKUrooPsg4VTrAO5Oa+2OplbCYmnBW8briC1YaWvZ4lxt3jVamXFIcsdk9K2I0r/bmp0/ruGc/PYX+675+RNHg/1iv4L9UtN+X3+l9XfNwF9K+j3DmcIrZMm9kpl3ffd7+Na6/frulHAcdFSeMd9TQN3qLXuXrd7wh9Qz3zV/W1JQEJRcuq3hqKY/Tto/tW6JQKGm3eemAgCug0ABAIDjSXSmHhmBAgAAEDGZEQAAiAgUAACA6L+9K4Cv9/v9fr+VUvf7/X6/71yb7ZXu7/v9fj6fcfr9fk+mn5Ld04vs7/ydvexnw32J1DU+FaX7e9kPRn0GHbja32XB/j4eDz7DSqnb7bZ3LRqZv7MX/Gy8Xq94Zx+Px9712sqy/b3gB2Mj9Cjsz86Pvd1u7hLq8/ksW+fgEJbtr72SeDweV+huSbrUjhft7AU/G//+/VPel+j9fv/9/f39/amTdi0s298LfjC2snekcnUuUvYTbcrr9dqpUhtavL+32+2sxyTv8XjYfbfO3aOwbGev9tmw+xscnOQ36xwW7+/VPhjbOeGn6ljsRznoQ7M9ZqdsEhbv71lPgpP8hvOsnwpn2c5e7bMhBdZnvcBYvL9X+2Bsh7sedvb5fFTU0Wo70+xLJ1Nlf92Epit4v9/2u5oZcD2N9Tt7hc+GDaeu051eZX+v8MHYDnMUunCd77y1eH+DRbJfr9fVDh0k1/lsJNs8l3i+vV65v9f5YGyHHgUcgH+muN1ubiT7379/p5y6hfn4bCil3u+3ne53hW4nNW9/+WDUtN+oB4yRR9Fs+lmHG6V0aX+TkxhOPHsr6cQzV2Lzd5bPhmssL/XZmNxfPhgVcbx2lm8429dna3X3Nzk18qwIFIpc4bPhX0+f76IiVmV/r/DBqI6hBxyYHWtkmhJip/9s3O93u5CAbfNOP+5ea39P/8HYApMZu/B+v0//PfddbX+Buu73u71LyJx0WbbA1fa3N/Qo7Mz2gwXhrZ1rE9xTfg7L9vd+v2ut44sA97SIqnXEkVzws2EXM7XjMnvXpYVl+3vBD8aG9hrzgMXKjCa1v6+B/VUasb7UZ5g5ChafjRPvWtLM/eWDsR2O1/7cAqX2U+4urPeu11bm7K9N8b/kLsW+y01rus6kJAIFi89G/trvfBcYM/eXD8Z2TtsaHUvc6753jbY1ub/xd/5qj8uLEShYF/9sJHc22XCew/z9vfgHY1OnfUTh4bhHrd/v9yuMny3b36sdJczHZwNJfDDWI1AAAAAi7noAAAAiAgUAACAiUAAAACICBQAAICJQAAAAIgIFAAAgIlAAAAAiAgUAACAiUAAAACICBQAAICJQAAAAIgIFAAAgIlAAAAAiAgUAACAiUAAAACICBQAAICJQAAAAov/2rgAAAFBKqff7/e/fP2OMn/J+v+3P9/v9fr+3rxWBAgAAXXAxgRqCBv/Vv7+/2+3m52mDoQcAALrw9/f3eDzszzZKuN1ur9fLGGPTP59P+04F7XdxAACAvWitX6/X/X5/Pp9/f39KqWAYwkYPjRtuehQAAKjj/X7f73ettW3slVLP59P+MMlmsx0GdnzB9S5Yri+h8egDPQoAAFRwv98/n0/ypTlNrdb68Xi4cOHz+djehSCPUipO3xQ9CgAArPV8Pm2UYKcUuFkFRfw+A2OMFA00nqZAoAAAwFp2SoF/rf98Pm+328y3++MOSe/323YnzC+zFgIFAABW8Zc68NNnzk6wJWQigOfz6W6CaH97JOsoAABQQdzSzx8jsDMS4nR/NQU3g6ExAgUAAFZZeZUvjTu42ZF7hQgWQw8AAKyycnZhctzBRQmv12vHKEHRowAAwEo2UIjvjZzZ0xCPO7h7KHpYwoAeBQAA6ggu/YOHNSTZYCLok3D3UNSr2nIsuAQAwFpu1uHtdrPhgh8lZJpafylGx94JKWHBJQAADuZ+v9sOgM/n8+/fPxslzFlzaZfnPBWhRwEAgGre77d94oNt/m3fQKapnf8wiL0QKAAAsJUgUHBjCrssnbQMQw8AALRg+xher9fj8fh8Pp13JDj0KAAAsBXXo2BnO7p5iHaZhEM0wfQoAACwlWAlJTdvsfMJjD4CBQAAtmIfGK2iGyCTd0X2iUABAACICBQAAGgtuSBjnwgUAADYXDDWcIhBB4tAAQCAzdlAwd4S+X6/P5/PnHUbe8DtkQAAtOCeB6EOteASgQIAAO0caHaCRaAAAABEzFEAAAAiAgUAACAiUAAAACICBQAAICJQAAAAIgIFAAAgIlAAAACi/wPDq6Fd7T4kggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ROOT import gROOT \n",
    "gROOT.GetListOfCanvases().Draw()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}