{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Working with the public 10-year IceCube point-source data\n",
    "=="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial shows how to use the IceCube public 10-year point-source data with SkyLLH."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Disclaimer**\n",
    "\n",
    "    The released 10-year IceCube point-source data can reproduce the published results only within a certain\n",
    "    amount of uncertainty due to the limited instrument response function binning provided in the data release.\n",
    "    The IceCube collaboration is able to reproduce the published results using detailed direct simulation\n",
    "    data, as done for the publication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting the datasets\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we import the dataset definition of the public 10-year point-source data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.datasets.i3.PublicData_10y_ps import create_dataset_collection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The collection of datasets can be created using the ``create_dataset_collection`` function. This function requires the base path to the data repository. It's the path where the public point-source data is stored. The public point-source data can be downloaded from the [IceCube website](http://icecube.wisc.edu/data-releases/20210126_PS-IC40-IC86_VII.zip)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsc = create_dataset_collection(base_path='/home/mwolf/projects/publicdata_ps/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ``dataset_names`` property provides a list of all the data sets defined in the data set collection of the public point-source data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['IC40',\n",
       " 'IC59',\n",
       " 'IC79',\n",
       " 'IC86_I',\n",
       " 'IC86_II',\n",
       " 'IC86_II-VII',\n",
       " 'IC86_III',\n",
       " 'IC86_IV',\n",
       " 'IC86_V',\n",
       " 'IC86_VI',\n",
       " 'IC86_VII']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dsc.dataset_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The individual data sets ``IC86_II``, ``IC86_III``, ``IC86_IV``, ``IC86_V``, ``IC86_VI``, and ``IC86_VII`` are also available as a single combined data set ``IC86_II-VII``, because these data sets share the same detector simulation and event selection. Hence, we can get a list of data sets via the ``get_datasets`` method of the ``dsc`` instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets = dsc.get_datasets(['IC40', 'IC59', 'IC79', 'IC86_I', 'IC86_II-VII'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting the analysis\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The analysis used for the published PRL results is referred in SkyLLH as \"*traditional point-source analysis*\" and is pre-defined:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.analyses.i3.publicdata_ps.time_integrated_ps import create_analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function create_analysis in module skyllh.analyses.i3.publicdata_ps.time_integrated_ps:\n",
      "\n",
      "create_analysis(datasets, source, refplflux_Phi0=1, refplflux_E0=1000.0, refplflux_gamma=2, ns_seed=10.0, gamma_seed=3, kde_smoothing=False, minimizer_impl='LBFGS', cap_ratio=False, compress_data=False, keep_data_fields=None, optimize_delta_angle=10, tl=None, ppbar=None)\n",
      "    Creates the Analysis instance for this particular analysis.\n",
      "    \n",
      "    Parameters:\n",
      "    -----------\n",
      "    datasets : list of Dataset instances\n",
      "        The list of Dataset instances, which should be used in the\n",
      "        analysis.\n",
      "    source : PointLikeSource instance\n",
      "        The PointLikeSource instance defining the point source position.\n",
      "    refplflux_Phi0 : float\n",
      "        The flux normalization to use for the reference power law flux model.\n",
      "    refplflux_E0 : float\n",
      "        The reference energy to use for the reference power law flux model.\n",
      "    refplflux_gamma : float\n",
      "        The spectral index to use for the reference power law flux model.\n",
      "    ns_seed : float\n",
      "        Value to seed the minimizer with for the ns fit.\n",
      "    gamma_seed : float | None\n",
      "        Value to seed the minimizer with for the gamma fit. If set to None,\n",
      "        the refplflux_gamma value will be set as gamma_seed.\n",
      "    kde_smoothing : bool\n",
      "        Apply a KDE-based smoothing to the data-driven background pdf.\n",
      "        Default: False.\n",
      "    minimizer_impl : str | \"LBFGS\"\n",
      "        Minimizer implementation to be used. Supported options are \"LBFGS\"\n",
      "        (L-BFG-S minimizer used from the :mod:`scipy.optimize` module), or\n",
      "        \"minuit\" (Minuit minimizer used by the :mod:`iminuit` module).\n",
      "        Default: \"LBFGS\".\n",
      "    cap_ratio : bool\n",
      "        If set to True, the energy PDF ratio will be capped to a finite value\n",
      "        where no background energy PDF information is available. This will\n",
      "        ensure that an energy PDF ratio is available for high energies where\n",
      "        no background is available from the experimental data.\n",
      "        If kde_smoothing is set to True, cap_ratio should be set to False!\n",
      "        Default is False.\n",
      "    compress_data : bool\n",
      "        Flag if the data should get converted from float64 into float32.\n",
      "    keep_data_fields : list of str | None\n",
      "        List of additional data field names that should get kept when loading\n",
      "        the data.\n",
      "    optimize_delta_angle : float\n",
      "        The delta angle in degrees for the event selection optimization methods.\n",
      "    tl : TimeLord instance | None\n",
      "        The TimeLord instance to use to time the creation of the analysis.\n",
      "    ppbar : ProgressBar instance | None\n",
      "        The instance of ProgressBar for the optional parent progress bar.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    analysis : TimeIntegratedMultiDatasetSingleSourceAnalysis\n",
      "        The Analysis instance for this analysis.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(create_analysis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As source we use TXS 0506+056."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.physics.source import PointLikeSource"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = PointLikeSource(ra=np.deg2rad(77.35), dec=np.deg2rad(5.7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[==========================================================] 100% ELT 0h:00m:14s[                                                          ]   0% ELT 0h:00m:00s\n",
      "[==========================================================] 100% ELT 0h:00m:13s[                                                          ]   0% ELT 0h:00m:00s\n",
      "[==========================================================] 100% ELT 0h:00m:13s[                                                          ]   0% ELT 0h:00m:00s\n",
      "[==========================================================] 100% ELT 0h:00m:13s[                                                          ]   0% ELT 0h:00m:00s\n",
      "[==========================================================] 100% ELT 0h:00m:12s[                                                          ]   0% ELT 0h:00m:00s\n",
      "[==========================================================] 100% ELT 0h:01m:36s\n",
      "[==========================================================] 100% ELT 0h:00m:00s\n"
     ]
    }
   ],
   "source": [
    "ana = create_analysis(datasets=datasets, source=source)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initializing a trial\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the `Analysis` instance was created trials can be run. To do so the analysis needs to be initialized with some trial data. For instance we could initialize the analysis with the experimental data to \"unblind\" the analysis afterwards. Technically the `TrialDataManager` of each log-likelihood ratio function, i.e. dataset, is initialized with data.\n",
    "\n",
    "The `Analysis` class provides the method `initialize_trial` to initialize a trial with data. It takes a list of `DataFieldRecordArray` instances holding the events. If we want to initialize a trial with the experimental data, we can get that list from the `Analysis` instance itself:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "events_list = [ data.exp for data in ana.data_list ]\n",
    "ana.initialize_trial(events_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maximizing the log-likelihood ratio function\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After initializing a trial, we can maximize the LLH ratio function using the `maximize_llhratio` method of the `Analysis` class. This method requires a ``RandomStateService`` instance in case the minimizer does not succeed and a new set of initial values for the fit parameters need to get generated. The method returns a 4-element tuple. The first element is the set of fit parameters used in the maximization. The second element is the value of the LLH ration function at its maximum. The third element is the array of the fit parameter values at the maximum, and the forth element is the status dictionary of the minimizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.random import RandomStateService\n",
    "rss = RandomStateService(seed=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "(fitparamset, log_lambda_max, fitparam_values, status) = ana.maximize_llhratio(rss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "log_lambda_max = 6.572529558548655\n",
      "fitparam_values = [14.58039149  2.1685849 ]\n",
      "status = {'grad': array([-2.09454353e-06,  2.13693588e-04]), 'task': b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH', 'funcalls': 15, 'nit': 9, 'warnflag': 0, 'skyllh_minimizer_n_reps': 0, 'n_llhratio_func_calls': 15}\n"
     ]
    }
   ],
   "source": [
    "print(f'log_lambda_max = {log_lambda_max}')\n",
    "print(f'fitparam_values = {fitparam_values}')\n",
    "print(f'status = {status}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the test-statistic\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the maximum of the LLH ratio function and the fit parameter values at the maximum we can calculate the test-statistic using the `calculate_test_statistic` method of the `Analysis` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TS = 13.145\n"
     ]
    }
   ],
   "source": [
    "TS = ana.calculate_test_statistic(log_lambda_max, fitparam_values)\n",
    "print(f'TS = {TS:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unblinding the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After creating the analysis instance we can unblind the data for the choosen source. Hence, we initialize the analysis with a trial of the experimental data, maximize the log-likelihood ratio function for all given experimental data events, and calculate the test-statistic value. The analysis instance has the method ``unblind`` that can be used for that. This method requires a ``RandomStateService`` instance in case the minimizer does not succeed and a new set of initial values for the fit parameters need to get generated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.random import RandomStateService\n",
    "rss = RandomStateService(seed=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method unblind in module skyllh.core.analysis:\n",
      "\n",
      "unblind(rss) method of skyllh.core.analysis.TimeIntegratedMultiDatasetSingleSourceAnalysis instance\n",
      "    Evaluates the unscrambled data, i.e. unblinds the data.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    rss : RandomStateService instance\n",
      "        The RandomStateService instance that should be used draw random\n",
      "        numbers from.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    TS : float\n",
      "        The test-statistic value.\n",
      "    fitparam_dict : dict\n",
      "        The dictionary holding the global fit parameter names and their best\n",
      "        fit values.\n",
      "    status : dict\n",
      "        The status dictionary with information about the performed\n",
      "        minimization process of the negative of the log-likelihood ratio\n",
      "        function.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(ana.unblind)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ``unblind`` method returns the test-statistic value, the best-fit fit parameter values, and a status dictionary of the minimizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "(ts, x, status) = ana.unblind(rss=rss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TS = 13.145\n",
      "ns = 14.58\n",
      "gamma = 2.17\n"
     ]
    }
   ],
   "source": [
    "print(f'TS = {ts:.3f}')\n",
    "print(f'ns = {x[\"ns\"]:.2f}')\n",
    "print(f'gamma = {x[\"gamma\"]:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the corresponding flux normalization "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default the analysis is created with a flux normalization of 1 GeV$^{-1}$s$^{-1}$cm$^{-2}$sr$^{-1}$ (see `refplflux_Phi0` argument of the `create_analysis` method). The analysis instance has the method `calculate_fluxmodel_scaling_factor` that calculates the scaling factor the reference flux normalization has to be multiplied with to represent a given analysis result, i.e. $n_{\\text{s}}$ and $\\gamma$ value. This function takes the detected mean $n_{\\text{s}}$ value as first argument and the list of source parameter values as second argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flux scaling factor = 1.423e-15\n"
     ]
    }
   ],
   "source": [
    "scaling_factor = ana.calculate_fluxmodel_scaling_factor(x['ns'], [x['gamma']])\n",
    "print(f'Flux scaling factor = {scaling_factor:.3e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, our result corresponds to a power-law flux of:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.423e-15 (E/1000 GeV)^{-2.17} 1/(GeV s cm^2 sr)\n"
     ]
    }
   ],
   "source": [
    "print(f'{scaling_factor:.3e}'' (E/1000 GeV)^{-'f'{x[\"gamma\"]:.2f}'+'} 1/(GeV s cm^2 sr)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluating the log-likelihood ratio function\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes it is useful to be able to evaluate the log-likelihood ratio function, e.g. for creating a likelihood contour plot. Because SkyLLH's structure is based on the mathematical structure of the likelihood function, the `Analysis` instance has the property `llhratio` which is the class instance of the used log-likelihood ratio function. This instance has the method `evaluate`. The method takes an array of the fit parameter values as argument at which the LLH ratio function will be evaluated. It returns the value of the LLH ratio function at the given point and its gradients w.r.t. the fit parameters.\n",
    "\n",
    "In our case this is the number of signal events, $n_{\\mathrm{s}}$ and the spectral index $\\gamma$. If we evaluate the LLH ratio function at the maximum, the gradients should be close to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method evaluate in module skyllh.core.llhratio:\n",
      "\n",
      "evaluate(fitparam_values, tl=None) method of skyllh.core.llhratio.MultiDatasetTCLLHRatio instance\n",
      "    Evaluates the composite log-likelihood-ratio function and returns its\n",
      "    value and global fit parameter gradients.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    fitparam_values : (N_fitparams)-shaped numpy 1D ndarray\n",
      "        The ndarray holding the current values of the global fit parameters.\n",
      "        The first element of that array is, by definition, the number of\n",
      "        signal events, ns.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    log_lambda : float\n",
      "        The calculated log-lambda value of the composite\n",
      "        log-likelihood-ratio function.\n",
      "    grads : (N_fitparams,)-shaped 1D ndarray\n",
      "        The ndarray holding the gradient value of the composite\n",
      "        log-likelihood-ratio function for ns and each global fit parameter.\n",
      "        By definition the first element is the gradient for ns.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(ana.llhratio.evaluate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "llhratio_value = 6.573\n",
      "grad_ns = 0.001\n",
      "grad_gamma = -0.027\n"
     ]
    }
   ],
   "source": [
    "(llhratio_value, (grad_ns, grad_gamma)) = ana.llhratio.evaluate([14.58, 2.17])\n",
    "print(f'llhratio_value = {llhratio_value:.3f}')\n",
    "print(f'grad_ns = {grad_ns:.3f}')\n",
    "print(f'grad_gamma = {grad_gamma:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the `evaluate` method of the `LLHRatio` class we can scan the log-likelihood ratio space and create a contour plot showing the best fit and the 95% quantile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [],
   "source": [
    "(ns_min, ns_max, ns_step) = (0, 80, 0.5)\n",
    "(gamma_min, gamma_max, gamma_step) = (1.5, 4.0, 0.1)\n",
    "\n",
    "ns_edges = np.linspace(ns_min, ns_max, int((ns_max-ns_min)/ns_step)+1)\n",
    "ns_vals = 0.5*(ns_edges[1:] + ns_edges[:-1])\n",
    "\n",
    "gamma_edges = np.linspace(gamma_min, gamma_max, int((gamma_max-gamma_min)/gamma_step+1))\n",
    "gamma_vals = 0.5*(gamma_edges[1:] + gamma_edges[:-1])\n",
    "\n",
    "log_lambda = np.empty((len(ns_vals), len(gamma_vals)), dtype=np.double)\n",
    "for (ns_i, ns) in enumerate(ns_vals):\n",
    "    for (gamma_i, gamma) in enumerate(gamma_vals):\n",
    "        log_lambda[ns_i,gamma_i] = ana.llhratio.evaluate([ns, gamma])[0]\n",
    "\n",
    "# Determine the best fit ns and gamma values from the scan.\n",
    "index_max = np.argmax(log_lambda)\n",
    "ns_i_max = int(index_max / len(gamma_vals))\n",
    "gamma_i_max = index_max % len(gamma_vals)\n",
    "ns_best = ns_vals[ns_i_max]\n",
    "gamma_best = gamma_vals[gamma_i_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.5, 4.0)"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAF5CAYAAAC2tqKTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZxcZZX/8c/pLd2dlSRkJRsQCMsoakQURxFwBGVxGRBcCCMadcQFFUFEQWUJoyOooBIRDQOyyIAQB0QmijtgQH7KnhC2QEL2kK338/ujKmMISXX3uV237q36vl+venW6Uk/dp6tv16lzn+c5j7k7IiIikh91le6AiIiI9I+Ct4iISM4oeIuIiOSMgreIiEjOKHiLiIjkjIK3iIhIzqQSvM3sNDN7yMweNLNrzazZzEaa2Z1mtqj4dZc0+iIiIpJ3ZQ/eZjYR+BQw0933B+qBE4AzgQXuPh1YUPxeREREepHWZfMGoMXMGoBW4HngWGBe8f/nAe9MqS8iIiK5Vvbg7e7PAd8EngGWAevd/VfAWHdfVnzMMmBMufsiIiJSDRrKfYDiWPaxwDRgHfAzM/tAP9rPBmYDDB48+DUzZswoSz9FRKR/7rvvvlXuvutAP+/b3jLYV6/pDre/72/td7j7EQPYpcwpe/AGDgeedPeVAGZ2E/AG4AUzG+/uy8xsPLBiR43dfS4wF2DmzJm+cOHCFLosIiK9MbOny/G8q9Z0c88du4XbN45/YvQAdieT0hjzfgY4yMxazcyAw4BHgFuBWcXHzAJuSaEvIiIiuVf2zNvd7zGzG4H7gS7grxQy6SHADWZ2CoUAf1y5+yIiInngdHtPkicYbmZzgfnuPn+AOpUpaVw2x93PAc7Z7u52Clm4iIjI/3Ggh0TbVa9399kD1J1MSiV4i4iI9EcPiTLvqqfyqCIiIjmjzFtERDLFcbo90WXzqqfgLSIimZNwzLvqKXiLiEimONCt4F2SgreIiGSOMu/SNGFNRESqzXAzm2tmR1e6I+WizFtERDLFIemENa3zFhERSZtWeZem4C0iIpniuCas9ULBW0REssWhW7G7JE1YExERyRll3iIikimFjUmkFAVvERHJGKMbq3QnMk2XzUVEJFMc6PH4Da3zFhERyR2t8xYREUmbLpuXpuAtIiKZUtiYRMG7FAVvERHJnB5X8C5FwVtERDJFmXfvNNtcREQkZ5R5i4hIpjhGt3LLkhS8RUQkczTmXZqCt4iIZIrGvHun4C0iIhljdLsum5eiV0dERCRnFLxFRCRTCruK1YVvqLa5iIhI+hKOeau2uYiISJrcNebdG706IiIiOaPMW0REMqdHS8VKUvAW2YkZ516cqP2j5542QD0RqS2Fdd66MFxK2V8dM9vbzB7Y5vaimX3GzEaa2Z1mtqj4dZdy90VERPKgMOYdvdWCsv+U7v6Yux/g7gcArwE2AzcDZwIL3H06sKD4vYiI1LgBWCpW9dK+bH4Y8IS7P21mxwKHFO+fB9wFnJFyf0R2qm1iV6W7ICKyQ2kH7xOAa4v/HuvuywDcfZmZjdlRAzObDcwGmDx5ciqdlOqy/+nBseu9kh13r/Nix338bI2Vi3RrY5KSUgveZtYEHAN8sT/t3H0uMBdg5syZXoauiexQ/bCOeNtFrQPYE5Haoi1Be5dm5n0kcL+7v1D8/gUzG1/MuscDK1Lsi4iIZFhPjUw8i0ozeJ/IPy6ZA9wKzALmFL/ekmJfpIZs3LM71G5oSzzzfnFcU7itSK3TUrHepfLqmFkr8Fbgpm3ungO81cwWFf9vThp9ERERybtUMm933wyM2u6+1RRmn4uUVf2I9lC7lqbO8DE3D41n7SK1zjFNWOuFKqyJiEjm1Mp67SgFb6l6rcGx69bGePa8qVlj3iJR7iStlDbczOYC8919/gB1K1MUvEVEpNpoP2+RvBs8KJZBtzTGx7xbmzTmLRJn2lWsFwreIhnz5qO/EWr32/mnD3BPRCrDSXzZvOopeEvVG9wUm23e2hDPvKMz1e2SHVYJFqk5WuddmoK35MY+X4nVCp/wlgHuiIiUlWP0aKlYSQreUvWiGXRTXXxXsZbgMRe/uT58TBGpHQreIiKSObpsXpqCt+TG5qmxbLa5PpZBN9fHx7yjx+waqT3ERRxtTNIbBW8REckYo1tLxUpS8JbcGDQsNms8mkEPSjDm3VQX28msYXA82xepFsq8e6dXR0REJGeUeUtuDGmNZd6DguPPjdYTagfxbL+5WZm3CKDL5r1Q8BYRkUxxN10274WCt+RGtF54dL12NGMHGNQTnOGeYA9xkWqi8qil6dURERHJGWXekhtDgzXKo7PGGyw2YzxJ2yQ7mYlUCwftKtYLBW8REckY02XzXih4S25Ea5RHM+/GBJl39JhJxtlFqkVhnbcy71IUvEVEJHNU27w0BW/Jjeja6WgGnSTzbqyLrRGP1kQXkdqi4C1SJd5ad1y47Z09PxvAnogko/28e6fgLbnRknbmnaC2eUNP7JjNDbFjvvjGlaF2IlnVo8vmJSl4S+r2uuDiULs3/csAd0REMskdupV5l6TgLbkRrZQWzaATjXkH20Z/xqe//oZQO5GsytJlczMbDHwP6ADucvdrKtwlXZcQEZHaY2ZXmtkKM3twu/uPMLPHzGyxmZ1ZvPvdwI3u/hHgmNQ7uwMK3iIikimFCWt14Vsf/QQ4Yts7zKweuAw4EtgXONHM9gV2A54tPix+SW4A6bK5pK5jt3TLnEYvYdfjoXZJjtlUF2vXMToT7yciA6bcW4K6++/MbOp2dx8ILHb3JQBmdh1wLLCUQgB/gIwkvQreIiKSKQNQYW20mS3c5vu57j63D+0m8o8MGwpB+3XAd4BLzewdwPwkHRsoCt6Supahecm8Y4VWAOot1rYxmHnbEG1oIrKNVe4+M9BuR58Y3N03Af+WsE8DSsFbREQyxvozdj2QlgKTtvl+N+D5SnSkN6kEbzMbAVwB7E/hisiHgMeA64GpwFPA8e6+No3+SGUNbYll3mln0HWW/ph3dCvRphZl3lJdKrQl6F+A6WY2DXgOOAF4XyU60pu0Ptp8G/ilu88AXgk8ApwJLHD36cCC4vciIlLjthZpid6A4WY218yO3tkxzOxa4M/A3ma21MxOcfcu4FTgDgpx6gZ3fyiNn7m/yp55m9kw4E3AyQDu3gF0mNmxwCHFh80D7gLOKHd/pPKGDmoLtWuuS7c8anTcGqAu2Da8lWijNjSR6pLwsvl6d59d6gHufuJO7r8NuC3JwdOQRua9O7AS+LGZ/dXMrihWqxnr7ssAil/H7Kixmc02s4VmtnDlStVvFhERSWPMuwF4NfBJd7/HzL5NPy6RF6f3zwWYOXNmfBBSMmNIY0eoXdoZdF0F1nmHtxJV5i1VRLuK9S6NzHspsNTd7yl+fyOFYP6CmY0HKH5dkUJfREQkB3qw8I0+jHnnXdkzb3dfbmbPmtne7v4YcBjwcPE2C5hT/HpLufsi2dDakG7mnWSDkajwDPdgu+ZGzTaX6jEARVp6HfPOu7TWeX8SuMbMmoAlFBa71wE3mNkpwDPAcSn1RUREJNdSCd7u/gCwo2o3h6VxfMmW5vrorPHYuG40m00iukY8epVgUL3GvKW6VKhIS26owpqIiGSLa8JabxS8JXUtwcw7ms3WJ6iUFtUUvEoQrW0e3Y1MJIucilVYyw1dlxARkczpKWbfkRuabS4y8KJVxKLZbJLdwaKia8Sj7ZrqlXmLbEOzzUVERNI0AEvFqp6Ct6Qu7X25o7PUuxOMKsX38w5eldCYt1QZBe/SFLwlZPqFF4fbvucdA9gRGRBHjv9EqN3tyy4b4J6IqDxqXyh4S+qiu4NF12uHZ5t7gl3Fon2NjnkHM3aAVe8cFG4rUi6abV6aZpuLiEi10WxzkR3p3K093Da+O1h0Bncwg7YkY97BvgbHyhsSjHkv+fc9w21FysJV27w3Ct4iIpIpmm3eOwVvCRk8rC3cNu1Z49GZ3wm28w6vLY+/NvHx+fZRmqku2aPgXZrGvEVERHJGmbeEDGuJZ96DgrPNo9lsuMJagg/+0bHraF+TjHkzRJm3ZIuWivVOwVtERDLHFbxLUvCWkGGD4pl3OIMOzzaPDl7Hx5GjP2Pa+4ADNDTHroSIlFPCdd7DzWwuMN/d5w9QlzJFwVtERDLFtVSsVwreEjKksSPcNpolhtdrB8Uz9rjwLPW6+GvT1BSvziYilaHgLSIimaMx79IUvCWktT79zDu6XjvartvTr7AW7WuSqxJNjZptLlmj2ea9UfAWEZHMUeZdmoK3hLTUx2cohyulpbxTFwmqltUFdySLjrMnmW3e3KAxb8kWlUftnSqsiYiI5IwybwlpSTDmnfZ67fh4cJIx73TH56PrwwGalHlL1nhhuZjsnDJvERHJnB4sfEP7eYvs2KC6eLaW9u5gedpVLHqVIPwzAo1J6qKLlIGTeMJa1RdpUeYtIiKSM8q8JSRJ5p36rPEcic4HSDLbvKlembdkjdZ590bBW0REMkcT1kpT8JaQ5uCe3JBktnme1nmn+84T3T8coElj3pJBKtJSmoJ3jZtx7sWhdh88boA7IiJS5K7g3RsFbwkZlCDzTntGdXgNdJLZ5tGZ8RWYD6DMWyR/UgneZvYUsAHoBrrcfaaZjQSuB6YCTwHHu/vaNPojIiLZpglrpaWZeb/F3Vdt8/2ZwAJ3n2NmZxa/PyPF/gjQtnt7qF2i2c3h2ubVP4MlOnadZMy7QZm3ZJAmrJVWyXXexwLziv+eB7yzgn0REZEMcbfwrRaklXk78Cszc+Byd58LjHX3ZQDuvszMxuyooZnNBmYDTJ48OaXu1o6hI7aE2iXJvKM1yqOi48hYgtrmwbShEmPeDQmydpFycGonCEelFbwPdvfniwH6TjN7tK8Ni4F+LsDMmTN1IUUkQ/b7Ymy1wkMXnjbAPRF5ieFmNheY7+7zK92ZckgleLv788WvK8zsZuBA4AUzG1/MuscDK9Loi7zU8NZY5h0dt4b0Z41Hs9I8jQRXYsz73pteET6mSG8SZmqqbZ6UmQ02s6Fb/w38C/AgcCswq/iwWcAt5e6LiIjkgGvMuzdpZN5jgZvNbOvxfuruvzSzvwA3mNkpwDOAyn6IiEiBBklLKnvwdvclwCt3cP9q4LByH19KG9qU/lKx6GXs8MSzoCTHixeiSf+1aQxect8yVu+uIpWiCmsiIpI5tXL5O0rBu8YNa4xl3tHMMknbaMYenejWneDNIzwpL/jahEvAEp/s1j0kT1P6JG9UpKU0BW8REckUR5l3bxS8a1xrQyzzbkoy5p2ToiB5KsdaiTHvuub4ckGRkhxQ8C6pkuVRRUREJECZd41rqY9t7dmYpEhLymPX0eP1kGDMOzqjPrqVaAWKtDQM0pi3lI/GvEtT8BYRkexR8C5JwbvGtQYz7+h6ZEg/gw6XR03wM0bfeKJj10k2e4n+LhsbNeYt5VI7ldKiFLxFRCR7lHmXpOBd41rqOkLtkqzzrgVJ1l2nLfq7bGzQmLdIpSh4i4hItrjWefdGwbvGNdfFxrwTbQmaeoW1UDPqKzDdNTp2nWS2efQqQVO9Mm8po2R/ftrPW0REJH2JMu+q389bwbvGRXcHSzS7WTNRBlySOQjRc6BJY94iFaPgXSVe/77/DLV719kD3BERkYGgz/glKXjXuEHBMe9KjLFGx66jklwhCI/rV6Due/SYjcHKbCJ9ouBdkoK3iIhkizYm6ZWCd5V4/vBYFhQd76zEuHXax0xS2zxtiSreBV/XRs02lzJSbfPStKuYiIhIzijzrhKDR24OtYuu104yuzkvs82TjLFH14jH67cn+H0Ex7wb6lRlT8ooH28TFaPgLSIi2aMx75IUvKvELoO3hNqFx7wTzIoOV0oLHzGmVkZ066KZdwVmxkvtyNH2ABWh4C0iItniVP1lczMbDLS5eyhPUPCuEiMGpZx5V/tfVkLhtezR/bwrcCVEtc1F+s7M6oATgPcDrwXagUFmthK4DZjr7ov6+nyabS4iIhljhTHv6C2bfgPsAXwRGOfuk9x9DPDPwN3AHDP7QF+fTJl3lRjW1BZq1xidbZ5kXXF0l6+U/yZ7tNC0JI15S1lV35/f4e6+o5KW+wCHuvt7zKyxr0+m4C0iqTv8TeeH2v3v7740wD2RzKqy4L1t4DazA4D3AccDzwA/2/4xvVHwrhKDGzpC7eJ7a1f/Ou88SfL7iI6XR9ut/MKUUDuRPDOzvSiMeZ8IrARuBN7g7s9Hnk/BW0REsqf6PuM/CvwP8DZ3fybpkyl4V4nB9e2hdpptXh7hSmmV2FUs2NfomPfzb2oNtZMaUp0bk7yHQub9ezP7FYVL5QuiS8U021xERDLHPH7LIne/2d3fC+wL3AV8ClhqZleY2RH9fb5+Zd5m1ujunWa2B7Da3df1sV09sBB4zt2PMrORwPXAVOAp4Hh3X9uvnstLtNSnuy93otnmwXb1Ke/y1Z2jqwuJdhVL+RxoG5Wf11UqqEpPE3ffBFwDXFOMhccBXwB+2Z/n6e/76AVmNhE4D7i0H+0+DTyyzfdnUrhcMB1YUPxeRESkKpnZy7IPd1/j7pe7+6E7e8zO9HfMexhwLHAh0KfF5Ga2G/AO4Hzgs8W7jwUOKf57HoVLCGf0sy+yjea6WObdGKzgrTHv8qjMPumxzLuxLnbu9LRqfbjUpN+Y2X8Dt2w7Yc3MmoA3ArMoFHL5SV+erL/B+y5gorv/zcz6WsbtEgqXBIZuc99Yd18G4O7LzGxMP/shIiJVLKtj1wkcAXwIuNbMpgHrgGagHvgVcLG7P9DXJ+s1eJvZoRRqsa4DHqTw6WGQu/+wD22PAla4+31mdkhfO7Xdc8wGZgNMnjw58hQ1YVAw867I7Obg0HV0rLwSeV7aGXSS/dXjxwz+jM2qiS59kIPZ5ma2O/AlYLi7/2upx7p7G/A94HvFSmqjgS19nTu2vb68H14N/IJC7dXdga8AD/Xx+Q8GjjGzp4DrgEPN7GrgBTMbD1D8umJnT+Duc919prvP3HXXXft4WBERyS1PeOsDM7vSzFaY2YPb3X+EmT1mZovNrOR8LHdf4u6n9OMn29qu092XRQM39O2y+WJ3v7n475/158nd/YsUirBTzLw/7+4fMLNvULi+P6f49Zb+PK+8XHOwRnl4PXKOxrzDs9uTfPDPz8sTnjXeEBzzrmtS5i2Z8BMKE6+v2npHcWXUZcBbgaXAX8zsVgqXti/crv2H3H2niefOmNlnd3D3euC+Ab1sDvzWzE4DLnEfsJ0a5gA3mNkpFOq6HjdAz5tre3/t4nDbj+kVFJFqkizajDazhdt8P9fd577k6d1/Z2ZTt2t3IIWEdQmAmV0HHOvuFwJHJerRP8ws3uYXv38H8BfgY2b2M3f/j748SV+C937A/sAZZnYf8ADwgLv3Nwu/i8KEN9x9NXBYf9pLadEx7/j+0aFmQPrrtfMkmgUnmbsQr28fzNgblXlL7xJOWFvl7jMD7SYCz27z/VLgdTt7sJmNorCS6lVm9sVikO/NKODV7r6x+BznUKhz/ibgPmBggre7v7t4gBb+EchfRz8voYuIiPRZZYaedpRZ7LQnxUT0Y/08xmRg252kOoEp7r7FzPpc57rPS8XcfQuFKmkLe3usxHRO3xJu2xQc826qSG3zWOZd3/f6BS/RrX25yyJama2hQZm3ZNZSYNI23+8GhHb9KuGnwN1mdguFN8OjKCwfGww83NcnUW1zERHJnmSzzYeb2VwzO7qfR/0LMN3MphWLp5wA3Jr4Z9mGu38d+AiF5dfrgI+5+9fcfZO7v7+vz6NdxTJklxGbwm2ju4MlqVEelfYnxmjG3pkgYw/PCQgeMtH+6inXNm+oV4U1KW0ANhhZ7+6zSx7D7FoKlT5Hm9lS4Bx3/5GZnQrcQWGG+ZXu3tel0f3RRaEEhVO4bN5vCt4iIpI9ZS7S4u4n7uT+24DbynVcM/s0hcz7vylcNr/azOa6+3f78zwK3hkysmVzuG18X+5g1hVqlUxdcKy8J7yWXbPiS4nuZNaoMW/pi+qdqnIK8Lri7mKY2UXAn4F+BW+NeYuIiKTH4CU7QnUTmMWrzDtDRjS1hds2plxhLYnoGHQtSDJ2HRWtUR6tEdBYpzFv6V3CMe/hZjYXmO/u83t9dLp+DNxjZjdTCNrvBK7s75MoeIuISPaUecJapbj7t8zsLgp7fxgwqz9lUbdS8M6QIY3xzDu6Xjte0SvUTMokvMNXkmNGzx1l3tKb5LPNM8fMNvDSjyS2zf+5uw/rz/MpeIuIiJSZuw8dyOdT8BYRkeypssx7oCl4Z8iQ+o7eH7QT0clDUUmWUUWXfKW/dCv+7pGnLVOjRVqik+sag1uJSo2p3glrA0LBW0REMqfcFdbyTsE7QwY39HlDmZeJbkySpwwxbSqCUB6asCaSnN6fREREckaZd4Y014Xq0wPxpULRzFuf+sqjEldCovMlwhuTKPOWvtBFwZIUvEVEJFuqcJ33QFPwzpDWuvhs83iRlvAhw7Thx8CLzhivxDErsQ2t5JBOk5J09VNERLLHE9yKS8XM7OjU+50SZd4ZMijRmHd0rW60PGr6n/ui68Oj4r+N2hCdZ9FQgasEUnO0VEz6b+r3vxlqd9pbBrgjIlXmVZ+4ONz2r5edNoA9kXIyNObdGwXvDGm2eK5XiTHPqLQz6ErI03hUfcrvktHZ5k/9dM8B7olkmoJ3SQreIiKSLZpt3isF7zJoHb8x1K4xOGO80DY48zd4vFrInustwc/o1f/Ok/Zs883jQs1EqpKCt4iIZE/1f/5NRMG7DEYP3RRq1xisTw7xHZ7qqz+BzpVoVpqnGvXRWeqdQ/LzM8oA0K5iJSl4i4hI5mhXsdIUvMtgl0FbQu2iVdIg/cwrSZW0RGPJAT01MP6cRPSqTV1wzLshuJ93T3N+VlTIANCfbUl5WtEiIiIiKPMuixFNscw7yZh3fKw0O4PebZvreOJvg+loN7o66ujsMLo6jc72Oro6jM7O4teOwv9v/b/ODmNQSw/v/sQyho18+WtYCzPj8yR6rvogZd414x9lTmUnFLwlE559vJkLP7QXzz/Z3Oc2DU09NDQ6jU09bN5Qz72/GsE5P32MsZPiG7yISDZonXdpZQ/eZtYM/A4YVDzeje5+jpmNBK4HpgJPAce7+9py9ycNwxraQu2SrPPO877cv71pFN8/YyrNg7s5/QeLGDmuk8Ymp6Gxp/C1aevXQqBuaHIaGp1th84fvncI55+8F1969z6cd+OjjJvSXrkfqEbEz7lYO2tU5l1TFLxLSiPzbgcOdfeNZtYI/MHMbgfeDSxw9zlmdiZwJnBGCv2RjGjfYvzonCn86pox7Pu6F/ncZU8wanysROy+B27k6zc8ylfeO4MvHzeD8258hLGTlYGL5JUy79LKHrzd3YGtJccaizcHjgUOKd4/D7iLKgneLfWxoNGUaJ13umd60nHkxf9vMJd8eneWLmrhPac+z/tOX0p9wrNx9/0387XrCwH8nBNncNGtDzN8VPw1rRXRWePx4wUz9gZl3iJbpXLV1MzqzewBYAVwp7vfA4x192UAxa9jdtJ2tpktNLOFK1euTKO7UkZdnca1/zmRLxy9L1s21HPOTx/lg19MHri32n3/zXz5qsdYs7yJ82ftRfuWLAwMiEi/aT/vklKZsObu3cABZjYCuNnM9u9H27nAXICZM2fm4kLKkPrYeGuSMe+6YCKc5prrh+8ZytyzpvLsY60c8u7VzP76MwwZ0U1fP0P29HE98t6v2cRnL32Ciz6yJ9/7wlQ++e0nSHlpebxyXfAMj+7nnkQ0Y4+2q1fmXTuSzzZXkZaB5O7rzOwu4AjgBTMb7+7LzGw8haxcqtD61Q1cfcEkfn39GHbdrZ2zf7yIA/9lfVmPedCRaznx88/x02/sxl4zN3DESTq9RPLCijfZuTRmm+8KdBYDdwtwOHARcCswC5hT/HpLufuSlua62KSrPNWn7kvG3r6ljvlXjOWmy8bTvqWOd3/ieY7/zPO0tsaOWdfPUZ7jP7Wcx+4byo/PncK+B25g8oz+rb/X+vDyiO4fXl+vzFtkqzQy7/HAPDOrp3B99AZ3/4WZ/Rm4wcxOAZ4BjkuhL/0y/cKLQ+0+cuwAdyRnurvhNz8bzbXf2I3Vy5t47VvXctKXnmXS9K1L6NIZh66rg89c8iSfeMt+fPe0Pbjw1odoaMzPBySRmqY/1ZLSmG3+N+BVO7h/NXBYuY9fCa11sTHvJOOW8TW3AxtIt2yq45z3zuDxvw5h+qs28tnLnmC/gzYM6DH6Y/ioLj564ZN8Y/ZezJ87jnd9YlnF+tIXubr6kvYKhzpl3rVES8VKU4U1GTDd3XDxqXuw+P8N5lMXL+Etx61KfaLYjhz09rW89l/WcuO3J/Kmd68KryUXkRQpeJek4F2C7bWx9wftQHTWeJJdxSq9IModfvzVydz7q1348Nef5tDjV1W4Ry918jlP85lDX8FPL5rEJy9ZUunu1LRohbUGjXnXFgXvkir9ni9V4uc/GMcvfjSOo05ZzlEfeqHS3XmZcVPaOWLWC/zuptEsXdz3+ukiIlmkzLuEXYfFMu/obPNo5SmIryuO7iq27Vj5gutHMe+8ybzx6DV8+NylvY6jR2dx9wQ/im+dGf+eTyzjzqvHcNN3JnLad3vPvjtr5JN/dD/vqPD68Loa+YUIuMa8e6PMWxL5w/xd+O7np3LAm9Zz2refpC7DZ9SI0V289X0r+cOtI1m9rLHS3RGRUpJVWKt6yrxLGN2yKdQuOuaddgaU1N2/HMF/njqNGTM3ctaPnqBxUPb/ao46ZTn/c+VYfnnVGN5/xnOV7k6uxSusBa+gaLZ5TUmYeQ83s7nAfHefPzA9ypYM50mSZX+6bQQXfXR39nzFZr48bzHNrfl4Yx07uYNXHbKeX9+wK93x+YEiUm7JMu/17j67WgM3KPMuaURT/ypybdVssV3FkqybjY9d97/d728exXdP253pB2zm3Ksfp3VoPgL3Voe9dyX/8dHp/P2PwzjgTS9WujvSR9HKbCLVSJm39MttV47lkk/uyb6v3ci51+QvcJlYMlEAACAASURBVAPMPHwdza3d/Pl/Rla6KyKyE+bxWy1Q5l3CsIa23h+0A9Ex7ySzzcv9KaynB679xm7c9N2JvPZtazjje0/S1JyPPcS3/3jR1Oy8+tB13HPHCD56IZmeZFeNoleYNOZdQ2po4lmU3rakVx1txiWn7slN353I4Seu4PTLF6UeuAfazMPWs25lE08/0lLprojIjmi2eUnKvEsY3BCrUd5kXaF2lahr3dvuYGuWN3LRR6bz2P1DOOlLz/Cujy/HzHK/49Yr3ljYkvRvfxjOtP1icxskXcq8Rf5BwVt26tG/DOGi2XuyZWM9X5i7iDe8Y22luzRgRk/oZOzkNh67b0iluyIpefXHY7sE3v/90wa4J9Ibo3bGrqMUvEsYWh8b847Wbk4023yAdwB5dOEQzj5uBqMndnDutY8xpZ97YZdSb7HRmm4f2Mxr+gGbeHShgnfa0l4fDrDxxgnhtlIhCt4l1UTw3uNb3wq1+9DbBrgjOfHimga++fE9GDWhg2/84iGG7lKdC6J3338zf7h1FBvX1zNkeHX+jCJ5Za7oXUpNBO+o1rrYeu3o7mB1CZLngRqD7umBSz69O+tWNXLRrQ9XbeAGmLRX4WrC0sdbmPHaWB37WpZ6TfQEqdiW0QPYESm/Gpp4FqXZ5vIS139rIvf/egQf/urT7PFPmyvdnbKaMK0wLLL8mUEV7omISP/URObdNG1DqN2g4O5gjTmabb7tDmB33zGC6y+eyGHvXcWRH1yNVflnu9ETC6sJVi5tqnBPakslapt3DQ43lQrRhLXSaiJ4S++efLiFb506jekHbOJj5z/NAM9/y6RBLU5zazcb1unPQCRzFLxLqol3rbHDYpl3s8Uy72gGnSTPjdY2B1i7ooHz/21PBg/r5qwfLWZQS+/9j84Yz5rWYd1s3lBf6W5IHySZbd7dh3NaskWZd2k1Ebxl57ZsquPrs6azfnUDF970GKPGxT6w5FVzSw/tmxW8RTJHwbukmgjeuwb35W4OjnmHx/RSvlTd1Wn858d3Z8mDrXzpysXs+YrqnqC2I43NPbS3VcdVhGqXJPPuaVR1thpT9ft510Twlpfr6YHvfW53Fi4Ywb/PeYrXvnV9pbtUEfX1Tk/1roYTyafku4Otd/fZA9SbTFLwrkHu8JNzpvD7m3blA194jiM+uKrSXaocA7wGZueJ5I0um5dUE8F7l6bY5eDokq8mYqlckklnfS3S4g7/dd5kbv/xeI6evYz3fmo5lvNNRpLwHjQzJmX1wdc7yWVzGvU7zhPVNu+dBvtqiDtcfcEkbr18PEeevJxZX36mJpaEldLVaTQ26V1CRPKlJjLvYY2xDUbiZU6zFwx6euDH507htivH8baTXuCUr9fGWu7etG+pp6lZk5mqXkP2/ialF6ptXlJNBO9a190Fl585jQXXjeGojyzj5K8o496qbVMdLUM0Y00kazKYA2VKTQTvIcGtPdMuc5qoSMtOonFHm/GtU/fg7ttH8t7TnuOEzz2HKXIDhQ/2m16sZ/AwBe88SLIxCfW6upIr2pikV7kK3n9fvZypV83pd7uTX12GzuSAO3z/jGnc88td+PDXn+aoD71Q6S5lyoa1DXR31TF8dG0VphHJg+C27zUjV8E7Krq1Z2N41ni0SMvAZsQLrh/Nb24czXtPe26ngTta5rQuwXWCnpS3ktyZVc8XNiQZPSF2fki6kswlsXqlcVJdaiJ416KnHmlh7pem8k8Hr+f4056rdHcyafnTha1Ax01pr3BPRORl9HmrpLIHbzObBFwFjAN6gLnu/m0zGwlcD0wFngKOd/e1pZ6robGbMWNe7HcfWutib87x2eahZolsmwlvXF/PnA9PZ/CwLj5/2ZM01mtF4I48v6QZgHFTYnMiyqG7htfcl1OdMu/c0YS10tJ4V+8CPufu+wAHAZ8ws32BM4EF7j4dWFD8XhLq6YGLPzWNlUubOGPuEnbZNTbprhY8+3gLoye00zo0G5fxRaTIKUzaid5qQNkzb3dfBiwr/nuDmT0CTASOBQ4pPmwecBdwRqnnaqzrZuzg/m/vGd1gJP3Z5sk/S13zHxP5y/+O4GPnP82+r92Y+PnKIfpzdg/wWPlTj7QyecaWAX1O6V1d8PeYaMy7Th/Q8kaZd2mpXk81s6nAq4B7gLHFwL41wI/ZSZvZZrbQzBZ2rNMbbSm/vXkkP/vueN72/pUcOWtlpbuTae1bjGcfb2GP/WM7zomIVFJqE9bMbAjw38Bn3P3Fvq41dve5wFyA0fuM9pGD+l+nfJDFMu+012snqW3+8L1D+PZnp7LfQRuYfZ6KsPRmyYOD6ek29nhleYJ3t7KGTKmr0y8kd/QrKymVzNvMGikE7mvc/abi3S+Y2fji/48HVqTRl2r0/JJmzv+3PRk7qZ2zrlisWt198MhfhgAwY2Y2hxZEatnWjUmit1qQxmxzA34EPOLu39rmv24FZgFzil9v6e256q2HEY39z7zjY97Bdd4pZr1rljdy3vv2oa7eOfeqxQzfpQc0Y7lXD989jAm7b2HE6J3Pa+jRR/9MSTbmrd9lrtTQxLOoNC6bHwx8EPi7mT1QvO8sCkH7BjM7BXgGOC6FvlSVjevqOe/9+7BhbQNzbnyc8VNVbKQvujqNB+8eyiHvruF9zKXP3jbz3FC7OxbG2on0RRqzzf/AzlPBw/rzXPXmDGvo/5rc6HrtcKW0YObb1z25ATZvqOeCD8xg2ZPNfOmqx9jzFbE9ywdihnvePHb/YNo21fOKf+5/zQDJp2jWPvl8ZX+VUiuXv6NUYS2H2jbXccGsvVjyYCufv3wxr3jji+hX2Xf3LRhBfUMPr1TwFskuBe+ScvWOX2c9DA3sENZsscvJWdyXe8vGOs6ftTePLRzKZy5bzIFvK1mUToDu7cbO7r1zBDNeu5HmoV1VMyu8pwauoCQa8w62Xfma4eFjSjIZfPvNlOr/i68im16s5+vvn1EI3Jcu5uCj11S6S7mzdFEzSx9v5aAj9dqJZJYDPR6/1YB8Zd54qE55dNZ4U7QSVKhV6V3FXlzTwLnv25unH2nh9B8s5vVvX8u2Uwmiu4PVmj/9YhQAr8to8O7x9FcKdNfAZ/ho3YP2EQPbD5GBkqvgXas62oyvvn9vli5q4Ys/WsTMw9dXuku55A6/v3kU+x30IqPGaw9vkUzLQQJtZu8E3kGhQuhl7v6rtI6dq+Bdh4fWbIcrpQU/rQ/kvtzucPlZU3nib4M568ePV0XgrtR+3ov+Opjnl7RwzEeXVeT4UjnR8fKOwQPcEemzco95m9mVwFHACnfff5v7jwC+DdQDV7j7nJ09h7v/HPi5me0CfBNQ8JaCX141hgXX78rxn3mOA/9lXaW7k2sLrhvDoJZuDj4mm5fMRWQb5S/S8hPgUgpbVgNgZvXAZcBbgaXAX8zsVgqB/MLt2n/I3bdWBj272C41uQreZh6qU56X3cG2b/fgn4dwxVcmM/PQdbzvc8tqck32QNm8oZ4/3DKKNxy9htahsTkQ/VELe1jlaaw8Otu8p2mAOyJ9Vu7M291/V9wsa1sHAovdfQmAmV0HHOvuF1LI0l/ax0IF0TnA7e5+f3l7/FK5Ct61ZPnTTVz4kT0YP7Wdz132JHX5eZ/MpLtuHE3bpnqOOOmFSndFRMpvtJkt3Ob7ucVNrnozEXh2m++XAq8r8fhPAocDw81sT3f/Qf+7GpOr4B0d8260dGeNJ7VxfT1fnzUdd+PsHy9m8LDyZ4rVrKcHbv/xWPZ85Ub2PKDvu4htvz48y7pzVM++EmvS45l3fs6BquIknbC2yt1nBtrt6A9ppz1x9+8A3wkcJzHlcxnT2WFc+OE9WPbUIM6c+wQTdu//0jh5qYV37sLzS1o4erYmqonkQWFXMQ/fElgKTNrm+92A55M8YbnkL/MOjHnXB0cgo7uDRWub9/TAdz83lb//aRifueRJXnHwhlgHqlRkly93+Pn3xrPrbu28/h3Zn6hWiQy6x6v/M3x0AUhPozLviqnMxJG/ANPNbBrwHHAC8L6K9KQX1f9XmyPzLtiNu24axQe+8ByHHre60t2pCn/7w1Aeu28o7/z356nP1UdVEUlguJnNNbOjd/YAM7sW+DOwt5ktNbNT3L0LOBW4A3gEuMHdH0qny/2Tq7czw2kKzByPjnlHM+iIW743npu/P453nLyC935qOdbPYw/ETPRvXLaGmQc085aDW3f6mN/8cTMLH2jj9E+MTHy8cnOHa781gV3GdvCW965IdX/uPI1B14LoOu8auCiRWQkvf69399mlHuDuJ+7k/tuA25IcPA06NTPgjnljufqCKbzp2DXM/tqz4Ut8Sc08oJkTZi/nN3/c8faiv/njZk6YvZyZBzSn3LOY+349nIfvHcp7Pv0cjYN0+VMkNzzhrQbkK/M2D9Upj67XjurPvtwLrtuVK740jZlvXctnv/0k9fVl7Fgv3nJwK9fNHccJs5dz3dxxL8nAtwbu7e/Pqu5uuOrCiYyb0sahJ6zovUHOdVcgRayFsXLqayQSZI6nUaQl12rgry+7fn39aL5/+jQOePM6Pvf9RTQ0VrpHLw3gWzPwvAVugAXXj+apR1o56cznaNRyH5HcMY/fakG+Mm+gkfQy73J+srnzml25/MxpvPJN6znjR4/T1JydM27bAP6xWcP5wbz1uQrcG9fX818XTWSf127g4KPXEtvNvTKTXStRtSztY3ZXYOe0MGXeeTXczOYC8919fqU7Uw65Ct7V4ta545j3tSm8+tC1fP7yRZkK3Fu95eBWPjZrOOddvIazTxuZm8ANcM03JrBhTQOzry7OH8jeyysivSnzhLW8y1XwLsw273/mnfbuYKXaXX/xBK795m684R1rOO3SJ2hsgq1FfbK0J/dv/riZH8xbz9mnjeQH89ZzyMEtuQjgi/5fK7fPG8ORs1awxz/teOJdGnpq4ANDnmbURyuseV0N/CKzyCG4SKhm5Cp4592dP92Va7+5G4e8ZxWf/NaSzK473n6M+5CDW3Ix5t3ZYXznc1MZMaaT95+eyaJIUkP++Z3fCLf9/c9PH8Ce5JQmrJWU0fCxc5H1mmnvDrYjD909lMvPmsIBb16fq8ANpWehp6m3fcB/dukEnn6klbN+/Ditwzv/79HdNXDdvCdBFpyrMeiURS+GTbpFr6mUV3au01axZU8NYs4p0xk3pZ3Pf39xrgL3VjuahZ4lix4YzA2XTORN71qlfc9FqkGydd69VljLu4yGkR2rw2kKzDZP+xPKthn7xnX1nD9rLwDO/slihg33CvSobxY+0FYys94awBc+0JYo++72gR3M2rKxjos/uTsjx3Yw+/ynB+x5k2Xswfr2OcqCo+u88/QzEhzz3jQuA+s+c67cFdbyLlfBO2862ws7hC1/ZhBfu/ZxJkzL9g5hfSl5+paDWzM17t3ZYVw0ezrLn27mq9c9ypDh2j5VpCpozLuk3AXv0Jh3yruDQXGHsM9P5e9/HsbnLl3C/gdtDD+X7FhPD3z3s9N44LfDOfU/l/BPb8jOLmxpz8ROsla7EvtrR1QkYw8esnPIwHaj5jiV2lUsN/LxV5tDV1048f92CHvzu7K/FWXeuMPcs6bwu5tH84Ezn+XwE1ZVuksiIqnJVeZtQFPg41iau4MB3HrFGG763niOPGkFx31qWarHTktvM7/Lobt4Gc0d/uv8Sfzyv8byzn9/nned+jzdJS7IRHcTK/Wc5ZKntdPRvuarJnrw3MnH3j2ZZXjSMe+ql6e/olz4/U2juOKcyRx05Fpmn/dMxXYIq1bucPUFk7jlBxM4YtZyPvDFZyvdJREpB/f4TbPNs8WIVUuLfkLpz+5gAPctGM5ln92DV7zhRb5w6ZM0Rgfbq1w0E94auH/+/Qm87aQX+PB5T2f2w1HqY94Jstm0dySLZ+z5GfPubhrYbtQkzTYvKVfBO8se+vNQvjl7L6bss4Wzr3wik/XK86ynB3705Snc/pNxxcD9VGYDt4gkpAlrvcpZ8PZQtbRy1wxf9NfBXHDy3oyZ0saXr3mU1qFlPVzN6e6CS78wlQXXj+aYjy7jpLP7NxzRHfwEX5FdxaJrpxPNNq+Fset0eX2leyDVrux/fWZ2pZmtMLMHt7lvpJndaWaLil93KXc/yuXJh1r5+gdmMHx0J1+55lGGjeyqdJeqSvuWOi44ZU8WXD+a409b2u/ALSL5ZO7hWy1II/P+CXApcNU2950JLHD3OWZ2ZvH7M3p7IiPdGXa97Sr2zGMtfO3EGbQM7uZr1z/KmAldhAfJ5GXWr27g/H/bk8fuH8zHL3yaQz+4bJvqh9kWHZ/Ny5prSHCVIPraJBjz9mDb6G5kPTm7pplJNRKEo8r+TuHuvwO2X+h8LDCv+O95wDvL3Y+B9uyiZr58/AzqG52vXf8oYyd1VLpLVWXp4mZOP3oGSx5q5YzLn+DIk1ZWuksikpoEM81rJOhX6vPhWHdfBuDuy8xszM4eaGazgdkAEyZmIyt54ZkmvnL8DMyc8254lAm7Z7vsad488LuhXPSxPWhodM674TFmvGZTpbskkqrD33xBqN3//vasAe5Jbg03s7nAfHefX+nOlEPmL+64+1xgLsArXtnokdVX/V3yVUpnh/HNj+9Je1sdF93yMBP3bBuw566EaLGVJJuL7OyY7nDL5eO46vxJ7DZ9C1+a9zhjJ3UknjgWbZ9kuVf6S8WSbAkanSSX7kS3PG1o4vXx7G/a/Hy/pwwIR0vFelGp4P2CmY0vZt3jgRUV6ke//dcFk1j0wBDO+OEiJu2lP7KBsvK5Jn549hTu/dUuvP7ta/jUJUtoGay1IiI1S3/+JVUqeN8KzALmFL/e0pdGhqVe6nRb9/9mOLf+cBxvP/kFXv/2tRXrRzXp7DDm/3Ac1188AXf4t688wzGzl+9wRnmyLTrzIZqxJ9mYJNo29Yw9yUZBaWftCUb4Nk5SbVVIvCVo1St78Daza4FDgNFmthQ4h0LQvsHMTgGeAY4rdz+SWr+6ge+cNo3JMzYz6+xnKt2dqvC3Pwxl7tlTWbqohdcdsYZTvvoMY3bTxD8RoWYmnkWVPXi7+4k7+a/Dyn3sraLZet02H59/+KWpbFzfwNeuXURLi5G1JWFpbxQSLXEKsOyZRn7ytSncc/tIxkxu44s/eYyZh68DyrMZSDRjT5KtpZ7NJiiYkvbSrSTj83mRpH5Nx9Dqf30kucxPWMuCP/3PCP4wfyQfOGMpU/fZUunu5NaWTXXc9L1x3Pz9cVidc+IXnuWY2ctUSlZEXsqBHr0vlJK74J32YrG2zXX88CuTmbbfZt798RdSPnr5JZk13udjdMH/Xj+aa74xkXUrG3njsas46UvPMmpC/y6Rp721Z5IZ4/EiLTka8w62i14l6OqJ/4ypX4FNkDx3tQxcN/KrdtZrR+UueKftxkvHsXp5E6f/YAkNjTqZ+sMdFi4YzrwLduOZx1rY57UbOOtHi9nj1S9WumsiknUK3iXlKngbvZcs3ZHoOu+N6+q59Ydj+edj1rDvazeGnqNa9TbG/uCfh3L1nN14dOFQxk1t4wtzF/H6t6/FDDqDf5TRDUaiKrHOOz5WnuQqQT7KnFaiPGpYsKwqQE+DxryBpMFbRVpq2e3/tSttm+v511OXVborudDRZtx7xy7ccfUY/v6nYYwc18HH5zzJYSes0lULkRQcduiFoXYLfv3FAe5JxalISy377c2j2P/1G5i2X/YnqaUxdr0j7vDE31tZcN2u/P6WUWxc18CuE9s5+cvPcOSsFxjUUvmgHX1lEs02D2azldgStDO4f2X0mPHKbOln3uENTUKtCqKbmoy/u4qKRmnCWq8UvHdi7QuNPPNYCyef/Wylu5JJ61Y28LubR7Hghl15+pFWGgf1cNCRaznsvSv5p4NfpF77GYtImEOFEpK8yF3wHsg65aU8u6gw5XPPV2xO5XhQmey5PzO4O9qMe381gl/fOIr77xpOT7ex5ys3MvvCJ3njMasZPLz7/x5baoZ3eN11qFWS8ecE48jhrDTdjB3iP2dncPZ3V0/sk12Scetw1h5M/pIMsQcvhLB5bFP8oFmkCWsl5S54p2XrGG1Pdy8PrHI9PfDwvUO466ZR/HH+Lmx6sYFR4zo49mPLePN7VjFpr+wPKYiIVBsF750YsWthDfLjDwzmgDdtqHBvehddA72jWeMb19fz8D1D+dsfhvGn/xnJmuVNNLd28/q3r+WQf13F/m94kZ66nmL79PqadqW0ZGun060+Fh23BugMDrKG12un3A6gpyeaCqc/8zv6Y3a2VtEsdY1590rBeycm7N7OzEPXccvl43jHySsZPKx6U/D1qxt4+J6hPHR34fbUw624G42Denj1W9bxz8esYeZb19Hc+o9Qrb8rkerx1jecV+kuvJwum5ek4F3CiZ9/ntOP2ocz3703Z13xBOOntvepXZK632lY+VwTD90zpBCs7xnC0sWF8f2m5m72fs1GTvjcc+x30Ab2etXGspQuja7XTrtSWiVmm3d67E8yyVWC6OzvzpQrpXUnqrCW7ph3gmXe4TKS3YNiP+PoBzbFDlhuCt4lKXiXMP2Vmznn6kV8899357Nv34dPfuMpDjpyHXVp12hNoLsLnnuimUcWFoL1w/cOYcXSQQAMHtbFjJkbOfT4Vex74Ab2eOUmGpv0ByMilabyqL1R8C6hhx5e+eZ1fPP2B5nz4enMmb0nE6a1ccRJL3Do8asYMmLHl9LrUq7A3tlhLH+2kWVPDmLZU80se6qZ5U8NYtmTzaxY2kR3V6E/w0d3st/rNnDMR5ez74EbmLLPZqh76R9IX7PbJHtrh9ddB9tVYrZ5J7Ex6I7g2HWyMe9Y267gMSuReYevokTHyiuwc1p0fXjbWO0fnkcK3n0wdnIH/zH/Yf5020h+OW8MV351CvPOn8T4qe1M2L2NCbu3MXGPNkZPbGf4qC52Gd3N8FFdA1JVzB062422zXWsW9nIsqcLQXnZU4OKt2ZWLm16yYScliHdjJ/Wxu7/tImDj17DxD22sNerNzFh9za2ry5bji04RUQScQpLXWSnFLz7qHGQ8+Z3rebN71rNkw+18Mf5o1i6uJnnlzTz198Op7P95VnBkBFdDB/VyfDRnQwf1cXw0Z00D+6hsbGHts31tG2po31zHW2b64tf62jfwf07mik7eHgXE6a1sfdrNnDIe9oZM6WN8dPaGDe1jWEju14WpKGYufr290VnqcejfjRr7wzPxE53/DlJ2+g670SZd7Bte3fsZ4zOGk9Uv707WmEteMAKfCiOngKZnaWu2uYlKXgHTNtvC9P2W/p/33d3w6rnmli9vIn1qxpZv6qB9asbC7dVDaxf1cjSxc08dPdQ2jbX09luDGrtobm1u/i1h+aWHga1FjL25sHF+1te+phhozoZN6Wd8VPbGLrLSy/Zp71ph4hIWSV7T1Nt82oQzhL72qwORk9qZ/Skwmz0Ujuf9fSAGTvMjPsTgLe/3J322unozmCFtrF24UpgwZQkOv4M6a/Xbu9pDLUrtE03g+7oDo6VdyeYbR4cL7dgxm6VuOIbTKC7m7KYebvWo/aiJoJ3luRpprqIiGSTgncJZc/YB+p4JMmgYylCktnm8Qw6OFYaPV5wxjhAm8fqTLcFM+gkY97RzLsj2C48uz2YsUN8zDs82zxJ5p3u20429w93cG1MUpKCt4iIZI8um5eUu+Cd9eplkCwrjYpm0NHPttFx60Lb2Cf9jmAm3ObRbDbJbPN0+xrN2CGeebd1BTP2YLtEY97Rtl3BMe9opk98vDx8xAwm3oCKtPRCI7AiIiI5k6vM2wkuicrqJ8vtJFnuFc2gO4LHjGbPhbaxz4xt0THWaBacIJvd1DMo1C46a3xLd7yvbcG2bcF13uHZ5p3xMW8PZtD1wf2I6hLsY2TBthWZ4V4u7irS0otcBW8REakRumxeUu6Cd2j/6JycBEnGyqNj0NEMui3B7Oa24FhydDw4mgVH2wFs7onNNt/QHaszvak7QV+7Yn3d3Blr194RrD7XGR/ls47gOu/O4Jh3V6hZorbxjD2b74+uzLuk3AVvERGpdtpVrDe5Ct7uHh6jTVPau19B+nW/o9kzwKbgGujNnm4GHc2eATZGM+iuWF83dMZ3hnqxI/j6dMSuhHREM++O+NWe+o7gmHd7sF1HqFmhbWesXV20XYKrBFI5uQreIiJSAxyt8+5FroJ3D0Z7ivvkpr0PdHjPYRKsgY6u8Q2OP0M8g34xmM2u724NtVvTNSTUDmBtZ/CY0XbtLaF2AC+2x17XLW2xKxPdW2LnXN3meOZdtyV2ham+LXa8+i2xdgANwbYNbbFgV9+e0SCpCmsl5Sp4i4hI9XPAlXmXVNHgbWZHAN8G6oEr3H1Oqcd3Y6xLMAO4v8KZd3T/6AS1tOMVvWLZU5KZ2NEZ1dEMen1XLCtd3Tk41A5gdXssa1/VFsy8N8faAWzcFPt9dG2KXX2p3xg7Vxs2xmebN2yMtWvcHDzepli7wjFjGWfDlmDm3ZbBDNc9aeat/bzLxczqgcuAtwJLgb+Y2a3u/nCl+iQiIlVB+3mX0YHAYndfAmBm1wHHAjsN3m3eyKMd41PqXoIMOuW61hCvzBWdUR3NZgE2dAUz747YMaPjwevaEvyMm2M/Y0cwm2Vj/E+5YVPsPG/ZGLsy1bgh1IymYPZcaBvL4po2BLPgjfEp3PWbYtPG67a0h9rZlgRT48tIl81Lq2Twngg8u833S4HXVagvIiKSJZqwVlIlg/eOPra/7KOWmc0Gtl7+aD95rz8/WNZe5ddoYFWlO5Fhen12Tq/Nzum1KW3vcjzpBtbe8b9+4+gET1H1v7NKBu+lwKRtvt8NeH77B7n7XGAugJktdPeZ6XQvX/TalKbXZ+f02uycXpvSzGxhOZ7X3Y8ox/NWk0puCfoXYLqZTTOzJuAE4NYK/knyqQAABEFJREFU9kdERCQXKpZ5u3uXmZ0K3EFhqdiV7v5QpfojIiKSFxVd5+3utwG39aPJ3HL1pQrotSlNr8/O6bXZOb02pen1qRDzHGz0ISIiIv9QyTFvERERCchc8DazK81shZntcEmYmR1iZuvN7IHi7Stp97FSzGySmf3GzB4xs4fM7NM7eIyZ2XfMbLGZ/c3MXl2Jvqatj69NLZ87zWZ2r5n9v+Lr89UdPKZWz52+vDY1e+5AoSKmmf3VzH6xg/+ryfOm0rK4MclPgEuBq0o85vfuflQ63cmULuBz7n6/mQ0F7jOzO7crKXskML14ex3wfWqj+E1fXhuo3XOnHTjU3TeaWSPwBzO73d3v3uYxtXru9OW1gdo9dwA+DTwCDNvB/9XqeVNRmcu83f13wJpK9yOL3H2Zu99f/PcGCn9ME7d72LHAVV5wNzDCzNKrKVshfXxtalbxfNhaYLSxeNt+wkutnjt9eW1qlpntBrwDuGInD6nJ86bSMhe8++j1xUtct5vZfpXuTCWY2VTgVcA92/3XjsrO1lQQK/HaQA2fO8VLnw8AK4A73V3nTlEfXhuo3XPnEuALwM7qldbseVNJeQze9wNT3P2VwHeBn1e4P6kzsyHAfwOfcfcXt//vHTSpmSyil9emps8dd+929wMoVDM80Mz23+4hNXvu9OG1qclzx8yOAla4+32lHraD+2rivKmk3AVvd39x6yWu4jrxRjNLUgM3V4pjcv8NXOPuN+3gIX0qO1uNenttav3c2crd1wF3AduXoKzZc2ernb02NXzuHAwcY2ZPAdcBh5rZ1ds9pubPm0rIXfA2s3FmZsV/H0jhZ1hd2V6lo/hz/wh4xN2/tZOH3QqcVJwBehCFfW2XpdbJCunLa1Pj586uZjai+O8W4HDg0e0eVqvnTq+vTa2eO+7+RXffzd2nUihh/Wt3/8B2D6vJ86bSMjfb3MyuBQ4BRpvZUuAcChNIcPcfAP8KfNzMuoAtwAleO5VmDgY+CPy9OD4HcBYwGf7v9bkNeDuwGNgM/FsF+lkJfXltavncGQ/MM7N6CoHnBnf/hZl9DGr+3OnLa1PL587L6LypPFVYExERyZncXTYXERGpdQreIiIiOaPgLSIikjMK3iIiIjmj4C0iIpIzCt4iIiI5o+AtIiKSMwreIikws+lm9pSZ7Vn8vrG4ycVule6biOSPgrdICtx9ETAXeFvxrlOBW9x9aeV6JSJ5lbnyqCJV7EHgcDMbCZwCvK7C/RGRnFLmLZKex4G9gXOBb7r7psp2R0TySrXNRVJS3LL0eeAJ4A3u3lPhLolITinzFkmJu3cCLwJnKnCLSBIK3iLpagR+W+lOiEi+KXiLpMTMpgJP1/I+0CIyMDTmLSIikjPKvEVERHJGwVtERCRnFLxFRERyRsFbREQkZxS8RUREckbBW0REJGcUvEVERHJGwVtERCRn/j+8ZFcsoAyg+wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.colors import LogNorm\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.pcolormesh(gamma_edges, ns_edges, log_lambda, norm=LogNorm())\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label(r'$\\log(\\Lambda)$')\n",
    "plt.contour(gamma_vals, ns_vals, log_lambda, [np.quantile(log_lambda, 0.95)])\n",
    "plt.plot(gamma_best, ns_best, marker='x', color='black', ms=10)\n",
    "plt.xlabel(r'$\\gamma$')\n",
    "plt.ylabel(r'$n_{\\mathrm{s}}$')\n",
    "plt.ylim(ns_min, ns_max)\n",
    "plt.xlim(gamma_min, gamma_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the significance (local p-value)\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The significance of the source, i.e. the local p-value, can be calculated by generating the test-statistic distribution of background-only data trials, i.e. for zero injected signal events. SkyLLH provides the helper function ``create_trial_data_file`` to do that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.analysis_utils import create_trial_data_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function create_trial_data_file in module skyllh.core.analysis_utils:\n",
      "\n",
      "create_trial_data_file(ana, rss, n_trials, mean_n_sig=0, mean_n_sig_null=0, mean_n_bkg_list=None, bkg_kwargs=None, sig_kwargs=None, pathfilename=None, ncpu=None, ppbar=None, tl=None)\n",
      "    Creates and fills a trial data file with `n_trials` generated trials for\n",
      "    each mean number of injected signal events from `ns_min` up to `ns_max` for\n",
      "    a given analysis.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    ana : instance of Analysis\n",
      "        The Analysis instance to use for the trial generation.\n",
      "    rss : RandomStateService\n",
      "        The RandomStateService instance to use for generating random\n",
      "        numbers.\n",
      "    n_trials : int\n",
      "        The number of trials to perform for each hypothesis test.\n",
      "    mean_n_sig : ndarray of float | float | 2- or 3-element sequence of float\n",
      "        The array of mean number of injected signal events (MNOISEs) for which\n",
      "        to generate trials. If this argument is not a ndarray, an array of\n",
      "        MNOISEs is generated based on this argument.\n",
      "        If a single float is given, only this given MNOISEs are injected.\n",
      "        If a 2-element sequence of floats is given, it specifies the range of\n",
      "        MNOISEs with a step size of one.\n",
      "        If a 3-element sequence of floats is given, it specifies the range plus\n",
      "        the step size of the MNOISEs.\n",
      "    mean_n_sig_null : ndarray of float | float | 2- or 3-element sequence of\n",
      "                      float\n",
      "        The array of the fixed mean number of signal events (FMNOSEs) for the\n",
      "        null-hypothesis for which to generate trials. If this argument is not a\n",
      "        ndarray, an array of FMNOSEs is generated based on this argument.\n",
      "        If a single float is given, only this given FMNOSEs are used.\n",
      "        If a 2-element sequence of floats is given, it specifies the range of\n",
      "        FMNOSEs with a step size of one.\n",
      "        If a 3-element sequence of floats is given, it specifies the range plus\n",
      "        the step size of the FMNOSEs.\n",
      "    mean_n_bkg_list : list of float | None\n",
      "        The mean number of background events that should be generated for\n",
      "        each dataset. This parameter is passed to the ``do_trials`` method of\n",
      "        the ``Analysis`` class. If set to None (the default), the background\n",
      "        generation method needs to obtain this number itself.\n",
      "    bkg_kwargs : dict | None\n",
      "        Additional keyword arguments for the `generate_events` method of the\n",
      "        background generation method class. An usual keyword argument is\n",
      "        `poisson`.\n",
      "    sig_kwargs : dict | None\n",
      "        Additional keyword arguments for the `generate_signal_events` method\n",
      "        of the `SignalGenerator` class. An usual keyword argument is\n",
      "        `poisson`.\n",
      "    pathfilename : string | None\n",
      "        Trial data file path including the filename.\n",
      "        If set to None generated trials won't be saved.\n",
      "    ncpu : int | None\n",
      "        The number of CPUs to use.\n",
      "    ppbar : instance of ProgressBar | None\n",
      "        The optional instance of the parent progress bar.\n",
      "    tl: instance of TimeLord | None\n",
      "        The instance of TimeLord that should be used to measure individual\n",
      "        tasks.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    seed : int\n",
      "        The seed used to generate the trials.\n",
      "    mean_n_sig : 1d ndarray\n",
      "        The array holding the mean number of signal events used to generate the\n",
      "        trials.\n",
      "    mean_n_sig_null : 1d ndarray\n",
      "        The array holding the fixed mean number of signal events for the\n",
      "        null-hypothesis used to generate the trials.\n",
      "    trial_data : structured numpy ndarray\n",
      "        The generated trial data.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(create_trial_data_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first we will generate 10k trials and look at the test-statistic distribution. We will time the trial generation using the ``TimeLord`` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.timing import TimeLord\n",
    "tl = TimeLord()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[==========================================================] 100% ELT 0h:07m:31s\n",
      "TimeLord: Executed tasks:\n",
      "[Generating background events for data set 0.]   0.002 sec/iter (10000)\n",
      "[Generating background events for data set 1.]   0.003 sec/iter (10000)\n",
      "[Generating background events for data set 2.]   0.003 sec/iter (10000)\n",
      "[Generating background events for data set 3.]   0.005 sec/iter (10000)\n",
      "[Generating background events for data set 4.]   0.024 sec/iter (10000)\n",
      "[Generating pseudo data.                     ]   0.030 sec/iter (10000)\n",
      "[Initializing trial.                         ]   0.030 sec/iter (10000)\n",
      "[Create fitparams dictionary.                ] 1.0e-05 sec/iter (593990)\n",
      "[Calc fit param dep data fields.             ] 2.9e-06 sec/iter (593990)\n",
      "[Get sig prob.                               ] 1.8e-04 sec/iter (593990)\n",
      "[Evaluating bkg log-spline.                  ] 2.6e-04 sec/iter (593990)\n",
      "[Get bkg prob.                               ] 3.2e-04 sec/iter (593990)\n",
      "[Calc PDF ratios.                            ] 6.2e-05 sec/iter (593990)\n",
      "[Calc pdfratio values.                       ] 8.2e-04 sec/iter (593990)\n",
      "[Calc pdfratio value product Ri              ] 3.5e-05 sec/iter (593990)\n",
      "[Calc logLamds and grads                     ] 2.9e-04 sec/iter (593990)\n",
      "[Evaluate llh-ratio function.                ]   0.004 sec/iter (118798)\n",
      "[Minimize -llhratio function.                ]   0.052 sec/iter (10000)\n",
      "[Maximizing LLH ratio function.              ]   0.052 sec/iter (10000)\n",
      "[Calculating test statistic.                 ] 3.5e-05 sec/iter (10000)\n"
     ]
    }
   ],
   "source": [
    "rss = RandomStateService(seed=1)\n",
    "(_, _, _, trials) = create_trial_data_file(\n",
    "    ana=ana,\n",
    "    rss=rss,\n",
    "    n_trials=1e4,\n",
    "    mean_n_sig=0,\n",
    "    pathfilename='/home/mwolf/projects/publicdata_ps/txs_bkg_trails.npy',\n",
    "    ncpu=8,\n",
    "    tl=tl)\n",
    "print(tl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After generating the background trials, we can histogram the test-statistic values and plot the TS distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU9fX/8dcxqOACKotEAgJClbAIGI2orUhFoYpgf4hrK7UVxa31Z6vUqoj9tVK16letC1ZxB60tIIpWvwKuNChIlUWWAmLYZBFxQ7bz+2Myt8MwM5kkM7kzyfv5eOSRzJ25956EkDOf9Zi7IyIiArBH2AGIiEjuUFIQEZGAkoKIiASUFEREJKCkICIigQZhB1ATzZo187Zt24YdhohIXpk1a9Z6d2+e6Lm8Tgpt27bl/fffDzsMEZG8YmafJHtO3UciIhLIy6RgZgPMbMwXX3wRdigiInVKXiYFd5/s7sOaNGkSdigiInVKXo8piNSWbdu2UV5ezpYtW8IORSRtDRs2pKioiD333DPtc5QURNJQXl7O/vvvT9u2bTGzsMMRqZS7s2HDBsrLy2nXrl3a5+Vl95FIbduyZQtNmzZVQpC8YWY0bdq0yq1bJQWRNCkhSL6pzu9sXiaFms4+GjV5Hmc/NINnylZkODIRkfyWl0khE7OPypZtZNKclRmMSiS7li9fTpcuXWp0jenTp3P66adnKKLMGjp0KM8//3xGrrVw4UIWLlyYkWvVN/VyoHnkgM7MX7U57DBE8oq74+7ssUdevpeUNOlfVySPbN++nQsvvJBu3boxePBgvvnmG2655RaOPvpounTpwrBhw4hWU1yyZAknn3wyRx55JD179uQ///nPLtd677336NGjB0uXLmXdunX07duXnj17cskll3DooYeyfv16li9fTqdOnbjsssvo2bMnn376KePGjaNr16506dKF6667LrjefvvtF3z9/PPPM3ToUCDSArjqqqs47rjjaN++fdAacHeuuOIKiouLOe200/jss8+y/NOTdNTLloJITYyaPC/jLc3iQxozckDnSl+3cOFCHnnkEY4//nguuugi7r//fq644gpuuukmAH7yk5/w4osvMmDAAM4//3xGjBjBmWeeyZYtW9i5cyeffvopAO+++y5XXnklkyZNok2bNlxxxRX06dOH3/72t7zyyiuMGTNml3uOHTuW+++/n1WrVnHdddcxa9YsDjzwQE455RQmTpzIoEGDUsa9evVq3n77bT7++GPOOOMMBg8ezIQJE1i4cCEfffQRa9eupbi4mIsuuqgGP0XJBLUURPJI69atOf744wG44IILePvtt5k2bRqlpaV07dqVqVOnMm/ePL788ktWrlzJmWeeCUQWMe2zzz4ALFiwgGHDhjF58mTatGkDwNtvv80555wDQL9+/TjwwAODex566KEce+yxQKR10bt3b5o3b06DBg04//zzefPNNyuNe9CgQeyxxx4UFxezdu1aAN58803OPfdcCgoKOOSQQ+jTp0+GfkpSE2opiFRROu/osyV+iqGZcdlll/H+++/TunVrbr75ZrZs2RJ0ISVSWFjIli1b+OCDDzjkkEMAUr5+3333Db5O9brY2OLnxu+9994Jr6FpvrknL1sK2hBP6qsVK1YwY8YMAMaNG8cJJ5wAQLNmzfjqq6+C/vrGjRtTVFTExIkTAfjuu+/45ptvADjggAN46aWXuP7665k+fToAJ5xwAs899xwAr776Kp9//nnC+5eWlvLGG2+wfv16duzYwbhx4zjxxBMBOPjgg1mwYAE7d+5kwoQJlX4vP/jBDxg/fjw7duxg9erVTJs2rZo/FcmkvEwK2hBP6qtOnTrx+OOP061bNzZu3Mjw4cO5+OKL6dq1K4MGDeLoo48OXvvkk09yzz330K1bN4477jjWrFkTPHfwwQczefJkLr/8csrKyhg5ciSvvvoqPXv25OWXX6awsJD9999/t/sXFhZy6623ctJJJwUD2AMHDgRg9OjRnH766fTp04fCwsJKv5czzzyTjh070rVrV4YPHx4kFwmXpWoO5rqSkhKvbpGdsx+KvNt69pJemQxJ6qgFCxbQqVOnsMPImu+++46CggIaNGjAjBkzGD58OHPmzAk7rGqLrlE4/PDDQ44kfIl+d81slruXJHq9xhREhBUrVjBkyBB27tzJXnvtxcMPPxx2SBISJQURoWPHjnzwwQdhhyE5IC/HFEREJDuUFEREJKCkICIiASUFEREJ5GVS0OI1EZHsyMukoMVrUt9s2LCB7t270717d1q2bEmrVq2Cx6NGjaJz585069aN7t27U1ZWFpw3ePBgli5dSmlpKd27d6dNmzY0b948OPejjz7isMMOY/HixQBs27aNrl27Btf4wx/+kPTaURs3bqRv37507NiRvn37Bquhly9fTqNGjYJ7XXrppcE5s2bNomvXrnTo0IGrrrpql60vnnvuOYqLi+ncuTPnnXde2j+jZcuWUVpaSseOHbn66qvZunUrEKkh0aRJkyCOW265JThn06ZNDB48mCOOOIJOnToFq8V//etfM3Xq1LTv/fHHH9OrVy/23ntv7rjjjuD4li1bOOaYYzjyyCPp3LkzI0eOrNL5UTt27KBHjx671MK4+eabd/k9mDJlStrxphTdIz0fP4466iivriEPvutDHny32udL/TJ//vywQwiMHDnSb7/9dnd3f/fdd/3YY4/1LVu2uLv7unXrfOXKle7uPnfuXB80aNAu544dO9Yvv/zyXY49++yz3rdvX3d3/+Mf/+jDhg2r9NqxfvOb3/itt97q7u633nqrX3vtte7uvmzZMu/cuXPC7+Hoo4/2d99913fu3On9+vXzKVOmuLv7okWLvHv37r5x40Z3d1+7du1u544dO9ZHjhy52/GzzjrLx40b5+7uZ599dvCaadOm+WmnnZYwjp/+9Kf+8MMPu7v7d999559//rm7uy9fvjz4maRj7dq1PnPmTL/++uuDfxt39507d/qXX37p7u5bt271Y445xmfMmJH2+VF//vOf/dxzz93l+4j9PUgl0e8u8L4n+buqdQoi1dC7d++MXi+6B1FVrV69mmbNmgUbzjVr1ix47umnnw62oEhlyJAhPProo9x22208+OCDwXqFVNeONWnSpCD+Cy+8kN69e/OnP/0pZcybN2+mV6/IbgI//elPmThxIv379+fhhx/m8ssvD3ZpbdGiRaXxQ+TN7dSpU3nmmWeAyK6s9913HzfffHPSczZv3sybb77JY489BsBee+3FXnvtBUR2ht2wYQNr1qyhZcuWld6/RYsWtGjRgpdeemmX42YW1JnYtm0b27ZtS7gJYLLzAcrLy3nppZf43e9+x5133llpLDWVl91HIhJxyimn8Omnn/K9732Pyy67jDfeeCN47p133uGoo45K6zp333031113HTfccAMHHXRQpdeOtXbt2mCvo8LCwl2K5SxbtowePXpw4okn8tZbbwGwcuVKioqKgtcUFRWxcmWkNO6iRYtYtGgRxx9/PMceeyyvvPJKWvFv2LCBAw44gAYNIu9zW7ZsuUscM2bM4Mgjj6R///7MmzcPgKVLl9K8eXN+9rOf0aNHD37xi1/w9ddfB+f07NmTd955B4Crr7466KaJ/Rg9enSlse3YsYPu3bvTokUL+vbtS2lpaVrfU9SvfvUrbrvttoQV7+677z66devGRRddlHQTw6pSS0GkGqr7zj7T9ttvP2bNmsVbb73FtGnTOPvssxk9ejRDhw5l9erVNG/ePK3rvPLKKxQWFjJ37ty0rp2OwsJCVqxYQdOmTZk1axaDBg1i3rx5Cbffjr573r59O4sXL2b69OmUl5fz/e9/n7lz57Jjxw5++MMfApExjK1btwY7wD755JMp38337NmTTz75hP32248pU6YwaNAgFi9ezPbt25k9ezb33nsvpaWl/PKXv2T06NH8/ve/ByLv3letWgXAXXfdldb3nEhBQQFz5sxh06ZNnHnmmcydOzftWtsvvvgiLVq04Kijjtrtd2748OHceOONmBk33ngj11xzDY8++mi144xSS0EkzxUUFNC7d29GjRrFfffdx9///ncAGjVqtFtdg0RWrVrFPffcw8yZM5kyZQoffvhhpdeOdfDBB7N69Wog0jUU7fLZe++9adq0KQBHHXUUhx12GIsWLaKoqIjy8vLg/PLy8qCuQ1FREQMHDmTPPfekXbt2HH744SxevJimTZsyZ84c5syZwy233MKll14aPO7atSvNmjVj06ZNbN++HYA1a9YEcTRu3DjowvnRj37Etm3bWL9+PUVFRRQVFQXv3AcPHszs2bODuLZs2UKjRo2AmrUUog444AB69+6ddusHIq29F154gbZt23LOOecwdepULrjgguDnXlBQwB577MHFF1/MzJkz075uKkoKInls4cKFwcwhgDlz5nDooYcCkW22lyxZUuk1rr76aq6//nqKioq48847ufzyy3H3lNeOdcYZZ/D4448D8PjjjwfjGOvWrWPHjh1ApKtm8eLFtG/fPtiW+1//+hfuzhNPPBGcM2jQoKCuwvr161m0aBHt27ev9HswM0466aSgnsTEiRODlsWaNWuC1snMmTPZuXMnTZs2pWXLlrRu3TrYUfX111+nuLg4uOaiRYuCd/R33XVXkIRiP0aMGJEyrnXr1rFp0yYAvv32W/73f/+XI444otLvJ+rWW2+lvLyc5cuXM378ePr06cNTTz0FECRigAkTJqTd+qiMuo9E8thXX33FlVdeyaZNm2jQoAEdOnQI6iufdtppTJ8+nZNPPjnp+a+99horVqzg5z//OQADBgzg4Ycf5oknnqBLly5Jrx1rxIgRDBkyhEceeYQ2bdrwt7/9DYiU27zpppto0KABBQUFPPjgg8F4xQMPPMDQoUP59ttv6d+/P/379wfg1FNP5dVXX6W4uJiCggJuv/32oLVRmT/96U+cc8453HDDDXTo0IHBgwcD8Pzzz/PAAw/QoEEDGjVqxPjx44PuqnvvvZfzzz+frVu30r59e8aOHQtEBoWXLFlCSUnC3aV3s2bNGkpKSti8eTN77LEHd999N/Pnz2f16tVceOGF7Nixg507dzJkyJBgWumDDz4IwKWXXpr0/MaNGye957XXXsucOXMwM9q2bctDDz2UVqyVUT0F1VOQNORjPYVvv/2Wk046iXfeeYeCgoKww6lVNa2nMGHCBGbPnh2ML+SzqtZTUPeRSB3VqFEjRo0aFczskfRt376da665JuwwQqHuI5E0uXveFZo/9dRTww4hL5111llhh5AR1ekJUktBJA0NGzZkw4YN1fpPJhIGd2fDhg00bNiwSueppSCShug0ynXr1oUdiqRhzZo1AOzcuTPkSMLVsGHDXRYKpiNnkoKZdQJ+CTQDXnf3B0IOSSQQnTcv+WH48OFA7iwyzCdZ7T4ys0fN7DMzmxt3vJ+ZLTSzJWY2AsDdF7j7pcAQIL15YCIiklHZHlN4DOgXe8DMCoC/AP2BYuBcMyuueO4M4G3g9SzHJSIiCWQ1Kbj7m8DGuMPHAEvcfam7bwXGAwMrXv+Cux8HnJ/smmY2zMzeN7P31b8rIpJZYYwptAI+jXlcDpSaWW/gx8DeQNJqEe4+BhgDkcVr2QtTRKT+CSMpJJro7e4+HZie1gXMBgADOnTokMGwREQkjHUK5UDrmMdFwKqqXMBVjlNEJCvCSArvAR3NrJ2Z7QWcA7wQQhwiIhIn21NSxwEzgMPNrNzMfu7u24ErgH8CC4Dn3H1eFa87wMzGfPHFF5kPWkSkHsvqmIK7n5vk+BRSDCancd3JwOSSkpKLq3sNERHZnfY+EhGRgJKCiIgE8jIpaExBRCQ78jIpaEqqiEh25GVSEBGR7MjLpKDuIxGR7MjLpKDuIxGR7MjLpCAiItmhpCAiIgElBRERCeRlUsjUQPP81Zs5+6EZPFO2IkORiYjkt7xMCpkYaB7YvRXFhY2Zv3ozk+aszGB0IiL5Ky+TQiacV9qGZy/pRXFh47BDERHJGfU2KYiIyO6UFEREJJCXSUErmkVEsiMvk4JWNIuIZEdeJoVM09RUEZGIrJbjzAcDu7cCIokBIrOSRETqq3rfUoidmqoWg4jUd/W+pRClFoOIiFoKAS1mExHJ06SgKakiItmRl91H7j4ZmFxSUnJxNq4fHVuIGti9lbqTRKReyMukkE3RsYUojTGISH2ipBDnvNI2uySAsx+aEbQc1GIQkbpOSaESmpUkIvVJpQPNZna8mb1mZovMbKmZLTOzpbURXC7QrCQRqU/SaSk8AlwNzAJ2ZDec3KZuJBGp69JJCl+4+8tZjyTHqRtJROqDdNYpTDOz282sl5n1jH5kPbIcE9uNVLZso7bCEJE6KZ2WQmnF55KYYw70yXw4uW9g91aULdvI9RM+YtKclepKEpE6pdKk4O4n1UYgVWFmA4ABHTp0qPV7RxPApDkr1ZUkInWOuXviJ8wucPenzOz/Jnre3e/MamRpKCkp8ffffz+0+0fXMBQXNlaLQSSH9O7dG4Dp06eHGkeuMrNZ7l6S6LlULYV9Kz7vn/mQ6obo4HPZso2ULduo7iQRyXtJk4K7P1TxeVTthZNfoqufnylboe4kEakTKh1TMLP2wP8AxxIZYJ4BXO3u9WYBW2WiySG6iV40ScRSC0JE8kE6U1KfAZ4DCoFDgL8B47IZVL6LbTVAZG1DfJIQEclF6SQFc/cn3X17xcdTRFoMkkJxYWOevaSXtsgQkbyStPvIzA6q+HKamY0AxhNJBmcDL9VCbCIiUstSjSnMIpIErOLxJTHPOfD7bAWVz6LdRmodiEg+SjX7qF1tBlIXxBboiS/WIyKSD1RPIYPiC/TESrTDanSWkmYmiUiuUFKoBfE7rEJkhlLZso3BYyUFEckFKZOCmRlQ5O6f1lI8dVLsOob5qzcHyaC03UG7JAYRkbClTAru7mY2ETiqNoIxs0HAaUAL4C/u/mpt3Le2xI85xC54ExHJBel0H/3LzI529/eqcwMzexQ4HfjM3bvEHO9HZKV0AfBXdx/t7hOBiWZ2IHAHUKeSQrIxh+h4Q5TGGEQkLOksXjuJSGL4j5l9aGYfmdmHVbjHY0C/2ANmVgD8BegPFAPnmllxzEtuqHi+zhvYvdUu01e1+llEwpROS6F/TW7g7m+aWdu4w8cAS6L7J5nZeGCgmS0ARgMvu/vsRNczs2HAMIA2bfL/3XR86yE67qBa0CIShnSK7HxiZicAHd19rJk1B/ar4X1bAbGD1+VEKrxdCZwMNDGzDu7+YIJ4xgBjIFJPoYZx5JxkM5Vin4t9rKQhIpmUzi6pI4mU4jwcGAvsCTwFHF+D+1qCY+7u9wD31OC6eS/ZTKX9G/73nyo2YSgpiEgmpdN9dCbQA5gN4O6rzKymhXfKgdYxj4uAVemeHGY5ztoSP1MpdufV2DGI2G261XIQkZpKJylsrZia6gBmtm9lJ6ThPaCjmbUDVgLnAOele7K7TwYml5SUXJyBWHJSqtXR0SQBu27THa3+Fn2NEoSIVFU6SeE5M3sIOMDMLgYuAh5O9wZmNg7oDTQzs3JgpLs/YmZXAP8kMiX1UXefV+Xo65H4JJGo5VC2bKO6lkSkRtIZaL7DzPoCm4HvATe5+2vp3sDdz01yfAowJd3rxKoP3UeVSdS9BNqdVURqJt29jz4CGhHZMvuj7IWTnvrQfVSZZN1LsQlCRKSqKl28Zma/AGYCPwYGE1nIdlG2A5OqOa+0Dc9e0ktdRiJSI+m0FH4D9HD3DQBm1hR4F3g0m4Glou4jEZHsSCcplANfxjz+kl0XntU6dR+lR9NVRaSq0kkKK4EyM5tEZExhIDDTzP4vgLvfmcX4pAZiZyiBZiOJSOXSSQr/qfiImlTxuaYL2KQWJFvoFqUWhIjESmdK6qjaCKQqNKZQdfHFfaLHQC0IEfmvvCzHqTGFysUubEtU3AdQgR8R2U1eJgVJLVESSNUaULeSiEQpKdRBlSWBWOpWEpFY6WydfRvw/4BvgVeAI4FfuftTWY4tVUwaU8gAdSuJSLx0WgqnuPu1ZnYmkTULZwHTiNRUCIXGFDIjVYsitm50oq4krYEQqZvSSQp7Vnz+ETDO3TeaJaqRI3VFbAsiWVeS1kCI1E3pJIXJZvYxke6jyyrKcW7JblgSptgWRKquJO3IKlL3VLohnruPAHoBJe6+DfiayKpmqSeiXUnPlK0IOxQRybKkLQUz+3GCY7EP/5GNgNKhgebaE+1KKlu2MajsFtu9JCJ1S6ruowEpnnNCTAoaaK490a6k6MBy7DiCiNQ9SZOCu/+sNgOR3BZNDsnGGDQbSaRuSGvxmpmdBnQGGkaPufst2QpK8o9mI4nUDeksXnsQ2Ac4CfgrkeprM7Mcl+Sw2H2VEj0WkfyVTkvhOHfvZmYfuvsoM/szIY4nSLjiV0HHP5eJ+tDqihIJTzpJ4duKz9+Y2SHABqBd9kKSXJZoFXTs40wkBXVFiYQnnaTwopkdANwOzCYy8+ivWY2qEpqSWvepK0okHOksXvu9u29y978DhwJHuPuN2Q8tZUyT3X1YkyZNwgxDRKTOSbV4rY+7T022iM3dNa4gaYuOE2iMQCS3peo+OhGYSuJFbKEuXpPcFrvDKvx3ADpat0FJQSR3pVq8NtLM9gBedvfnajEmyWPxM5Ki22Ps31D1nETyQcr/qe6+08yuAJQUJC3xs5Nip5eWLdsYtCLUjSSSm9J5+/aamf0aeJbIDqkAuPvGrEUldUZskogmiPjN9ZQcRHJHOknhoorPl8ccc6B95sORuizZ5npKCiK5I52k0MnddymqY2YNk71YpDLJNteLTRZapyASjkrXKQDvpnlMpEZiE0Kimg3PlK3g7IdmqOCPSBalWqfQEmgFNDKzHkC0wk5jIhvkhUYrmuuu4sLGPHtJL2D3LTO0/YVI9qXqPjoVGAoUAX/mv0lhM3B9dsNKTUV26i91K4lkV6p1Co8Dj5vZ/6nY4kJEROq4tAaao1+Y2d7u/l0W4xHJCG2/LVI9SQeazexaM+tFpKhOVOJajCI5Jjr+MH/15oxs5y1SX6RqKSwEzgLam9lbwAKgqZkd7u4LayU6kRrQ+INI1aVKCp8TGVDuXfHRicjg84iKxHBc1qOTeqE66xPKlm3kmbIVu62Wht33XxKR9KVap9APeAk4DLgTOAb42t1/poQgmVTZ+oR40dfEdgupu0gkM1LNProewMz+DTwF9ACam9nbwOfunmhLbZFqiV2fECt2A72o80rbJPzDn6yVkWgTPg1EiySWzuyjf7r7e8B7Zjbc3U8ws2bZDkzqh+hitER/0KOJIHb77aqOE0SvEb/PkhbCiSRWaVJw92tjHg6tOLY+WwFJ/RH77j9Rt1H8BnrJXpdKsn2WQAPRIolUqfKJu/87W4FI/RNfe6Gmr6tMoq4oEdmVymFJvRDfjSQiiaWzS2qtMLP2ZvaImT0fdiySH6Lv/NPZMfW80jY8e0mvhF1G0emtldEurVIfZDUpmNmjZvaZmc2NO97PzBaa2RIzGwHg7kvd/efZjEfqjoHdW1Fc2LjGU1ATTW9NRtNepT7IdkvhMSLrHQJmVgD8BegPFAPnmllxluOQOibVO/+qXqe03UFpv764sLEGqKVOy+qYgru/aWZt4w4fAyxx96UAZjYeGAjMz2YsUn+kmuaajmyuYdD6CMl1YQw0twI+jXlcDpSaWVPgD0APM/utu9+a6GQzGwYMA2jTRv+hZFeVTXNNRzbXMGh9hOS6MJKCJTjm7r4BuLSyk919DDAGoKSkxDMcm+S56k5fjQ5aR7/OZheRup8kl4Ux+6gcaB3zuAhYVZULmNkAMxvzxRdfZDQwqZ+ig9ZR6e7BJFIXhdFSeA/oaGbtgJXAOcB5VbmAynFKJiVrXWiGkdRHWU0KZjaOyLbbzcysHBjp7o+Y2RXAP4EC4FF3n5fNOKRuq+nAcnUl2n5j0pyVGkCWvJbt2UfnJjk+BZhS3eua2QBgQIcOHap7CakjMjGwXF3xg8YQWQgHGkCW/JWX21yo+0iiMrUvUnVp0FjqmrxMCiK1qaqV4RLVb8hkHKA1DpI9ObP3UVVo9pHUpqpWhsvE9hup4tA2G5JNeZkU3H2yuw9r0qRJ2KFIPRGtDFfZu/PSdgdlZPuNVHGoy0qyKS+TgoiIZIfGFERSyNR012TjErHjBFEaL5Aw5WVLQWMKUhuiK50zscI52bhE/LRWjRdI2PKypaApqVIbMj3dNTouAbuulo49nqiWtEhtysuWgoiIZIeSgkgWRUt41qQ2dKJrVKUUqUhV5GX3kba5kHxR1TUOVblGNEloUFoyKS+TgsYUJJ/Ejhlk6hrnlbbR+INkhbqPREQkoKQgIiIBJQUREQnk5ZiCSE1kYpVyWIV9Ekm2e2qq1dKV7biqHVnrr7xMCpp9JNWViaI88dcIewVy/Kro6B/w+G01YmcrJTunsmtK3ZeXSUGzj6S6MrFKOf4aYScFSN5iSbVaurJWTi60gqT2aUxBREQCSgoiIhJQUhARkYCSgoiIBPJyoFkkH2Ri2uozZSsoW7aR0nYHZSSmsmUbg831khX7STUrK/Z132zdwT57FQTn1OYMpcqm4X7Zohv7f/ZhrcVTl+RlUtCUVMl1mZj6CqT1h7qqYjfXi14/fgpqqniir/tyy3b2b/jfPyG1mRRSTcMtW7aRvZt1UlKoprxMCpqSKrkukwV6StsdVGvFftI9HyKtjjCnrWrKbHZoTEFERAJKCiIiElBSEBGRgJKCiIgElBRERCSgpCAiIoG8nJIqkmuqs1At2TnzV2/m7IdmpLU2If4aVTm3KnFFF4XFL3pLR3VqM6ieQ3jyMilo8ZrkkuosVEt2TvTrdBaSJbtGsnOjr0kVY7JaEbEJoao1JKpTm0H1HMKTl0lBi9ckl1RnoVqyc6LH42sfpHuNVOemE2eqWhHJFr2lozoLzbQ4LRwaUxARkYCSgoiIBJQUREQkoKQgIiIBJQUREQkoKYiISEBJQUREAkoKIiISUFIQEZGAkoKIiASUFEREJKCkICIigZzZEM/M9gXuB7YC09396ZBDEhGpd8oxuBoAAAdaSURBVLLaUjCzR83sMzObG3e8n5ktNLMlZjai4vCPgefd/WLgjGzGJSIiiWW7pfAYcB/wRPSAmRUAfwH6AuXAe2b2AlAEfFTxsh1Zjksk51WncE91zq3JfWKvEbtld7QuQ6LCPPEFdKKvi79eosJBibb+TlRYaOs+LVhdfDbPlK0IzklWuCf2eGz88c+nKvaT7jXij1f2XKp7FR/SmJEDOqd8bXVkNSm4+5tm1jbu8DHAEndfCmBm44GBRBJEETCHFC0YMxsGDANo00aFN6Ruqk7hnuqcW937pCrYE1scJ1FhnvgCOrGvi4qeE3/N+D+YyYoSHdPxEOav3o9Jc1YG5yQr3BOfuOLvNWnOSsqWbUx4/6h0rpGsaFBVCwpF4yk+JDv1JsIYU2gFfBrzuBwoBe4B7jOz04DJyU529zHAGICSkhLPYpwioalO4Z7qnFvd+8SfF/t1bIshWWGe+FZJ7OsS3auywkGxz5e2O4hnL+mV8JxkraHY+6dT4Kg610jVEqtqK6203UFZaSVAOEnBEhxzd/8a+FlaF1A5ThGRrAhjSmo50DrmcRGwqioXcPfJ7j6sSZMmGQ1MRKS+CyMpvAd0NLN2ZrYXcA7wQghxiIhInGxPSR0HzAAON7NyM/u5u28HrgD+CSwAnnP3eVW87gAzG/PFF19kPmgRkXos27OPzk1yfAowpQbXnQxMLikpubi61xARkd1pmwsREQnkZVJQ95GISHbkZVLQ7CMRkeww9/xd/2Vm64BPqnFqM2B9hsPJBMVVNYqr6nI1NsVVNTWN61B3b57oibxOCtVlZu+7e0nYccRTXFWjuKouV2NTXFWTzbjysvtIRESyQ0lBREQC9TUpjAk7gCQUV9UorqrL1dgUV9VkLa56OaYgIiKJ1deWgoiIJKCkICIigXqXFJLUhw6VmbU2s2lmtsDM5pnZL8OOKZaZFZjZB2b2YtixRJnZAWb2vJl9XPFzS1yhpZaZ2dUV/4ZzzWycmTUMKY7d6qOb2UFm9pqZLa74fGCOxHV7xb/jh2Y2wcwOyIW4Yp77tZm5mTXLlbjM7MqKv2PzzOy2TN6zXiWFmPrQ/YFi4FwzKw43KgC2A9e4eyfgWODyHIkr6pdEdrTNJf8DvOLuRwBHkgPxmVkr4CqgxN27AAVEtoYPw2NAv7hjI4DX3b0j8HrF49r2GLvH9RrQxd27AYuA39Z2UCSOCzNrTaSe/IraDqjCY8TFZWYnESlh3M3dOwN3ZPKG9SopEFMf2t23AtH60KFy99XuPrvi6y+J/IGrWmHeLDGzIuA04K9hxxJlZo2BHwCPALj7VnffFG5UgQZAIzNrAOxDFQtIZYq7vwlsjDs8EHi84uvHgUG1GhSJ43L3Vyu21Af4F5HCW6HHVeEu4FoglBk5SeIaDox29+8qXvNZJu9Z35JCovrQOfHHN8rM2gI9gLJwIwncTeQ/xc6wA4nRHlgHjK3o1vqrme0bdlDuvpLIu7YVwGrgC3d/NdyodnGwu6+GyBsRoEXI8SRyEfBy2EEAmNkZwEp3/3fYscT5HvB9MyszszfM7OhMXry+JYWE9aFrPYokzGw/4O/Ar9x9cw7EczrwmbvPCjuWOA2AnsAD7t4D+JpwukJ2UdFHPxBoBxwC7GtmF4QbVf4ws98R6Up9Ogdi2Qf4HXBT2LEk0AA4kEhX82+A58ws0d+2aqlvSaHG9aGzxcz2JJIQnnb3f4QdT4XjgTPMbDmRrrY+ZvZUuCEBkX/HcnePtqaeJ5IkwnYysMzd17n7NuAfwHEhxxRrrZkVAlR8zmi3Q02Y2YXA6cD5nhuLpw4jktz/XfH7XwTMNrOWoUYVUQ78wyNmEmnFZ2wQvL4lhZysD12R5R8BFrj7nWHHE+Xuv3X3IndvS+RnNdXdQ3/n6+5rgE/N7PCKQz8E5ocYUtQK4Fgz26fi3/SH5MAAeIwXgAsrvr4QmBRiLAEz6wdcB5zh7t+EHQ+Au3/k7i3cvW3F73850LPidy9sE4E+AGb2PWAvMriTa71KCpmoD50lxwM/IfJOfE7Fx4/CDirHXQk8bWYfAt2BP4YcDxUtl+eB2cBHRP5/hbJNQqL66MBooK+ZLSYyo2Z0jsR1H7A/8FrF7/6DORJX6JLE9SjQvmKa6njgwky2rrTNhYiIBOpVS0FERFJTUhARkYCSgoiIBJQUREQkoKQgIiKBBmEHIJLPzKwpkc3lAFoCO4hswQEwARhScWwncEnMgjuRnKQpqSIZYmY3A1+5+x0VW3nfCfR29+8qtl3ey91zYgW9SDJqKYhkRyGwPmYny4ytOBXJJo0piGTHq0BrM1tkZveb2YlhBySSDiUFkSxw96+Ao4BhRMYYnjWzoaEGJZIGdR+JZIm77wCmA9PN7CMim9A9FmZMIpVRS0EkC8zscDPrGHOoO/BJWPGIpEstBZHs2A+4t6II/XZgCZGuJJGcpimpIiISUPeRiIgElBRERCSgpCAiIgElBRERCSgpiIhIQElBREQCSgoiIhL4/9NYVDE+gbElAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "(h, be) = np.histogram(trials['ts'], bins=np.arange(0, np.max(trials['ts'])+0.1, 0.1))\n",
    "plt.plot(0.5*(be[:-1]+be[1:]), h, drawstyle='steps-mid', label='background')\n",
    "plt.vlines(ts, 1, np.max(h), label=f'TS(TXS 0506+056)={ts:.3f}')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('TS')\n",
    "plt.ylabel('#trials per bin')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the TS value of the unblinded data for TXS is rather large and 10k trials are not enough to calculate a reliable estimate for the p-value. Hence, we will generate a few more trials. SkyLLH provides also a helper function to extend the trial data file we just created. It is called ``extend_trial_data_file``: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.analysis_utils import extend_trial_data_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function extend_trial_data_file in module skyllh.core.analysis_utils:\n",
      "\n",
      "extend_trial_data_file(ana, rss, n_trials, trial_data, mean_n_sig=0, mean_n_sig_null=0, mean_n_bkg_list=None, bkg_kwargs=None, sig_kwargs=None, pathfilename=None, **kwargs)\n",
      "    Appends to the trial data file `n_trials` generated trials for each\n",
      "    mean number of injected signal events up to `ns_max` for a given analysis.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    ana : Analysis\n",
      "        The Analysis instance to use for sensitivity estimation.\n",
      "    rss : RandomStateService\n",
      "        The RandomStateService instance to use for generating random\n",
      "        numbers.\n",
      "    n_trials : int\n",
      "        The number of trials the trial data file needs to be extended by.\n",
      "    trial_data : structured numpy ndarray\n",
      "        The structured numpy ndarray holding the trials.\n",
      "    mean_n_sig : ndarray of float | float | 2- or 3-element sequence of float\n",
      "        The array of mean number of injected signal events (MNOISEs) for which\n",
      "        to generate trials. If this argument is not a ndarray, an array of\n",
      "        MNOISEs is generated based on this argument.\n",
      "        If a single float is given, only this given MNOISEs are injected.\n",
      "        If a 2-element sequence of floats is given, it specifies the range of\n",
      "        MNOISEs with a step size of one.\n",
      "        If a 3-element sequence of floats is given, it specifies the range plus\n",
      "        the step size of the MNOISEs.\n",
      "    mean_n_sig_null : ndarray of float | float | 2- or 3-element sequence of\n",
      "                      float\n",
      "        The array of the fixed mean number of signal events (FMNOSEs) for the\n",
      "        null-hypothesis for which to generate trials. If this argument is not a\n",
      "        ndarray, an array of FMNOSEs is generated based on this argument.\n",
      "        If a single float is given, only this given FMNOSEs are used.\n",
      "        If a 2-element sequence of floats is given, it specifies the range of\n",
      "        FMNOSEs with a step size of one.\n",
      "        If a 3-element sequence of floats is given, it specifies the range plus\n",
      "        the step size of the FMNOSEs.\n",
      "    bkg_kwargs : dict | None\n",
      "        Additional keyword arguments for the `generate_events` method of the\n",
      "        background generation method class. An usual keyword argument is\n",
      "        `poisson`.\n",
      "    sig_kwargs : dict | None\n",
      "        Additional keyword arguments for the `generate_signal_events` method\n",
      "        of the `SignalGenerator` class. An usual keyword argument is\n",
      "        `poisson`.\n",
      "    pathfilename : string | None\n",
      "        Trial data file path including the filename.\n",
      "    \n",
      "    Additional keyword arguments\n",
      "    ----------------------------\n",
      "    Additional keyword arguments are passed-on to the ``create_trial_data_file``\n",
      "    function.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    trial_data :\n",
      "        Trial data file extended by the required number of trials for each\n",
      "        mean number of injected signal events..\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(extend_trial_data_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[==========================================================] 100% ELT 0h:29m:56s\n"
     ]
    }
   ],
   "source": [
    "tl = TimeLord()\n",
    "rss = RandomStateService(seed=2)\n",
    "trials = extend_trial_data_file(\n",
    "    ana=ana,\n",
    "    rss=rss,\n",
    "    n_trials=4e4,\n",
    "    trial_data=trials,\n",
    "    pathfilename='/home/mwolf/projects/publicdata_ps/txs_bkg_trails.npy',\n",
    "    ncpu=8,\n",
    "    tl=tl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TimeLord: Executed tasks:\n",
      "[Generating background events for data set 0.]   0.002 sec/iter (40000)\n",
      "[Generating background events for data set 1.]   0.003 sec/iter (40000)\n",
      "[Generating background events for data set 2.]   0.003 sec/iter (40000)\n",
      "[Generating background events for data set 3.]   0.005 sec/iter (40000)\n",
      "[Generating background events for data set 4.]   0.019 sec/iter (40000)\n",
      "[Generating pseudo data.                     ]   0.027 sec/iter (40000)\n",
      "[Initializing trial.                         ]   0.032 sec/iter (40000)\n",
      "[Create fitparams dictionary.                ] 1.1e-05 sec/iter (2375320)\n",
      "[Calc fit param dep data fields.             ] 3.3e-06 sec/iter (2375320)\n",
      "[Get sig prob.                               ] 2.0e-04 sec/iter (2375320)\n",
      "[Evaluating bkg log-spline.                  ] 2.8e-04 sec/iter (2375320)\n",
      "[Get bkg prob.                               ] 3.5e-04 sec/iter (2375320)\n",
      "[Calc PDF ratios.                            ] 6.8e-05 sec/iter (2375320)\n",
      "[Calc pdfratio values.                       ] 8.5e-04 sec/iter (2375320)\n",
      "[Calc pdfratio value product Ri              ] 3.9e-05 sec/iter (2375320)\n",
      "[Calc logLamds and grads                     ] 3.1e-04 sec/iter (2375320)\n",
      "[Evaluate llh-ratio function.                ]   0.005 sec/iter (475064)\n",
      "[Minimize -llhratio function.                ]   0.054 sec/iter (40000)\n",
      "[Maximizing LLH ratio function.              ]   0.054 sec/iter (40000)\n",
      "[Calculating test statistic.                 ] 3.7e-05 sec/iter (40000)\n"
     ]
    }
   ],
   "source": [
    "print(tl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The local p-value is defined as the fraction of background trials with TS value greater than the unblinded TS value of the source. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-log10(p_local) = 2.93\n"
     ]
    }
   ],
   "source": [
    "minus_log10_pval = -np.log10(len(trials[trials['ts'] > ts]) / len(trials))\n",
    "print(f'-log10(p_local) = {minus_log10_pval:.2f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXzU9bX/8dcxiOCCyqaRgARBSwANmBqp3opeUVAR6A9BsI9KbQUXtO1tq5YuiF2ktbW91VqFiksVxNqLiHKtbZHiQoNgc2WTpUAxECGAgBv7+f0x+X47mUySSchkZjLv5+PBI5nvzPc7J+M4Zz7rMXdHREQE4KhUByAiIulDSUFEREJKCiIiElJSEBGRkJKCiIiEWqQ6gCPRvn1779q1a6rDEBHJKEuXLt3u7h3i3ZfRSaFr164sWbIk1WGIiGQUM/tXTfep+0hEREJKCiIiElJSEBGRUEaPKYg0lQMHDlBWVsbevXtTHYpIwlq1akVeXh5HH310wucoKYgkoKysjBNOOIGuXbtiZqkOR6RO7s6OHTsoKysjPz8/4fPUfSSSgL1799KuXTslBMkYZka7du3q3brNyKRgZkPMbOru3btTHYpkESUEyTQNec9mZFJw97nuPu7EE09s0PmT565g1COLmFGyqZEjExHJbBmZFBpDyYadzCndnOowRBK2ceNGevfuHfe+1atXs3r16jqvsWDBAq666qrGDq1RjB07lueeey7VYWS9rBxonjSkFyu37El1GCIZxd1xd446Kmu/S2YF/dcVySAHDx7k+uuv5+yzz2bEiBF88skn3HPPPYwYMYIhQ4Ywbtw4gmqK69at49JLL+Wcc86hX79+/POf/6xyrbfeeou+ffuyfv16KioqGDhwIP369WP8+PGcfvrpbN++nY0bN9KzZ09uueUW+vXrx3vvvcfMmTPp06cPvXv35s477wyvd/zxx4e/P/fcc4wdOxaItABuv/12Pve5z9GtW7ewNeDuTJgwgYKCAq688kq2bduW5FdPEpGVLQWRIzF57opGb2kWnNaGSUN61fm41atX8+ijj3LBBRdwww038NBDDzFhwgRGjRoFwI9+9CNefPFFhgwZwnXXXcddd93F8OHD2bt3L4cPH+a9994D4M033+S2225jzpw5dOnShQkTJnDJJZfwne98h5dffpmpU6dWec7HHnuMhx56iC1btnDnnXeydOlSTj75ZC677DKef/55hg0bVmvc5eXlvP7667z77rtcffXVjBgxgtmzZ7N69WqWLVvG1q1bKSgo4IYbbjiCV1Eag1oKIhmkc+fOXHDBBQB88Ytf5PXXX+fVV19l5MiRDBkyhPnz57NixQo+/PBDNm/ezPDhw4HIIqZjjz0WgFWrVjFu3Djmzp1Lly5dAHj99de59tprARg0aBAnn3xy+Jynn346559/PhBpXQwYMIAOHTrQokULrrvuOhYuXFhn3MOGDeOoo46ioKCArVu3ArBw4UJGjx5NTk4Op512GpdcckkjvUpyJNRSEKmnRL7RJ0vsFEMz45ZbbmHWrFnk5uYyc+ZM9u7dG3YhxZObm8vevXv5xz/+wWmnnQZQ6+OPO+648PfaHhcdW+zc+GOOOSbuNTTNN/2opSCSQTZt2sSiRYsAmDlzJhdeeCEAJ598Mh9//HHYX9+mTRvy8vJ4/vnnAdi3bx+ffPIJACeddBIvvfQSEydOZMGCBQBceOGFPPvsswC88sorfPDBB3Gfv7i4mL/97W9s376dQ4cOMXPmTC666CIATjnlFFatWsXhw4eZPXt2nX/L5z//eZ555hkOHTpEeXk5r776agNfFWlMaimIZJCePXvyxBNPMH78eHr06MHNN9/MBx98wNVXX02nTp347Gc/Gz7297//PePHj+cHP/gBRx99NH/4wx/C+0455RTmzp3L4MGDmT59OpMmTWL06NHMmjWLiy66iNzcXE444QQ++uijKs+fm5vLvffey8UXX4y7c8UVVzB06FAApkyZwlVXXUXnzp3p3bt3tXNjDR8+nPnz59OnTx/OPPPMMLlIalltzcF0V1RU5A0tsjPqkci3rVnj+zdmSNJMrVq1ip49e6Y6jBoFaxTOOuusBp2/b98+cnJyaNGiBYsWLeLmm2+mtLS0MUOUFIn33jWzpe5eFO/xadVSMLPjgIXAJHd/MdXxiGSLTZs2MXLkSA4fPkzLli2ZNm1aqkOSFElqUjCz6cBVwDZ37x11fBDw30AO8Dt3n1J5153As8mMSUSq69GjB//4xz9SHYakgWQPND8ODIo+YGY5wG+AwUABMNrMCszsUmAlsDXJMYmISA2S2lJw94Vm1jXm8HnAOndfD2BmzwBDgeOB44gkik/NbJ67H469ppmNA8YB4RxrERFpHKkYU+gEvBd1uwwodvcJAGY2FtgeLyEAuPtUYCpEBpqTG6qISHZJRVKIt1ol/HB398ebLhQREYmWisVrZUDnqNt5wJb6XEBFdkREkiMVSeEtoIeZ5ZtZS+Ba4IX6XOBIi+yIZJodO3ZQWFhIYWEhp556Kp06dQpvT548mauuuoqrr76awsJCSkpKwvNGjBjB+vXrKS4uprCwkC5dutChQ4fw3GXLlnHGGWewdu1aAA4cOECfPn3Ca/z4xz+mV69enH322dWuHdi5cycDBw6kR48eDBw4MFwNvXHjRlq3bh0+10033RSes3TpUvr06UP37t25/fbbq2x98eyzz1JQUECvXr0YM2ZMwq/Rhg0bKC4upkePHowaNYr9+/cDkRoSJ554YhjHPffcE56za9cuRowYwWc+8xl69uwZrhb/1re+xfz58xN+7nfffZf+/ftzzDHH8POf/zw8vnfvXs477zzOOeccevXqxaRJk+p1fuDQoUP07du3Si2Mu+++u8r7YN68eQnHW6tgj/Rk/ANmAuXAASIthK9UHr8CWAP8E/huA647BJjavXt3b6iRD7/pIx9+s8HnS3ZZuXJlqkMITZo0ye+77z53d3/zzTf9/PPP93feecffffddr6io8M2bN7u7+/Lly33YsGFVzn3sscf81ltvrXJs1qxZPnDgQHd3/8lPfuLjxo2rcu29e/e6u1e5drRvf/vbfu+997q7+7333ut33HGHu7tv2LDBe/XqFfdv+OxnP+tvvvmmHz582AcNGuTz5s1zd/c1a9Z4YWGh79y5093dt27dWu3cxx57zCdNmlTt+DXXXOMzZ850d/fx48f7Qw895O7ur776ql955ZVx4/jSl77k06ZNc3f3ffv2+QcffODu7hs3bgxfk0Rs3brVFy9e7BMnTgz/27i7Hz582D/88EN3d9+/f7+fd955vmjRooTPD/ziF7/w0aNHV/k7ot8HtYn33gWWeA2fr8mefTS6huPzgAanNXefC8wtKiq6saHXEDkSAwYMaNTrBXsQ1Vd5eTnt27enZcuWALRv3z687+mnnw63oKjNyJEjmT59Oj/72c94+OGHw/UKwbWDzeyirx1tzpw5YfzXX389AwYM4Kc//WmtMe/Zs4f+/SO7CXzpS1/i+eefZ/DgwUybNo1bb7013KW1Y8eOdcYPkS+38+fPZ8aMGWEcd999NzfffHON5+zZs4eFCxfy+OOPA9CyZcvwdTz99NPZsWMH77//Pqeeemqdz9+xY0c6duzISy+9VOW4mYV1Jg4cOMCBAwfibgJY0/kAZWVlvPTSS3z3u9/l/vvvrzOWI6UN8UQy2GWXXcZ7773H5ZdfzuTJk/nb3/4W3vfGG29w7rnnJnSdX/3qV9x5551873vfo23btlWufeaZZ3LLLbdUuXa0rVu3kpubC0T2RooulrNhwwb69u3LRRddxGuvvQbA5s2bycvLCx+Tl5fH5s2R0rhr1qxhzZo1XHDBBZx//vm8/PLLCcW/Y8cOTjrpJFq0aFHtmgCLFi3inHPOYfDgwaxYsQKA9evX06FDB7785S/Tt29fvvrVr/Lxxx+H5/Tr14833ngDgG984xthN030vylTplCXQ4cOUVhYSMeOHRk4cCDFxcUJ/U2Br3/96/zsZz+LW/HuwQcf5Oyzz+aGG26ocRPD+kqrbS4SZWZDgCHdu3dPdSiSpRr6zb6xHX/88SxdupSnn36akpISRo0axZQpUxg7dizl5eV06NAhoeu8/PLL5Obmsnz58mrXfu2113j11VerXDsRubm5bNq0iXbt2rF06VKGDRvGihUr4m6/HXx7PnjwIGvXrmXBggWUlZXxH//xHyxfvpxDhw7xn//5n0BkDGP//v3hDrC///3v436bD67Zr18//vWvf3H88cczb948hg0bxtq1azl48CBvv/02DzzwAMXFxXzta19jypQp/PCHPwQi3963bInMgfnlL3+Z0N8cT05ODqWlpezatYvhw4ezfPnyGmttx3rxxRfp2LEj5557brX33M0338z3v/99zIzvf//7fPOb32T69OkNjjOQkS0F10CzSCgnJ4fi4mJuv/12HnzwQf74xz8C0Lp162p1DeLZsmULv/71r1m8eDHz5s3jnXfeqXLtAQMGMHny5CrXjnbKKadQXl4ORLqGgi6fY445hnbt2gFw7rnncsYZZ7BmzRry8vIoKysLzy8rKwvrOuTl5TF06FCOPvpo8vPzOeuss1i7di3t2rWjtLSU0tJS7rnnHm666abwdp8+fWjfvj27du3i4MGD1a7Zpk2bsAvniiuu4MCBA2zfvp28vDzy8vLCb+4jRozg7bffDuPau3cvrVu3Bo6spRA46aSTGDBgQMKtH4i09l544QW6du3Ktddey/z58/niF78Yvu45OTkcddRR3HjjjSxevDjh69YmI5OCiESsXr06nDkEUFpayumnnw5Ettlet25dndf4xje+wcSJE8nLy+P+++/n1ltvxd1rvXa0q6++mieeeAKAJ554IhzHqKio4NChQ0Ckq2bt2rV069Yt3Jb773//O+7Ok08+GZ4zbNiwsK7C9u3bWbNmDd26davzbzAzLr744rCeRHQc77//ftg6Wbx4MYcPH6Zdu3aceuqpdO7cOdxh9q9//SsFBQXhNdesWRN+o//lL38ZJqHof3fddVetcVVUVLBr1y4APv30U/7yl7/wmc98ps6/J3DvvfdSVlbGxo0beeaZZ7jkkkt46qmnAMJEDDB79uyEWx91UfeRSAb76KOPuO2229i2bRs5OTn06tUrrK985ZVXsmDBAi699NIaz//zn//Mpk2b+MpXvgLAkCFDmDZtGk8++SS9e/fmtttuY9euXbRo0YLu3btXqd0cuOuuuxg5ciSPPvooXbp0Ces2LFy4kB/84Ae0aNGCnJwcHn744XC84re//S1jx47l008/ZfDgwQwePBiAyy+/nFdeeYWCggJycnK47777wtZGXX76059y7bXX8r3vfY++ffuGf9Nzzz3Hb3/7W1q0aEHr1q155plnwq6lBx54gOuuu479+/fTrVs3HnvsMSAyKLxu3TqKiuLuLl3N+++/T1FREXv27OGoo47iV7/6FStXrqS8vJzrr7+eQ4cOcfjwYUaOHBlOK3344YcBuOmmm2o8v02bNjU+5x133EFpaSlmRteuXXnkkUcSirUuqqegegqSgEysp/Dpp59y8cUX88Ybb5CTk5Oq0DLS7Nmzefvtt8PxhUxW33oK6j4SaaZat27N5MmTq8zCkcQcPHiQb37zm6kOIyXUfSSSIHfPuELzl19+eapDyEjXXHNNqkNoFA3pCcrIloJmH0lTa9WqFTt27GjQ/2QiqeDu7Nixg1atWtXrvIxsKYg0tWAaZUVFRapDiev9998H4PDhuDvOS5Zq1apVlYWCiVBSEElAMG8+XQXbOaTLojrJXBnZfSQiIsmRkUlB9RRERJIjI5OCBppFRJIjI5OCiIgkh5KCiIiElBRERCSkpCAiIqGMTAqafSQikhwZmRQ0+0hEJDkyMimIiEhyZHVSKNmwkxklm1IdhohI2sjapDC0sBMAc0q117yISCBrk8KY4i4U57dNdRgiImkla5OCiIhUl/VJYWX5HkY9skhjCyIiZGg9hcYqxxmMK6ws3wNEupRERLJZRrYUGmudwpjiLswa35+C3DaNFJmISGbLyKQgIiLJoaRQSWsWRESUFACtWRARCSgpoDULIiIBJYUo6kISkWynpFAp6EKaOHuZ1i2ISNbKyHUKyRCsUZhTulnrFkQka6mlECV63YK6kkQkGykpxKHZSCKSrdImKZhZTzN72MyeM7ObUxlLMBtJ+yKJSLZJ6piCmU0HrgK2uXvvqOODgP8GcoDfufsUd18F3GRmRwHTkhlXIoLWQsmGnZRs2AlojEFEmr9ktxQeBwZFHzCzHOA3wGCgABhtZgWV910NvA78Nclx1SkYX/jJ8D5AZFaSWgwi0twlNSm4+0JgZ8zh84B17r7e3fcDzwBDKx//grt/Driupmua2TgzW2JmSyoqKpIVemhMcZcqiUHdSSLSnKViSmon4L2o22VAsZkNAL4AHAPMq+lkd58KTAUoKiry5IX5b5quKiLZIhUDzRbnmLv7Ane/3d3Hu/tvar2A2RAzm7p79+4khVidpquKSDZIRVIoAzpH3c4DttTnAo1VT6EhNF1VRJqzVCSFt4AeZpZvZi2Ba4EX6nOBVLQUAto8T0Sas2RPSZ0JDADam1kZMMndHzWzCcCfiExJne7uK+pzXXefC8wtKiq6sbFjTlSwhiEwtLCTxhlEJOMlNSm4++gajs+jlsHkdBd0IQWCtQxzSjcrOYhIRqszKZjZBcDdwOmVjzciA8PdkhtarTENAYZ07949Jc8/prhLlQ/+GSWbNDNJRJqFRMYUHgXuBy4EPgsUVf5MmVQONMcTPTNJW2OISCZLpPtot7v/b9IjaQZit8ZQd5KIZJpEksKrZnYf8D/AvuCgu7+dtKjqkOruo5oE3UpBd5L2TRKRTJNI91ExkS6jnwC/qPz382QGVZd06z6KFbtvktY0iEimqLOl4O4XN0UgzdGY4i7hAPSoRxapK0lE0l6NScHMvujuT5nZf8W7393vT15YzYfGGUQkk9TWUjiu8ucJTRFIfaTrmEI8seMMmrYqIunM3Jtko9GkKCoq8iVLlqQ6jHoZ9cgiVpbvoSC3jVoM0mgGDBgAwIIFC1Iah2QGM1vq7kXx7ktk8Vo3IlXSzgccWAR8w93XN2qUWUIV3UQknSUy+2gG8CyQC5wG/AGYmcygmrN4Fd202E1E0kUi6xTM3X8fdfupyg3tUiaTxhRqEl24R4PQIpIuahxTMLNgf+g7gF1EymY6MAo4xt1/2CQR1iITxxTiiV7sBlCc31bJQepFYwpSHw0dU1hKJAkEldLGR93nQMqTQnOhGUoiki5qTArunt+Ugci/k8OoRxZRsmEnQx54nWNb5gCq1yAiTSOp9RSkYYYWdqJkw06Wbd7NCa3+/Z9ISUFEki0V5TilDtElPwty21CQ2ybFEYlItqi1pWBmBuS5+3tNFE9CmsPso7oE6xmGFnYKB6GDaauapSQiyVJrUnB3N7PngXObKJ6EpEON5mSLre4WTFkNftfCNxFJhkS6j/5uZimttJbtoruTomlLbhFpbIkMNF8M3GRmG4GP+XeN5rOTGZhUF0xVDRKEtuQWkcaWSFIYnPQopE7BGEPs7+pKEpHGlEiRnX+Z2YVAD3d/zMw6AMcnPzSJFjvGEBybUbKJibOXMad0s5KCiByxRHZJnUSkHOdZwGPA0cBTwAXJDU0SEVvdDbTQTUQaLpHuo+FAX+BtAHffYmZpV3gnm8XrTgoGoZUgRKQ+EkkK+yunpjqAmR1X1wnJlg3rFOojumsp2D8JIgPR2n1VROojkSmpz5rZI8BJZnYj8BdgWnLDqp27z3X3cSeeeGIqw0hLQb2GWeP7hyuhSzbsVN0GEUlIIgPNPzezgcAe4EzgB+7+56RHJo2qOL+tdl8VkToluvfRMuA1YGHl75JBivPbhi2H6O0yRERi1ZkUzOyrwGLgC8AIIiucb0h2YNL4ggHpOaWbmVGySd1JIlJNIgPN3wb6uvsOADNrB7wJTE9mYNL4gumrQJVKb+pOEpFAIkmhDPgw6vaHQFrtmirxRe+0KiKSiESSwmagxMzmECnDORRYbGb/BeDu9ycxPjkC8VZBw7/3UILIzCQtehORQCJJ4Z+V/wJzKn9qAVsGil3oBlWThJKCSHZLZErq5KYIRJpGdOshqAWtym4iElCNZglF75/0yf5DHNsyR11KIlkmrZKCmQ0DrgQ6Ar9x91dSHFKzVtNA9MryPXy492B4W0lBJHskPSmY2XTgKmCbu/eOOj4I+G8gB/idu09x9+eB583sZODngJJCEsUORMd2K4lI9klk6+yfAT8CPgVeBs4Bvu7uTyX4HI8DDwJPRl0zB/gNMJDIlNe3zOwFd19Z+ZDvVd4vKRC0HEo27IzbpRT9OLUiRJqXRLa5uMzd9xD5tl9GZP+jbyf6BO6+EIj92nkesM7d17v7fuAZYKhF/BT4X3d/O971zGycmS0xsyUVFRWJhiH1EGyq95PhfcJB6JXle1i2eXc4U2ll+R7ViBZphhJJCkdX/rwCmOnujdGv0ImqC+DKKo/dBlwKjDCzm+Kd6O5T3b3I3Ys6dOjQCKFITeLtuFqQ2ya8HbQitFWGSPORyJjCXDN7l0j30S2V5Tj3HuHzWpxj7u6/Bn5d58mqp5ByQReTdl4VaV7qbCm4+11Af6DI3Q8AHxNZ1XwkyoDOUbfzgC2Jnqx6Ck1vaGEnivPbhskgaEVojYNI81JjS8HMvhDnWPTN/zmC530L6GFm+US20bgWGHME15Mkq2nLDBFpXmrrPhpSy31OgknBzGYCA4D2ZlYGTHL3R81sAvAnIlNSp7v7isRCVveRiEiy1JgU3P3LjfEE7j66huPzgHkNvOZcYG5RUdGNRxKbNI6gcI9aEiKZL6HFa2Z2JdALaBUcc/d7khVUAvGopZAmhhZ2omTDTuaUbmZMcRdmlGxiTunmhNcwBI8PrqXEIpJaiSxeexg4FrgY+B2R6muLkxxXrdRSSB9B4Z6gtRBbvCc6SQDVEsac0s3apVUkjSTSUvicu59tZu+4+2Qz+wVHNsgszUzQWpg4exkntIq8peIlieA4VP3w1wwmkfSRSFL4tPLnJ2Z2GrADyE9eSHVT91F6if7WD5HtMJZt3l0lScQzo2QTJRt2UpzftkniFJG6JZIUXjSzk4D7gLeJzDz6XVKjqoO6j9JP7JTV6LGCmjbXix5L0JYZIukhkSI7P6z89Y9m9iLQyt13JzcsyXSxxXxqUpzfNhyXEJHUq23x2iXuPr+mRWzurnEFEZFmpraWwkXAfOIvYkt48VoyaEwhswUb6a0s36NBZpE0U9vitUlmdhSRbayfbcKY6qQxhcwWveNqbNU3rVsQSa1axxTc/XDldhRplRQkcxXnt2XW+P413h+7biE4BkoSIk0hkdlHfzazbwGziOyQCkAj1VUQCQWzlIIpqiUbdobHgqmtSgoiyZVIkZ0bgFuBhcDSyn9LkhmUND8ry/dUawFEi+1GihZdAS4wo2STCvyIJEEiSaGnu+dH/wMKkh1YbcxsiJlN3b1bM2MzwdDCThTktok7hhAYU9wl7iK2YMpqrGCltKayijSuRLqP3gT6JXCsyWigObMkWoshSBhazCaSOrWtUziVSN3k1mbWl3+X0GxDZIM8kUYVnTyUFERSo7aWwuXAWCKlMn/Bv5PCHmBicsMSqa5kw85wdXRNW2eIyJGpbZ3CE8ATZvb/3P2PTRiTSI1iB6uDhXCgKasijSGRMYWewS9mdoy770tiPAnRiubmL3p8IVowCyloKQS3o5NFfYr8iEhVNc4+MrM7zKw/kaI6gZp3NmtC7j7X3cedeOKJqQ5FkmRMcRdmje8ffrAPLexEcX7bKkkiWAg3a3z/MDloVpLIkamtpbAauAboZmavAauAdmZ2lruvbpLoRCo1dBC6vuVBRbJdbesUPiAyoLwOGAD8uvL4XWb2ZpLjEjkisZXf1HIQSUxtLYVBwCTgDOB+4P+Aj939y00RmMiRUiIQqb/aZh9NBDCz/wOeAvoCHczsdeADd4+3pbZIysVbGR09nVVdSSI1S2T20Z/c/S3gLTO72d0vNLP2yQ5MpL6CGUjxZiTF3lZSEIkvkXKcd0TdHFt5bHuyAhJpiOhZSbHTWGOnr4pIzRJpKYTc/f+SFYjIkYi3v1J0jWitgBZJTCK7pKYd7ZIqIpIcGZkUtHgtu8VbyNaYjxfJZvXqPhJJB4luxR3v8UG5z/rORNIiOMkWGdlSEGmooOAPRGYjRa9lqK2amxbBSbZQS0GySnSrIRiAHvXIorCwT1AXOvjwV8tAso2SgmS92JlJweK3YF2DkoJkEyUFkSjBzqtAOO4gkk2UFCRrBbORSjbsrLYauibaLkOaOyUFyVrB+EIwswiqr4aOR9tlSHOmpCBZr75TXKOrv6nVIM1N2iQFM+sGfBc40d1H1PV4kXSgVoM0N0ldp2Bm081sm5ktjzk+yMxWm9k6M7sLwN3Xu/tXkhmPSGMryG1T5zhEbesfRNJNshevPU6kWE/IzHKA3wCDgQJgtJkVJDkOkZTRwjfJJElNCu6+EIjdnvI8YF1ly2A/8AwwNJlxiDRUsCVGbG2GWMH4gloDkulSsc1FJ+C9qNtlQCcza2dmDwN9zew7NZ1sZuPMbImZLamoqEh2rJLForfESKQWg1oD0hykYqDZ4hxzd98B3FTXye4+FZgKUFRU5I0cm0godlaS6jJINkhFS6EM6Bx1Ow/YUp8LqJ6CpIK24JZskIqWwltADzPLBzYD1wJj6nMBd58LzC0qKroxCfGJxBW7BXc8wRjEJ/sPcWzLHKD63koi6SypScHMZgIDgPZmVgZMcvdHzWwC8CcgB5ju7ivqed0hwJDu3bs3dsgiCQlaC8HuqoGC3DasLN/Dh3sPckKrtFkGJJKwZM8+Gu3uue5+tLvnufujlcfnufuZ7n6Gu/+4AddV5TVJqTHFXZg1vn+VMYdgM73owelEBqhF0om+yoikqXjV3mL3adIqamlsGZkU1H0k2SBY9AZUKycaUFKQxpaR5TjVfSTZTN1SkkwZ2VIQyXQN6UziTrAAAAgLSURBVAaaUbKJkg07w8pw0dcpLxjF8dtXJS1eyR4Z2VLQOgXJFMHahlhBN9DK8j0Jr4KOV/MhuM7+YzvyUfuejRO0ZLWMTArqPpJMEcxSipcYGtINVJzftlqroiC3DS0/2XZEcYoEMjIpiIhIcmhMQaSJBKuda9smI3oaalOIN+01kfuk+crIpKApqZJpgg/5urbgjp6G2hTiTXtN5D5pvjKy+0hjCpJpgrEFTSWVdJeRSUFERJIjI7uPRDJZ0CUTzEgKxhqC36N/BuKtUYi9jkhjyMikoDEFSWfRO6jGuy96zCD2MQW5beJuux27RiH2OiKNJSOTguopSDqLrdgWe1/sAG5tg7hBCwKqrlGIvY5IY9GYgoiIhJQUREQkpKQgIiKhjBxTEMkmwUykutY47GvThRklm4DIwHT0gHXwe0NWJ0fPjgoGulXop/nKyKSg2UeSTmqbbXSkj49+TOzjo6+z7q8z2demS/hhHQxAB3WiP9x7MDyvvh/iQTKKniarQj/NV0YmBc0+knRS1wyiI3l8XTOZgvumbnuncuvsqmsWgg/0hs5SCupOQ9WZUFqZ3XxpTEFEREJKCiIiElJSEBGRkJKCiIiElBRERCSUkbOPRCQ9lGzYmZQ1DEHVt0D0tRtjbYSqytUsI5OC1imI1K04vy1DCztV+XBNhthtvhtjDcOc0s2sLN9DQW6bKtdrrEpwqipXs4zsPlLlNZHaBesLmuIDryC3TZV1C7G3j+S6qlbX9DIyKYiISHIoKYiISEhJQUREQkoKIiISUlIQEZGQkoKIiISUFEREJKSkICIiISUFEREJpc02F2Z2HPAQsB9Y4O5PpzgkEZGsk9SWgplNN7NtZrY85vggM1ttZuvM7K7Kw18AnnP3G4GrkxmXiIjEl+yWwuPAg8CTwQEzywF+AwwEyoC3zOwFIA9YVvmwQ0mOSyQrrSzfU6XWcuxmdnUJNpErzm9b530lG3Yyo2QTY4q71LrraXC7vvs0xV6zIdcJXo/6nJfqHVaD5y84rQ2ThvRq9OsntaXg7guB2Irh5wHr3H29u+8HngGGEkkQeXXFZWbjzGyJmS2pqKhIRtgiGWfBggUUnlMY976hhZ0ozm9Ln04nVttcrqGbzQUf6nXdDj60g11PIfJBPKd0c3gsuF1f0deMvm59BLuw1ue8YIfVZO8+W9fzJ0sqxhQ6Ae9F3S4DioFfAw+a2ZXA3JpOdvepwFSAoqIiT2KcIs3CmOIutX6jjW45JKI4vy1jirtU+VAMbkdvRx37oRnsehr9fEe6A2pwTWjY3xEbT6Yozm+blFYCpCYpWJxj7u4fA19O6AKqpyAikhSpmJJaBnSOup0HbKnPBVRPQUQkOVKRFN4CephZvpm1BK4FXqjPBcxsiJlN3b17d1ICFBHJVsmekjoTWAScZWZlZvYVdz8ITAD+BKwCnnX3FfW5rloKIiLJkdQxBXcfXcPxecC8ZD63iIjUX0Zuc6HuIxGR5MjIpKDuIxGR5MjIpCAiIslh7pm7/svMKoB/NfD09sD2RgynOdBrUp1ek+r0mlSXaa/J6e7eId4dGZ0UjoSZLXH3olTHkU70mlSn16Q6vSbVNafXRN1HIiISUlIQEZFQNieFqakOIA3pNalOr0l1ek2qazavSdaOKYiISHXZ3FIQEZEYSgoiIhLKuqRQQ33orGZmG81smZmVmtmSVMeTKvFqiptZWzP7s5mtrfx5cipjbGo1vCZ3m9nmyvdLqZldkcoYm5qZdTazV81slZmtMLOvVR5vFu+VrEoKUfWhBwMFwGgzK0htVGnjYncvbC5zrRvocWBQzLG7gL+6ew/gr5W3s8njVH9NAH5Z+X4prNzgMpscBL7p7j2B84FbKz9HmsV7JauSAjXXhxapqab4UOCJyt+fAIY1aVApVsNrktXcvdzd3678/UMiJQA60UzeK9mWFOLVh+5Uw2OziQOvmNlSMxuX6mDSzCnuXg6RDwOgY4rjSRcTzOydyu6ljOwmaQxm1hXoC5TQTN4r2ZYU4taHbvIo0s8F7t6PSLfarWb2+VQHJGntt8AZQCFQDvwiteGkhpkdD/wR+Lq770l1PI0l25LCEdeHbo7cfUvlz23AbCLdbBKx1cxyASp/bktxPCnn7lvd/ZC7HwamkYXvFzM7mkhCeNrd/6fycLN4r2RbUjji+tDNjZkdZ2YnBL8DlwHLaz8rq7wAXF/5+/XAnBTGkhaCD75Kw8my94uZGfAosMrd74+6q1m8V7JuRXPl9LlfATnAdHf/cYpDSikz60akdQCR8qwzsvU1qawpPoDINshbgUnA88CzQBdgE3CNu2fNwGsNr8kAIl1HDmwExgd96dnAzC4EXgOWAYcrD08kMq6Q8e+VrEsKIiJSs2zrPhIRkVooKYiISEhJQUREQkoKIiISUlIQEZFQi1QHIJLJzKwdkc3PAE4FDgEVlbdnAyMrjx0mMnWzpMmDFKkHTUkVaSRmdjfwkbv/3Mz6A/cDA9x9n5m1B1oGq8dF0pVaCiLJkQtsd/d9AO6+PcXxiCREYwoiyfEK0NnM1pjZQ2Z2UaoDEkmEkoJIErj7R8C5wDgiYwyzzGxsSoMSSYC6j0SSxN0PAQuABWa2jMgmaY+nMiaRuqilIJIEZnaWmfWIOlQI/CtV8YgkSi0FkeQ4HnjAzE4iUtN3HZGuJJG0pimpIiISUveRiIiElBRERCSkpCAiIiElBRERCSkpiIhISElBRERCSgoiIhL6/ymxUAt6QB7JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "(h, be) = np.histogram(trials['ts'], bins=np.arange(0, np.max(trials['ts'])+0.1, 0.1))\n",
    "plt.plot(0.5*(be[:-1]+be[1:]), h, drawstyle='steps-mid', label='background')\n",
    "plt.vlines(ts, 1, np.max(h), label=f'TS(TXS 0506+056)={ts:.3f}')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('TS')\n",
    "plt.ylabel('#trials per bin')\n",
    "plt.legend()\n",
    "pass"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}