{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Working with the public 10-year IceCube point-source data\n",
    "=="
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial shows how to use the IceCube public 10-year point-source data with SkyLLH."
   ]
  },
  {
   "attachments": {},
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "**Disclaimer**\n",
    "\n",
    "    The released 10-year IceCube point-source data can reproduce the published results only within\n",
    "    a certain amount of uncertainty due to the limited instrument response function binning \n",
    "    provided in the data release. The IceCube collaboration is able to reproduce the published \n",
    "    results using detailed direct simulation data, as done for the publication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import scipy.stats"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating a configuration\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we have to create a configuration instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.config import Config\n",
    "cfg = Config()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting the datasets\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now 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"
   ]
  },
  {
   "attachments": {},
   "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": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsc = create_dataset_collection(\n",
    "    cfg=cfg, \n",
    "    base_path='/home/mwolf/projects/publicdata_ps/')"
   ]
  },
  {
   "attachments": {},
   "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": 6,
   "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": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dsc.dataset_names"
   ]
  },
  {
   "attachments": {},
   "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 access operator ``[dataset1, dataset2, ...]`` of the ``dsc`` instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets = dsc['IC40', 'IC59', 'IC79', 'IC86_I', 'IC86_II-VII']"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting the analysis\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.analyses.i3.publicdata_ps.time_integrated_ps import create_analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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(cfg, datasets, source, refplflux_Phi0=1, refplflux_E0=1000.0, refplflux_gamma=2.0, ns_seed=100.0, ns_min=0.0, ns_max=1000.0, gamma_seed=3.0, gamma_min=1.0, gamma_max=5.0, kde_smoothing=False, minimizer_impl='LBFGS', cut_sindec=None, spl_smooth=None, cap_ratio=False, compress_data=False, keep_data_fields=None, evt_sel_delta_angle_deg=10, construct_sig_generator=True, tl=None, ppbar=None, logger_name=None)\n",
      "    Creates the Analysis instance for this particular analysis.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    cfg : instance of Config\n",
      "        The instance of Config holding the local configuration.\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",
      "    ns_min : float\n",
      "        Lower bound for ns fit.\n",
      "    ns_max : float\n",
      "        Upper bound for 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",
      "    gamma_min : float\n",
      "        Lower bound for gamma fit.\n",
      "    gamma_max : float\n",
      "        Upper bound for gamma fit.\n",
      "    kde_smoothing : bool\n",
      "        Apply a KDE-based smoothing to the data-driven background pdf.\n",
      "        Default: False.\n",
      "    minimizer_impl : str\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",
      "    cut_sindec : list of float | None\n",
      "        sin(dec) values at which the energy cut in the southern sky should\n",
      "        start. If None, np.sin(np.radians([-2, 0, -3, 0, 0])) is used.\n",
      "    spl_smooth : list of float\n",
      "        Smoothing parameters for the 1D spline for the energy cut. If None,\n",
      "        [0., 0.005, 0.05, 0.2, 0.3] is used.\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",
      "    evt_sel_delta_angle_deg : float\n",
      "        The delta angle in degrees for the event selection optimization methods.\n",
      "    construct_sig_generator : bool\n",
      "        Flag if the signal generator should be constructed (``True``) or not\n",
      "        (``False``).\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",
      "    logger_name : str | None\n",
      "        The name of the logger to be used. If set to ``None``, ``__name__`` will\n",
      "        be used.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    ana : instance of SingleSourceMultiDatasetLLHRatioAnalysis\n",
      "        The Analysis instance for this analysis.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(create_analysis)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As source we use TXS 0506+056."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.source_model import PointLikeSource"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = PointLikeSource(ra=np.deg2rad(77.35), dec=np.deg2rad(5.7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/5 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 43/43 [00:04<00:00,  9.11it/s]\n",
      "100%|██████████| 43/43 [00:04<00:00,  9.90it/s]\n",
      "100%|██████████| 43/43 [00:04<00:00,  9.70it/s]\n",
      "100%|██████████| 43/43 [00:04<00:00,  9.85it/s]\n",
      "100%|██████████| 43/43 [00:04<00:00,  8.89it/s]\n",
      "100%|██████████| 5/5 [01:17<00:00, 15.48s/it]\n",
      "100%|██████████| 220/220 [00:00<00:00, 2409.32it/s]\n"
     ]
    }
   ],
   "source": [
    "ana = create_analysis(cfg=cfg, datasets=datasets, source=source)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initializing a trial\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "events_list = [ data.exp for data in ana.data_list ]\n",
    "ana.initialize_trial(events_list)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maximizing the log-likelihood ratio function\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "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": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.random import RandomStateService\n",
    "rss = RandomStateService(seed=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "(log_lambda_max, fitparam_values, status) = ana.llhratio.maximize(rss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "log_lambda_max = 6.5725295599619535\n",
      "fitparam_values = [14.58017285  2.16856498]\n",
      "status = {'grad': array([ 2.22650525e-06, -7.55260353e-05]), 'task': 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH', 'funcalls': 27, 'nit': 22, 'warnflag': 0, 'skyllh_minimizer_n_reps': 0, 'n_llhratio_func_calls': 27}\n"
     ]
    }
   ],
   "source": [
    "print(f'log_lambda_max = {log_lambda_max}')\n",
    "print(f'fitparam_values = {fitparam_values}')\n",
    "print(f'status = {status}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the test-statistic\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "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": 17,
   "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}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unblinding the data"
   ]
  },
  {
   "attachments": {},
   "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": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.random import RandomStateService\n",
    "rss = RandomStateService(seed=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method unblind in module skyllh.core.analysis:\n",
      "\n",
      "unblind(rss, tl=None) method of skyllh.core.analysis.SingleSourceMultiDatasetLLHRatioAnalysis instance\n",
      "    Evaluates the unscrambled data, i.e. unblinds the data.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    rss : instance of RandomStateService\n",
      "        The instance of RandomStateService that should be used draw random\n",
      "        numbers from.\n",
      "    tl : instance of TimeLord | None\n",
      "        The optional instance of TimeLord that should be used to time the\n",
      "        maximization of the LLH ratio function.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    TS : float\n",
      "        The test-statistic value.\n",
      "    global_params_dict : dict\n",
      "        The dictionary holding the global parameter names and their\n",
      "        best fit values. It includes fixed and floating parameters.\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)"
   ]
  },
  {
   "attachments": {},
   "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": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "(ts, x, status) = ana.unblind(rss=rss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the corresponding flux normalization "
   ]
  },
  {
   "attachments": {},
   "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": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flux scaling factor = 1.422e-15\n"
     ]
    }
   ],
   "source": [
    "scaling_factor = ana.calculate_fluxmodel_scaling_factor(x['ns'], [x['ns'], x['gamma']])\n",
    "print(f'Flux scaling factor = {scaling_factor:.3e}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, our result corresponds to a power-law flux of:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.422e-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)')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluating the log-likelihood ratio function\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "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": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method evaluate in module skyllh.core.llhratio:\n",
      "\n",
      "evaluate(fitparam_values, src_params_recarray=None, 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 : instance of numpy ndarray\n",
      "        The (N_fitparams,)-shaped numpy 1D ndarray holding the current\n",
      "        values of the global fit parameters.\n",
      "    src_params_recarray : instance of numpy record ndarray | None\n",
      "        The numpy record ndarray of length N_sources holding the parameter\n",
      "        names and values of all sources.\n",
      "        See the documentation of the\n",
      "        :meth:`skyllh.core.parameters.ParameterModelMapper.create_src_params_recarray`\n",
      "        method for more information about this array.\n",
      "        It case it is ``None``, it will be created automatically from the\n",
      "        ``fitparam_values`` argument using the\n",
      "        :class:`~skyllh.core.parameters.ParameterModelMapper` instance.\n",
      "    tl : instance of TimeLord | None\n",
      "        The optional instance of TimeLord that should be used for timing\n",
      "        measurements.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    log_lambda : float\n",
      "        The calculated log-lambda value of the composite\n",
      "        log-likelihood-ratio function.\n",
      "    grads : instance of numpy ndarray\n",
      "        The (N_fitparams,)-shaped 1D ndarray holding the gradient value of\n",
      "        the composite log-likelihood-ratio function for each global fit\n",
      "        parameter.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(ana.llhratio.evaluate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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}')"
   ]
  },
  {
   "attachments": {},
   "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 68%, 90%, and 95% quantile assuming Wilks-theorem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "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",
    "delta_ts = 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",
    "\n",
    "        delta_ts[ns_i, gamma_i] = (\n",
    "            ana.calculate_test_statistic(llhratio_value, [14.58, 2.17]) -\n",
    "            ana.calculate_test_statistic(ana.llhratio.evaluate([ns, gamma])[0], [ns, gamma])\n",
    "        )\n",
    "\n",
    "# Determine the best fit ns and gamma values from the scan.\n",
    "index_max = np.argmin(delta_ts)\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": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine the delta lambda value for the 95% quantile assuming a chi-sqaure\n",
    "# distribution with 2 degrees of freedom (i.e. assuming Wilks theorem).\n",
    "chi2_68_quantile = scipy.stats.chi2.ppf(0.68, df=2)\n",
    "chi2_90_quantile = scipy.stats.chi2.ppf(0.90, df=2)\n",
    "chi2_95_quantile = scipy.stats.chi2.ppf(0.95, df=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.5, 4.0)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAF5CAYAAACof5IgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB3GUlEQVR4nO2dd5hcVd3HP2fa9t5Ldje9N0ioopFiV6yggmJFURCkBhAQpAQQsCEa0deooIKoWBBBIIgikN7LJluT7b3vtPP+MRNI2TJz7syduTvnk2ef2Zncc8+Zu3fmd77nV46QUqLRaDQajcZ8bLEegEaj0Wg0iYo2whqNRqPRxAhthDUajUajiRHaCGs0Go1GEyO0EdZoNBqNJkZoI6zRaDQaTYwwxQgLIb4phNglhNgphPitECJZCJErhHheCFEdfMwxYywajUaj0cQLUTfCQogy4BvACinlIsAOfBJYDbwgpZwNvBB8rtFoNBpNwmDWcrQDSBFCOIBUoAk4H1gX/P91wIdNGotGo9FoNHFB1I2wlPIw8F2gAWgGeqWUzwFFUsrm4DHNQGG0x6LRaDQaTTzhiHYHQV/v+cB0oAd4UghxcRjtLwUuBUhJSTm5qqoq7DE49uwJuw2Aa4ZSMwCGs9Xa7fYqdtiv2E5MU2wI+NPU2rkVb7thtWYADCi2U/17qGLkE5ms2C5FsZ2Rsaq2tfnU2gmPYjvF/gBQbetXbGegBLF0q7Wrbe+QUhaodzw2Z5xxhuzp6VFuv2fPnn9KKd8TuRFFj6gbYeBcoFZK2Q4ghPgjcAbQKoQokVI2CyFKgLaxGksp1wJrARYsWCB//etfhz2AghUrlAZecY9SMwC2na/Wblm3YofrFdvZr1FsCAydqtbusOKixxa1ZgC8ptiu00CfKuQZaDtbsd1SxXZGvnqzFduldqm1S25Ra2fvUWsHwKDJ7RQnGgDeOrV2n/5JvXqn49PT2cmvf/xj5fYrzj03P4LDiSpmGOEG4DQhRCoBLXMOsJHAnXYJsCb4+PSkZxodhbq6sAfgDLtFEFVlAXiFYkPVybNqO7tdsSEgFW8f1bEaUcKqbY30qYkfVO9VQx471c+WajsDRlgormpFC68XDChhKxF1IyylfF0I8QdgM4HFvS0ElG068IQQ4osEDPUnoj0WjUaj0VgAny/qRlgI8U3gSwTW8XcAnycQOPx7oAqoAy6QUqquT4aEGUoYKeVtwG3HvTxKQBWHjtutpIRdYbcIYkAJj6hOZlXdQartjMz0VdWF6liHFNsZQtUpHAO/t9lqPxbuUuV7TvWrzsBKkTJmK2jAlqre1oIclTq7QEo5LIR4gkDq7AICqbNrhBCrCaTO3hDNsZhihDUajUajCRkTlDBvpc56eCt19kZgVfD/1xGIttFG+E0s5BNWVsKKQYrK7ZINzJ5V1UUsfMKmq2iTFXQsUF59MdBWWdHGArOrAhvoL958wsaNcL4QYuNRz9cGg3yBQOqsEOJI6uww8JyU8jkhxDGps0KIqKfOWumO1mg0Gk0iYDwwq0NKOW5ajNHU2UhiLSOs6BOOhRJWTUs13yeso6MnxuRE4SEDH0nV9xgTX7vJKEdHG/EJq7bV++qYgKHU2UhiLSOs0Wg0mqlP9H3CkUudNYi1jLCiT1h5Lmvg6jSqNlT17SrLdgtFRxtCVdFaKDrabEVrpehoSxGD6OjYfCjHJ8pGOJ5SZxPhjtZoNBqNlTAhOjpiqbMGsZYRNtkn7FN2JkOvakPVCalydKNyFjX4TfZdWcp3aXbRaQMkQlUw0/OLAZvZOcZG+jOyrBEFEqhilo4A0Gg0Go0mRlhLCfv92AbCjzsWioXmjShh033CdtVq+gbmYaqT5xHFdoYUm0V8wkYU9LDJfmgr+YRVFa0hH7TqZysGtaPjDXOKdcQF1jLCGo1Go5n6aCMcnwgU55bFav1JAy4W6/iEDbxJs/N9LaWENVHB7DiEmKD6HmOx41OU8PtBYdXTiljKCGs0Go0mMUiE6RZY0AgrzdcUlXBMfMKq6lJ5FxQjtaMV21kqylmVGChvVZ+w6t8jFrWjVVH17RryCVspOjrO8oQTCMsZYY1Go9FMfUSsB2ASljLCAnOVsN9ACm2DakPVAEdVn7CRPEizfcKGFLRVoqMTBNV7x+x01pjs2qSjowXaCGs0Go1GEzO0TzhOMTM62ohPuE+1oepMXzVP2IjPKyGio83uz8A4VVcKVK+rar43QIaBtirEpOa02TWgp1DFrATCckZYo9FoNFMfvRwdh5jtEzaUJ6w6sTTbJ2xEIZgeHR2LnF0LVcwyGyMrE2b7hGNRMUu1T+V12CmUJ4w2whqNRqPRxATlwkwWxHJGWGm+lq3W12gs0u5UZ/qqStgXA59wTEiAilkx8bWbTEKksyaK+ZmYRFHC+q+t0Wg0Gk2MsJQStqG4+22yWn8jRpSw6m5IZlfM8looOtqQKrWKorXKODGWt50IFbOU28YiOjr+lhgSRQlbyghrNBqNZuqjfcJxigCUUncV3+WIkbvA9OjoTLV2sYiOjokStlKfqiiOdchCXwNWqphlem7y1MoT1kpYo9FoNJoYoY1wHKKshBV9wl4jSlhV7SlPSBUHG4va0THJE06AillmE4s8YVViUjFLFV07OpGw0p2p0Wg0mgRB+4TjEIHJ0dFG7gJVf6lykKLi7NmveHGMEJO81ETYRUlxrKr7EMcCKylo1VUm5UpbU6dilnJ1RAtioU+fRqPRaBIBbYTjFLN9wjHJE1Zt51JVwomSJ2yVPi3k97aST9ifCIubUytPOFGI+p0phJgrhNh61E+fEOIqIUSuEOJ5IUR18DEn2mPRaDQajTWwGfixElFXwlLKfcAyACGEHTgM/AlYDbwgpVwjhFgdfH7DROcyPTraSIy8aqCi8oRU8daLRZ6wMhZSiVYiEWpHqyrvWFTMUu5zVLFd/JFIy9FmTxrOAQ5KKeuB84F1wdfXAR82eSwajUajiVPsBn6shNk+4U8Cvw3+XiSlbAaQUjYLIQrHaiCEuBS4FKDCBiIv/E79SiHV0BOLilmqPuFkxTdpZKY/otjOSM1hZaySt2shtW+l2tGq6IpZMSGRylaa9j6FEC7gQ8CT4bSTUq6VUq6QUq4oSJS/ikaj0WgSAjOnau8FNkspW4PPW4UQJUEVXAK0TXoGJ1AcfseqSrhRrVkAq/iEjUSNmr5/rfbrRoVYrEyorqKoCrb4EnoTo6qgjah2W/xdIKstK6tiprb8FG8tRQP8Bbgk+PslwNMmjkWj0Wg0ccqRwKxo+YTjKWvHFCUshEgFzgO+ctTLa4AnhBBfBBqAT0x6IgdqSlhxqtGr1iyA8sRStYKV4rzRyATY9BrQ5vtL7YoXSP2yWsgnbKWoalUsVTFLrVm8Es23E8msHaOYYoSllENA3nGvdRKIltZoNBqNJla8mbUjhDgfWBV8fR2wnqlghDUajUajCZUI5AnnCyE2HvV8rZRy7TjHhp21E0msZYQTITDLXqDaUK2ZkeVo0wOszA/MSlWMlOu3UlhJLJaVzd7q0yopUUYwlBIVX0GPETDCHVLKFZP281bWzo3GulNninkRNBqNRjMVMKls5ZhZOwAhZ+0YJCGUsE+p1mUgWkwZ1Rm7SFNsqJqipNgdWCowRzXAyuz+fDowa2LMVrSWKltpACPpTdFAgDCyy2ro9+Z4WTtrMClrRythjUaj0SQcR2Xt/PGol9cA5wkhqoP/tyba44iz6c8kmKyEY+ITVlbCio7vBElRUvXtpiGV2sWkMqcyqtfVwNeH2RcoFrUozFbCRhS0iC+fMHYg3UD7EJRwvGTtWMsIazQajWbqY8OYEW6P1ECij7WMsB3IDr+ZanR0TIp1qEZHq/p0dHR0VFCPqraQTzgWct9sRWukrKvZGPHrijgzBUaVsIWIsyuv0Wg0moTHqBK2ENYywjaUqjqOKiacWconrFz0Xa0ZYKmNGFR9u6rtBhFK7WJDDHzCZkdWx2LjByuVrYxFRLYGsJoR1mg0Gs3URyvhOEVRCY+oll5RVbOgPoO2paq1U53JxiQ6WhV1Ba2qaM3uz1I+4VjkCeuKWeMzhSpmaZ+wRqPRaDSxQivhOEWgpIS9qu45IypRVUXbFOuF+2KghJUxd1tBI6hGOQ9aqXZ0LFSQVaqtGbnlrJQnHG8VsxIIfeU1Go1GE1/o5eg4RfuEx8ejWGg1QaKjVRWtKqo+YSNq3/R1AivVjo7Jio8iOlJZL0drNBqNRhMztBKOUxR9wiOq+XNuxXZgfp5wAlTMMlvNgtVqR1uoEplVimsbueXMVrRTqWJWAilhC9Vk02g0Go1mahFn059JMNsnbGQWrKwwFbd8slSesJryMpLra37FLLWbx4ja71duqUgs1OyIYjvFUAtDqNad1hWzEkoJx9mV12g0Gk3CI0gY62SttynU9gZWVsKx8AnbFQdrqdrRmuhgIZ+wVe6dhImqjrOKWdoIazQajUYTI7QRjk+kohI+rNphLPKEVZWwqh9J1ccGmK28YuETVvXRpinuomTkPZruE46FmlXtMxa1o5W/A3TFrERCX3mNRqPRxBdaCccpAvyu8Jv1qfYXi9rRLpOjow2pGXOVsJHI4VjkGJuP9glHHCv5hI0Qb9HR2ghrNBqNRhNDEsQ6WeptqvqEG1Q7jIVPWHX3HdWZbEwqF8VZJOYEqFfMMjcvOTYY+TuafL/GQtGq9qmc6WDg69xunc/kVMNSRlij0Wg0CYBejo5PVJVwo2qHRpSwX3FXI9WyNx7FdjHxCasRi+hoa2Ehn7BViEUogWqlralUMUsbYY1Go9FoYoQ2wnGKYnS0shI24kdS3Q0JhTcI6mONgRJW3TM3FmrWSrWjLYWqb1f1fo2FT9jsP6WumGVJTNlFSQiRLYT4gxBirxBijxDidCFErhDieSFEdfAxx4yxaDQajUYTL5g11/g+8KyU8uNCCBeBPU1uAl6QUq4RQqwGVgM3THQSv4BRheDhmERH2wsUG6ruvKLYnaHo6Km/L7CVUF1hUBeJMYiONhsrRVXbpljFrDgcUjSIuhIWQmQCbwd+DiCldEspe4DzgXXBw9YBH472WDQajUZjAY4sR6v+WAgzhjsDaAf+TwixFNgEXAkUSSmbAaSUzUKIwrEaCyEuBS4FKJumtiNSTCpmKfuEFfOELeQTViUWCtrs2tGxwUIVsxKhdrRi0bwphfYJRxQHcBLwiJRyOTBIYOk5JKSUa6WUK6SUK3LzozVEjUaj0cQTfpv6j5UwY65xCDgkpXw9+PwPBIxwqxCiJKiCS4C2yU7kB0YULnCv6ozUiE/YlqrYcOor4Vj4Z60SdWwt33WcRdRGAyOrYaptlXdfMmB9RILIzqMQQmQDjwKLAAl8AdgH/B6oAuqAC6SU3dEcR9TnDFLKFqBRCDE3+NI5wG7gL8AlwdcuAZ6O9lg0Go1GYwEESLv6T4gcCRieBywF9hAQiC9IKWcDLxDGqq0qZk1/rgAeC0ZG1wCfJzABeEII8UUCAcyfmOwkUqj5hJUVrVuxHYBtTBf35PgUK21ZKDpaVe1ZqWKWau1oI6iq/X4r+YRVUd03W7XwndWIs4pZUoRlTMPmqIDhz0EgYBhwCyHOB1YFD1sHrGeSrB2jmHLlpZRbgRVj/Nc5ZvSv0Wg0Gmth0LebL4TYeNTztVLKtUc9NxQwHEnia/ozCX4UlXAsIipVfcKqM9KY+ITVUI84tpK/VA1rvccYKOiY7PqliOr3h2o7I/5riwUzhUCHlHIs4XeEIwHDV0gpXxdCfB8Tlp7HYupdeo1Go9FYm+j7hMcKGD6JYMAwQKgBw0axlhIWatHRyr7dmPiEraSErRMda7bCVK0dbQTV99gf4XGEhNmKVvU+z4joKOKXOMzriaZPWErZIoRoFELMlVLu462A4d0EAoXXYFLAsKWMsEaj0WimPlKYMi+ISMCwUSxlhH1Aj8ofRlWUxCJP2KMYjqk6VgtFR1sl19cIsXmP1lnRiEUMgzJm5wkbMVpRVJ2qRFMJQ/wEDMffGoRGo9FoNAmCpZSwG8W9gVV9u4ZqR2eqtVPdzcRCPuFYRAA7TFeYiTC/tZCCVl3xicUuSrEg3t5nlPOE4wlLGWGNRqPRTH0kcRkrFhUsZYTdKO4NrOovNbSfsOIdpJonrCz0jKgZc5WQlXJoYzFW9T4tVDFL+4THZyopR62ENRqNRqOJHdoIxyE+FPcGVlW0RvwkdpdaO1UlrFobNwYqKBYVs6RiW2HyvsBWUvuW8glbSUHrilkJhaWMsEaj0WimPlKAL0EmBtoIazQajSaukIDX3EWomGEpI6wcmKW6TGMkMMulOI1T3cpQeblNfUnRrnhhrbTkqrqMHYv3qF7oIxbLyop9Dpn8lWUks81KtWXicKyJYoQTRPBrNBqNRhN/WE4JKxXriEWKkmq+gGpgVgy2eDO7xKKRght+HZgVBSwUmGUldNlKpABvgkhESxlhjUaj0Ux9tE84TvFgctlKFP2zADjVmqmWiYmBT1hVtVlJ7an6hM0vk2k1VO87xa8s1c+Hcuofyl8Byu2MEGdlK7UR1mg0Go0mRujl6DjFB/SqzNhUfbv2AsWGoOxkMX0jBvOVsKovWVWVAvgV+7RZKHZR9e+hGuUeE/FkpaIbqsSibGWc+YQTCUsZYY1Go9EkBno5Oh6RqKla1ZmlSFNsCMqOHdWxKkdHWyfC1YgSVsVsBW12xLkxrHPvKH8+jCjvVMV22iesfcIajUaj0cQKifYJxyeqSthKPmHVscbAJxyLjRhUMXsDh1ioditFnSvfd4ngE9YbOAQCsxJECcfZpddoNBqNJnGwlhL2o5bzq6ouDfmEdXR0pNup+mc1E6O6otEf4XFEFdXPh5FKdDkG2ppNnH20tE9Yo9FoNJoYon3C8YiqT1i1YpYRn7DqbkiqM9IYREdbyQcZCx+tCla6psaio1V3UTLQpVXQtaPxAQOxHoRJJMhcQ6PRaDSa+CMxlLCqurQVKjZEfTckC/mEVVGtq+yxkEqMhfK2Vo6xyah+PoxEY6vWndZ5woHqiLEehElYywhrNBqNZsrjA/piPQiTsJYRlqj5d1V9wjbVkjeAW9EnnAB5wolQ+UqVWPiEVfvst1LFLCuha0drJRxphBB1BDIafIBXSrlCCJEL/B6oAuqAC6SU3WaMR6PRaDTxizbC0eGdUsqOo56vBl6QUq4RQqwOPr9h0rOoCBrlilmZig0x3yecANHRVolwBp3TPDkWqZhlpD/VtgYW4JSJM59wIhHL5ejzgVXB39cB6wnFCGs0Go1mSqN9wpFHAs8JISTwUynlWqBIStkMIKVsFkKMGYoshLgUuBSAEsytmGU34A/0WyU6Wh1VJZwIijYW+xCbXcHMUlgpv1jXjtbL0VHgTCllU9DQPi+E2Btqw6DBXgsgFooE+LbQmMVw6gBd+a0kjaSQ3VmIw2etOEWNZqpihhGOl1glU751pJRNwcc2IcSfgFOAViFESVAFlwBtk58Ic/cTxqXaEPwmV8xS9LHZLeQMMuJnlcF/w2kDdBa20FXQwmDmWwtewm8jqyuP/NYSctuLcXgDyZqquyipoppDHRtiER2t2qeeYE1InN12Ji5HRyZWyQBRvzOFEGmATUrZH/z9XcAdwF+AS4A1wcenoz0WTeIhkfRnddNZ0EJXYQsjqYF1yYyebKr2LSCvrZiRlCG6ClrpKmzhwMLtCP8OsjsLyGstIa+9GIcvFtUTNBpNDDA9VsmM6WER8CchxJH+HpdSPiuE2AA8IYT4ItAAfGLSM5m9n7ARR4mqT1i1yo5iQyNVllTbRjty2Otw05PbQXdeO935bXiSRhF+QVZXPqX1M8htL8J1VB530mgKWT15VFXPZyCzl46iJjqLmukuaOOgbwc5nQUUNk8jp6MQmwztnrCS3zs2mKyiY7GLkmrbFMV2cebXNUIElqPzhRAbj3q+NujaPBrlWKVIEnUjLKWsAZaO8XoncE60+9ckBv0ZPXTnt9GT105/VjcIsHscAUXbXkx2R8GkilYgyOjLJqMvO2CQs3roKGqmo6iJrsJWnG4XBc1lFDZNI23QQPqaRqOZkAgY4Q4p5YpJjlGOVYok1nKUqFbMEnmKHRrwCavmCZtc+cpIZKzZ0dFjtZNI6ubsprmiDiSk92VRXjuL7M4CMvqyEUHlesQnHA7pvdlk9OZQVT2P7rx22koP0TytjqbKWtJ6syivm0Vee7HSe4kGqn8PXXNaE2/4ib5POGKxSgaxlhHWaI7Cb/NRvXAbnUXNFDdWMq1mDk6PgYnTOAhpI7ejiNyOIjxON+3Fh2ktr2ffkk3M3rWUgpbyiPep0WiiRzzFKlnLCKv6hJVrQBtwsvgUjUGnaofmK+FY4nG42bd0E305XVTtn09pwwxT+nV6XJQ2TqfocAV7lm2geuE2hN9OfltJRM5vLV9yAkRHx6Jilmo7xYSMeMSEFKXIxSoZxFpGWKMBhlMG2bNsA6Mpw8zesYyC1jLTx2D325m/dQW7l79B9aIt2LfZyemMegyHRpM4RDF7Mp5ilaxlhJWVsOqXo4GlTYvsC2zEH6iqov0GfMK9OZ3sX7wZECzYdAqZvbkhqUfViGz7BNvL2P0O5m9dya6TX2Pfkk0s2Hwqmb25b45VMxEWiY7WxAZJ3OUuR4spFNSumcpIJM3ldexe/gZOj4vFG05/0+DFEofPyfwtp5A0ksKeZRsYSuuP9ZA0mqmBz8CPhbCeElaJjrYXKHZoxCes2M7k3ZCM+IRVqzuF+yf023wcmLuD1rIGctoLmbVzqWkFNEJRtE6Pi3lbVrJz5f/YvfwNFm04nbTRDBNGdyw6yjkKxMInrIpyjQGmlD/ZamglrIlrhlMG2bryFVrLGiirmcXcbSfHZQWr5JFU5m9Zic/hZe+yjXjtyhViNBqNRCvhuEVlsq+shA34hFW/gy3ku1L1e4bqn+0saKZ6wTYEgvlbVpLZmaeU72sWaQOZzNm+nD3LNrJ38SYWbjvlzTzleEZ1NcRI3XGfsk9YsZ2VdlEyu9JWvJIgCzvx/w2hSTh8Nh8H5+5g79JNJA+nsfT1s8jtLIr1sEIiu6uAGXsX0Z3fRs3sXbEejkZjTbQSjlP8KFbMSlPsz8DlsUh0dCzyhCdSsoNp/VQv3sJQej+l9dOpODAPm7QZVr9mqueipmmMpA3QVFlLymA6xYcrQ24brypfo4CqolWNNzSyiqZ9wjHDWkZYM6XpKGqiesE2HF4H8zefQk6Xqhsh9lQdmM9Q2gA1c3eSOpARF5HcGo1lOKKEEwBrGWHVPGFVn7DPwPRQ1Z9hcnS0kYjaSKq2wfQ+DizYRnpfFnO3n4zLkxTV/kJFNb/YIR3M3bmcbaf8h71LNrHs9bOO2blJYxGM+JLNVrSqhQHBWGR1tNA+YY3GHLx2D/uWbMLudTJ3x9gG2Io4vC7mbVuBz+Fl3+LNSJEg3yoajVG0TzhOUVXCyYrbzhnxCav4rgG6VDs03yesqhKPVrQSSfXCbYwmD7Ng06k43a6oKF7VcwqEoX7TBjOZuWcx1Yu2Uj9zH1UH5hs634R9Kb7HWOy+pF7SRDE62kJZB8oYUe3xGFltMWOqilbCmphyuOog3YWtVByYN2X9poUt5RQdqgi817yo74ym0WgshPWUsIrCTFaca/invk84lrsodee10ThzP/nNpZQ0VE16fCx8wpHa+7hq/wL6s7rZv3Ary14b3z+surpgPWKxA5MCVtpFyYiajbeVggSqHW0tI6yZMgyl9VO9eCtp/ZnM2LPY8LLvmNjBmWXHmWvHmeMIPtpx5joCjzl2HFl2+rYM0/R4N96e6K1/2f125uxYzvZT/0P1wm0s2HJKdN6zRjNVSJDlaIsZYQF+laAdxcpXCZAnbARVlTjqHGHv0o3YfDZmb1uO8AtDCtCebiP/3AxcBccaWkemDWE71tBJr8TT48PT5WO02cNwvZucM9PIPjWV1j/20vb3XmSULmXqUAaV1fOpnbeL7oJWctuLo9PRlMYiCjoWGFGzRiKro4FOUdJoooPf5qN6yRbcSaPM33QKSaPGIkLS5iVR9Y0CXPmOgHHt9uLp9jFU48bT7cPT5cXT5Qv83u3F2+fn+LlD8xM9lH0mh9KLcshdlU7j2g4G9owaGtd4FB+uoHlaHQ0z9pPTXqTVsEYzHno5Oh6xKVa/Uiz4b8QnrFy/X22mr1rH18w8YYlk34It9Od0M2vHUjL6ssNqf4xaFlD8kWxKL8hhtNXLnhsOMVSrFpI+0urm4HdbyVyaQsWX8pl9ewkdL/Vz6NedMKB0ynER0sa02llUL9pGV2ELeW0lETt3LKKcrYOqgjbwFamqTGNRO9pKtbWnGBYzwhorUzdrD+3Fhyk/MJu8VnXj4yp0MP3yQtLnJdP5Sj8Nj3bgHzYetNW3bZhd1xyi5GPZFH8wm6zlKTT+tIu+TZGNWslvKePQ9AM0Tq8mt61Yq2GN5nj0cnS84lSsfmVX686IT1i5Ao25la/Mio4+PK2GQ1UHKD5USUnddGV/cv45GZRfkgd+qP1RG13/jqxUlW5J02+76X51kKqvFzDjhkI6Xxrg8C+7ImLoIZB7XH5EDRe0kmdR37CRe6ffKr7deIsaThS0EdZoIkdb0SFq5uwkr62YWXuX4FbIM3MVOii/tIDMJan07Rim7sdteDqj9ykdrnez98bDlH48l6IPZ5G+MJn677UzdEC1Csux5LeW0jijmkPTD5CrfcMazYkkgpcEXaxDE2W68lrZv3ALmT15zN15cvjGxgYF789g7ndLSJuVTMOjHVTf2RxVA3wE6YOW3/dSfWsLAph9RzEFH8wkEvZSSBtldTMZzOylN7fD+Ak1Go0lsZYSFg7F5egYBGZZZEtCI0uK3kmmqj05HexZsoHUgUzmb10RdipS2twkyr+YS0qVi95NQ9Q/2m6K8T0aiWSwepS91zcx7at5lH0mh7R5SdT/qH3C5elQltvzm99Sw1ld+SG3Gw+tpROcqbR0nkDL0VoJa6JCX1YXe5ZtIHk4jYVbTsXhC30i5MiyUfH1PGZ/pxh7uo3aB9qpvdd8A3w0viE/dQ+2c+j/Osk6KYW5d5eSVKo4uQtik3ZKGqroy+miP7MnMgPVaKYKegOHOEQ4wFao0FBVCRuYo3SqNjRXCTsMOF7GG2lfVje7l79B0kgyCzefhtNzbLGU8dSecAoK35dJ0UeyEC5By596aP1jL/7RwPGjMdhvLeW4Kgbt/+hnuM5N1dWFzPpWEXuuOWwoYKvocAWHph+gqbKGuTtOMjrcBEA1oCsGKUpmpxrFW8ENIyRQ2UqthDURxeMcZc/SN3C6XSzcfBoud4gVzmww47pCSi/KoX/XCHuvOUzzb3veNMDxxMCeUWrubcWZa6fsM8Y2nbD7HBQdrqCrsIWRZJ2sqdG8SZwrYSHESiFE8VHPPyuEeFoI8QMhRMhfDNZSwmanKBkJhDXZJxyLggtjKdqaubvwObws2nQ6rtHkkH2c5Z/PJXNZCg0/7aDzhQhXyIgCQwfctP21j6Lzs+h+dZCBneoqvbixiuaKWloq6sjYv1T5PKqrGrGIJ7BM+Ukj8yKzNwXTcziz+SlwLoAQ4u3AGuAKYBmwFvh4KCfRSlgTMbryW+kobqK8dhapgxkhtyt4XyYF786k9eleSxjgIzQ/0cNIs4eKS/MQTvWwqKTRZPJaS2grPYTXrlxqTaOZOhwJzIpjJQzYpZRHdoC/EFgrpXxKSnkLMCvUk4SlhIUQTimlRwgxE+iUUvaE2M4ObAQOSyk/EJTqvweqgDrgAill9+RnUixbqbTpA8Z8EhbZktBINO7RbT0ONwfn7yC1P4PSupkTnvfoCOmslamUfTaH7tcHOfT4xI70YQaVx5qCSrnTSbYW9EDD2nbm3FZK0cczafptd2jtxqCosYKOkiZaSxopOzRDaayaiYiB8rbSVobxSPz7hO1CCIeU0gucA1x61P+FbFvDVcJ3CyHKgDuBH4XR7kpgz1HPVwMvSClnAy8En2ssTN3c3XidbmbtXoJNhnZbpc1NYsaVhQweGKX2h20nbKxgBfp3jdD5cj9FH8w2FC2d3pdNWm8WTdNqY7JvskYTV1hDCf8WeFkI8TSBqdMrAEKIWUBvqCcJ1yecCZwP3ANcHEoDIUQ58H7gLuDq4MvnA6uCv68D1gM3TH4yxTxhn6LrOyY+YTXfYiyU8BE6C5vpKGmi/OBs0vqzQmqTPM3JrBuKcXd4OXBvC9IdXcOjqqKTmDxX/NBvusg6OZWKL+RRfWeLUj8ARYcqqFm4g+6cdnK6VWIfEgGTFa2Vcm+tNNZQiP9Uo3sJiMgS4Dkp5ZEvMRsB33BIhKuE1wOpUsrtQHWIbb4HXM+xiwtFUspmgOCjSt6RJg5wu0aombeTtN4syupmhtTGVeBg9s0l+N2S6rta8PXH/7rTRHh7fTT9vpvMJalkr1TPE8lrLcbhdtJSXhe5wWk0mmjxhpTyNSnln6SUb87ypZT7pZSbQz3JpBJRCHE2cBHQA+wEXhJCJEkpfxZC2w8AbVLKTUKIVaEO6rhzXMqRtfb8HJRyfn2Kla+M2AaT84RVo6PD9V0ejQ8fBxZsx2/3MXPXYpCTn8+RY2fmLcUIJ+y8rZ7h9ujs22s27c/3UXBeJuWfzaN3y5DSn9Hmt1PUXEHTtBrcrhFcbgMV2zQaK2NSnrDBeKWIFKkLRQn/Bvgb8BowA7gV2BXi+c8EPiSEqAN+B5wthPgN0CqEKAEIPraNdwIp5Vop5Qop5Qoy1YJrNNGhZVo9vfkdVFTPI2UofdLjHdk2Zt5WiCPTzp67Gxg+NDUMMAB+aFzXSVKRk8L3hrYkPxbFhyuRNklrSWMEB6fRWBBzfMJG4pUKhBBXj/cT6gBCcZYekFL+Kfj7k6GeGEBKeSNwI0BQCV8rpbxYCHE/cAmBvKpLgKdDO6MAXJMedQKqNaAtlCds1paER+jP6KFx9j6y2wspPDRt0uMd2TZm3lqEM9fO7rvqGTgYvu97yEAiZKoJ5YT6dwzTu3mI4o9m0/PSML6B8KfyqUPpZHbn0lLWQHn9LFN2V1K9d4zlplskT9gIqt8BZlfaikdMqB0dgXglO5COQUUcihF+WQjxTeB7RzmejbIGeEII8UWgAfhEhM6rMQGv3cOexRtwuF3M2L1oUkPhyLEz87ZCnLl2au9pp3/fVIsgeYtDj3Wy4P5yCs/PpPmxHqVzFDdVsH/hVvqyusjqzYvsADWaxCBfCLHxqOdrpZRrjzvmewTilY4uanBMvJIQYqJ4pWYp5R1GBxqKEV4ILAJuEEJsArYCW6WU4ari9QRmFUgpOwnkVYWJQMkn7FeMjjYy0TdZCasSbnS0RLJ//laGUwaZv+kUHB7nhOdwFdiZeUtgCbrmrjYG940ayvdVRVVFh5tfPNo4QvsrveS/N5P2Z/rxdoc/nc9vLeXg3B20ljWEZYQTI7XJ7NrRRj6PJhcknGpzW2M+4Q4p5Yrx/jMS8UqMo4CFEGcCn5ZSfj2Uk0zqE5ZSflRKOQeYDtxGICr61DAGqplCNJfX0l58mKqD88noyZnw2KQyJ7OCOyEdvKuVwX1TyAc8AY1PtCNsgqKPZCq1t/sd5LeV0lHYjM+WAMu2Gs3xRD9P2HC8EkcJSSHEMiHEfcHzPUiIKbwQxlRNSjlMIIps42THRg+BUh1oVZ+wEZ/EFKyY1ZfZxcE5O8ltL6K8bhajjG9UU2e5mLm6COmTVN/Wwkij8XKMsaiYpcJou4fOlwbIPSed1qd7w9qC8Uh0eX5zGa2ljXQUNFPQWhatoWriEVVFO9V2UYqiTzhC8Ur5QojLgU8DAwRiplZJKeuEELWhjkXXjtaEhNs1wp4lG0kaSWHOruUT+oEzl6cw+7ZifEN+9t8aGQNsNVr/1IsQUPghNTWc1Z2HaySZ9uLDER6ZRmMR/AZ+1FkDnCeEqAbOCz4fj70EArs+HszguVdKWRf8v5DVjcV2UTLZJ2zEdijOZu2K079o7qLkF372LN6I1+lh2Ya34fSOH6GeVOJg+tUFjBz2cPCeVry9J47LSJSzVfB0+Oj69yB556TT+sfeMa/DRAgEBS2lNFXU4nG6T9iT2eqo3ufxX0TpKMyOcp76H6uoYCBe6WPAJ4HnhRD/Ap4AnpVShmU5tBLWTIhEcmDedvpyupizexlpAxPkwNqg8vIC/B5Jzb1tYRueqUbrn3sRDkH+e9XUcH5rGdIm6SxQL4Wp0VgSC9SODlbKupDAjknPAl8BDgkh/o9AieeQSBAlrOgTVt8iVrmxqqJV9QlPVuGqaVoNrWUNlNfMJq+15Jjjj29b8tFs0mYncfDBVka742sJWt2frF7D2d3ipff1IfLflU7bn3vxj4T3N0rrzyR5KJWOoiaKmyqUx6GBhIiOnmpYZNkjWLLyMeCxYMWtTxCouBUSWglrxqUrv5XaObvJbSumombOhMemzU2i5OM5dL7cT89r5qcgxSttf+3DkW4nd9XkFcWORyDIby2lN6cTj9NI5RiNxmIcKVtpvk/YEFLKLinlT6WU7wy1jcWmaooVszyKcw3l+s9gdpSzajvPOO2GUwbZt2gzaf1ZzN65FBg/ktqeZmP6FYW427w0/KJj0j5VVWksoqNV+0wN9jd0wM3gvhEK3pdBxz/7w96uMa+thEPTD9BV0ELRJGpYNU/Y7NUXY1goZUvvJ2wMiyhho2glrDkBiaRm3k4EgnnbTsY+SWBb5WUFuHId1Hy/Df9wIhSMCI/2Z/pJKnaSsTR8t0hafyZJwyl0FrRGYWQajSbWWE8Jq0Q6q86oDFWgMVcJOxTVzFie246iJnrzOpi+dyFJo+NPr0cZoeS9ueSckkbtuha6D4a8j/WU52hV2vPGIJ7uHPLfk0Hf1vBuKoEgt72YlrJ6fHYvdtW9sTVTn6lUMcuE2tHxglbCmmPwON3Uzt1NWm8WxYcqJzw2Y04KlZ8pomtDP81/7zJphNZD+qDzxQEyl6XgzAu/2ExuexHS7qcnd/Klfo1mymBBn7AK1ppWS5tapLNqoG4MlLCqf07VH3h8u7rZu/E5PMzcc8qE53Vk2ZhxdTmjHR6qH9YFJSaj88UBij6SRd7ZGbQ82TPucWNd74yeHOweB135LeS2F0VxlFOZWERHa9QRQJKB9oZSW0xFK2HNm3TntdNeepjS+hmkDYyf5ibsUHF1Ho40O/seaMQ3ZLGpZwxwt3vp3z5C3qr0sDc+s0kb2Z0FdOe3J8gmDRpN4mAtJYwAn4ISVvUtWCg6WpUjub5eu4ea+TtIHkijtHbGhPnDZZ/PIW1+Eju/d5DO+vCXoVUjjmNTaUs9T/h4OtcPMP2qAtIXJjOwM7yZek5nAZ3FzQyl901cMMUiqK749FtKmSqOdchiX8tRwQbCSL136yhh/dfWANAwZy/upBEWbj8Nm398v2Xeu9PJf1cGbU/30fqq9gOHQ++GIXxDfnLPSgvbCGd3BiYDPbkdU8IIazQTYwObkR0pDCkoU7GWEZY2NSUcE5+w2kzMjN2Qjqc7r432ssOU1s4gvS973OPSlyRT9rkcejcO0fx4j3J/VkJVfeeO8Zr0SHpeHyT71DQaH+1CekL/m7ncySQPptGX00VZw0ylMUWa2OQJaxIDu0ElbB20TzjB8Tjd1C7YSUp/OmU1s8Y9LrnCSdXV+Yw0emj4QWfYRSc0AXpeG8KeaiN9XvhBJ5k9ufRld2m/sGbqI4LL0ao/FsJaShgB0sQ8YUMuSHOjoyerAT0WEknN/J14nR7mbl6BkGLML3hnrp3pNxbgG/ZTs6YN30igL1WVmAi7KI3HwO4R/B5JxtIU+neEt1qS2ZNLW1kjQ2n9pA2qbQoxFlrRTkQMfNB6P+GEwmJGWBNJmsvq6ClsY9r+uaQOZox5jD3NxoybC7Gn2Ki+rQVPV/xk0AshmLGkiuVnL6V8ThnP/uJ59m7YH+thTYh/VDK4b4SMJSlAd1htM3tyAOjP7o6oEdZo4g+jPmHrYC0jLBWjo1V9wobijsyNjg53iXIwrY+Dc3aQ2ZFHUUPFmO2FSzD9hgKSip0cvLuV4frYbiIgU33MWTaLeSvmsvRtizjpncvJzg8EKY0Oj/Khr76XF55Yz4+u/wkt9ceWeUyNI5nQv2OE0k/l4Miw4e0PfQUjaTgV56iL/qxuig+fWEhFZTUkcYhFVLVqn4pfy1NqgclodLR1sJYR1kQEn83LnsUbsXudzNi9CDFW4qoNZlxTQNqcJOoeamdgl7kh/06Xg5lLZzBnxSzmrpjNnBWzqFpYgd0eiNxubWzjv3/7H5te3MLGFzcz2DvIp665gIuuu5AzP3Aav7jjV/z2wSfxeuIvpaU/GBmdvjCZntdC/+YUCDJ6c+jP6onSyDSaOEFoIxynKFbMspBP2AwlfGDedobS+lm05XSc7rEDhAo/mEnm8lQa1nbQ8/rYF8LIrkZjMWvZDD741fcyd8Vspi+uxOkK7B3d097Lvg3VvPznV9izcR97N+6jo/nEFISf376Ov/78Gb75/cv52ppLed8l7+aeL3+X7f/dGdEo51CYSJUO1AzjG/aTtjCJrtcGwjpvel82XYWteBxunF6FHcXiBNX7vD/C49DEK3o5WjNFaSmtp7W0kWk1c8jpKmBkjFSqpDInJZ/Ioef1QTr/FZ6RUOXU963gtidXI/2SXf/by5MP/Jl9G6vZt6GatsZ2IDSj33aonRs/dhtnvO80rvnRN/jxy9/j12t+y2PffhKfN0782T4Y2DdCxvzw957L6M0GYCCzl5yuyBUS0Wg0scFaRlgCfoWsKlVFa6Ha0aEwkNHDgbk7yO7Mp7JmLjCGYhNQcVke/lE/9Y92RMXPeLwx/eAX3sf1P72a6q0HuOb9N9LdFl7A0li8+sxrbF2yjau+dzmfu/liTnvXKdz5qftpOths+NyRYGDPCGWfysWeZsM3GPo1TusP+MAHM3sS1AjHn3sh8pjsS4Y49CfrPGHNFMPjHGX3kg04PS7m7jx5bD8wUPShLNLnJNPwiw68vdFVjjabjSvu/yo3/fw6Nr6wmcvf+c2IGOAjDA0Mc/eX7mf1R2+ldGYJj2x4iJPPXRax8xthYF9gBSJtTnj5wg6vk+TBNAYy9baRmqmMzhOOTySgEqCrrGjVg5Hsio5oVV/ZRIpVCj97Fm/E7Rpl8cYzcHicbx4/etR7TK1MovTCXDr+10vLfycv+2bEJ5yWmcbtj93MmR84nT/86E9875sPR225+OU/vULDlsPc+ZdbWPPs7fzkmp/z1Pf/EpW+QmXo4CjSJ0mfk0zflrdu0FB8+2n9mQxk9eiiHaZgRHmbrGiNrNyF7xmJLiJxfMJaCScAtXN205vbyay9i8nozx7zGJtTMOeKMrz9XmoebYnqeGYtmcn/bfwJp77nFO677CEeuOIHUffXttS1cvnp1/Lq06/z9e9dypfvuSSq/U2Gf1QyXO8mbXb4lbPS+jMZTRnG61DNvdNoNPGC9ZSwyveOci1v9VlwgclK2DNOu/biQzRPq6OkfjoFzeXjqqfKiwtJrUhm9131ePujZxDP/dwqLvvxF+nr6uPrq65i+393htzWiPJOIY2RwRFu/8QavvGjr/Kp1Z8gLSuN73/9EaQc/5pHOgL8aAYPjpJzevhLZ6nBbSaH0vvJ7FGN3zZONOMXIk8i+JKnEjpFSTMFcLtGqZ27m4yeHKqq5497XM7J6ZS8N4+mv3fSsy06Ric9J42vPfJl3n7hGWx8cTO3furOiPp/Q8Xv9/O9r/04kFe8+hP4vD5++I2fmj4OCBjhgvMySSp2MNoSupFI6w9UNxtM74upEdZooocAxt/NbSphLSMsUcv5VfaVqM+eVRWtQ1FduMfor2buTnw2HzN3LwbG9jf6c93M+lopfTWD7H2sNiw/Y6i5t8vPXsrqdd8kpyibn924jkfuewS/P7Yq6mc3rsPusHPBtR+lfk8jf3nkGdPHMFQ7CkDq9KSwjLDTnYTd42Q4XWfNTk1iEB1tKBMkWiSGtzQx3mUC0lnYTFdRC+W1s0gZSh/7IDssumomwinY+f2DSG9kA33SstK4Zu0VPPDCXQwPDPP1067ht2uejLkBPsJPrvsF217ewUU3XYDdYf6se6TRjfRKUqrCK7ohEKQOpDOUZk4Ot0ZjPkeUsOqPdbCeEjbVJ6weHR1Lf5nHOUrNvF2k9WVSVj9j3OOKP51F9rwMdn7vIMPNoxEdw1kfOZ0rfvRVcoqy+d19T7Hu248zOhzZPsJhPN/uunsf48Fn1nDGBafw3OP/GuOI6OXiSi8MH3aTUhF+5auUwXS6CluQyHHTzULFWrsome3btZIv2chYrWUKphJRV8JCiGQhxBtCiG1CiF1CiNuDr+cKIZ4XQlQHH3OiPZZEQCKpmbcLn8PDrF1LEXLsP3HmaSkUfDCDxmdbaX3V0E4Vx2Cz2fjW49dx+x9vpru1h6+tvJq1N/xfTA3wRLz27BvU7annwqs+FpP+hxvcpFaGHyGdMpSO1+XB64ztphoaTfTQSjhSjAJnSykHhBBO4D9CiH8AHwVekFKuEUKsBlYDN0x4JlUlbCGfsGru55G83/biw3QVtTCteg7Jg2lj5g8nT3NS/rUcBveNsm3dPuU+x1KXl935Jc7+1Dv42W3/x7q7H4ufUpHjIKXkqYf/zDU/upJ5J89h7yZzt0IcOeQh76wMbCkC/3Dof4eUwUDk6HDqIM7e8I14rLGW+lYlBr7dKYMgUbylUX+XMsAR55Uz+COB84F1wdfXAR+O9limOqPJQ9TP20N6TzYl9dPHPMaebqPqunz8w5K6BzuQvsh9Gb7r0+fy2Rsv4s8//Su/uONXcW+Aj/CPXz/P8OAw51/6QdP7Hj4UULLJZeEtSScPBYzwSGr0Uqg0mtiilXDEEELYgU3ALOBhKeXrQogiKWUzgJSyWQhROE7bS4FLAcgqVYuOjkGecCyU8IGF25FIZu4cJxraDpVX5+PMc3Dg2y14uiPn71p02gJu+vl1bF6/lQeu+MGEx6ruaGSEFMbPORzsG+TFJ1/m3AtX8b2rfmTq0vloU2BpJ7nEydCB0PtNGkkBv2AkNe6K/sYpieDbNfJ1Hm/XJ7opSkKIZODfQBKBC/cHKeVtQohc4PdAFVAHXCCljGoupSl6X0rpk1IuA8qBU4QQi8Jou1ZKuUJKuYJUnRM5Hk3TDzKQ00PV3gUkjYxd7q38C7lkLEqm8aedDFVHzpdYPrOUe5++k/bgDkbxuIfvZPz9l8+SnpXOWR86w9R+R1s8SJ8kudQZVjubtJE0ksxIijbCGo0CR9ykS4FlwHuEEKcRcIu+IKWcDbwQfB5VTF10l1L2AOuB9wCtQogSgOBjm5ljmUr0ZHdweMZB8ppLyGspGfOYgg9mkn9eBq1/6qX735FbwszOz+LBf9yLzWbj6vetpq+rL2LnNpOt/95O++F2zr3wnab2K33gbveSVByeEQZIGkllNDkuEzw1mghgM/AzMfHkJo36crQQogDwSCl7hBApwLnAvcBfgEuANcHHpyc9memBWeanKIW7daDbOcquRW+QNJxK5d5AVazjl6GzT0+l9KJsul8dpOl3x66sGCnLmJqRygPPrKGwvIDLz76axupDyueKNVJKXnzyZT5y2YdITU9haMA84zba6iGpKPyPYtJwCj357W8+V3VlGEtw0sQPUylFyfBydL4QYuNRz9dKKdce04MBN2kkMUMJlwAvCSG2AxuA56WUfyNgfM8TQlQD5wWfa8JAItm7aBMep5tZ25dg9534QUqbk0Tl5QUM7hul/uEOIhWU6kp2ct/TdzJn+Wy+dcHt7Hxtd2ROHEPW//EVXEkuTnvPKab2O9ruxZWvYIRHk/G4RvGL+Ch+otFEDsPFOjqOuDGDP2uP78GImzSSRH36I6XcDiwf4/VO4JywTuYFVFzkUzQwq27GXrrz2pi9e+mbRf2PxpYkqLw8H0+Xl5r72pCeE8+tooQdLgfX/+EKlr9jKd/6zG08/7exilxYjx2v7qSno5e3fegMXvzDy4bONRrGKspQ+wgF2Zl4XKO43KFHSbtGUkCAO2mE5HHiADRGsdBWhholgqu06znKTRpUwaa4SRMjEWsK0pnfQsOMfRQ1TaO4qXLMY0ovzsFV6KD+xx34BiOjlhwuBzf94RpWvv8k7vzKPTzz+D8jct54wOfz8/o/N3Dqu1YghHmLtO7OgI/FlRPel69rNBkAT1J8FkLRaIwRPZ+wEKJACJEd/P2Im3Qvb7lJIVQ3qUGsNeXyg1p2i+qMND6V8FDqAHsXbiKtL4uZe5cAJ/qSM5emUPDuTFr/2kP/nsj4N13JTlY/eTWnfOBkHv7qz/jjo2r3p6ofeqI0o0jxxnMbefdF5zJryQyqtx2Men8Ao0eMcJ4T2Rp6O9dooEiHO0k9dsFq2BW3CLVGxnosibeMhqjvolQCrAv6hW3AE1LKvwkh/gc8IYT4ItAAfCKagwCrGWENXruH3UvfQEgbC7avxO4/8UZ1ZNio+loBww1uDv8uMiluyWlJ3PKXG1i8agE/vPSn/PNnL0TkvPHGppe2ALD8HctMM8KensAXoDPbQTiJY0530Ai7tBLWTEWiZ4Qj6iY1iLWMsA9QKnOsphRUZ90QnehoKfzsXbyRkZRBFm4+DddI8pvHH+2DnP7VadjTbey6u54Rz8Rf0KEUzsjKz+Q7f/sWc06exT2feZAXHl8fclur0drYRlNtM8vevoQnfvCUKX26g0bYlePAHYYicbpdIMHjik39aGuVnow3pTcRsRhrvJkCAf54G1N0SIx3OUWonb2b7vx2Zu5eTFZP3pjHFL87h9yVGdSua2Go3rhCKplexJpn76BwWj63fexuXv3L64bPGe/seHUXJ61aalp/vkE/0idxZtgJ5wtYIHB4XHoTB43GwljLCJvsEzayHWGkfcItZfU0V9RRUj+doqaKMY9LrUyi6jNFdG/up/nvxndGmnfKHO78yy3YHXauPedmdv1vr+Fzxopw/NBbXtvKuy86l4yyNIYPm1Ob2dPvw5HpIFDIJ3QcHifeGClhI8Ryq0+NBZACpLXMkyo6OtoCDKb3UTt3N9kdBVRVzx/zGHuKjblXl+Pp91H9cJPhPmcvn8lD6+9hZGCEb5x5vaUNcLjs3Rx4r3OXzTGtT++AD0da+D4wh8eJ12GlpVaNJhSCRlj1x0JYa7ReFH3Cal9ShTFQwscXBPPZfOxftAWHx8nsXcvG3cB9+mWFJBW62PztvfT3D4x5zFiMpRBT01O4+XfX0tvZx5fO+DrdbVGtXx4yhqp7EXoe7f7tB4CAEa79uzlVwLyDPhxp4c+JHV4nHr0cHadYKU84Didy2iesiQfqZu9mOH2ABZtPwekZu5BD3gfTKTwti+pfN9K7L3QDPB43/PQaymaWcsXZ18SNATaToYEhmuqamDG/yrQ+fcN+JSVs9zr0Jg4ajYWxlhFW9gmrRUcb8Vs5FNu6j1LQnQUttE5roKRuOlld+WP6gdMWJlF8USZN/2tn719rlMd7hA9f+kHe9elz+MnNj7Ll39smPFZVmcZznvARavfWM91kI5yUH/4mDnafHZ9ejo4isajHrKOj31yOTgAS411akJHkIQ4u2E5abxYVB+eOeYyz0E7F1bmMHvay9ZF9hvtcsHIe3/zB5bz27Bv86p7HDZ/PyjRUN7LkNPNKyfrdfmyu8Kt02XwO/DZdikIzxUigwCxrvUuTfcJG8iBVd7Tx48cv/OxfvBmAWTuWgjwxf9iWIqi8Pg8E1N7Xjm9E7Yv4iCrNLczlzj/cRkdzBzdcdDND0viytpU5XNtERnYGadlpDPZEP0La75bYXLawd9ESPoHf7gu7nUYT3+g8YU0MqZ+7h8GsPmZvWz52YX4bVF6ZT3KZk5q72nC3Glu+cjgdPPDUGrLzs/n8275Mb1evofNNBVobA/Uj88tzTTHC0isRDgUlLG1Im1Se9Gk08YlWwvGJD8UdkdR8wrFQwm0lh2grb6Skroqc9sIxz1N2SQ6ZJ6XQ+NNO+ncG3puRyOGbfnwDy9+2jOsvvIm9W4wva08F2psC+/TmleVSv7Mx6v1Jn8SmUKVP+AMR1VLIiG1TqUlEdFxBrLCWEZ7i9Gd0UzdvNxlduZQfnD3mMQUfzKTgvZm0/bWXzheMLxl/7vrP8tEvnc/a7/yc556YGlsSRoKuYFR4dmGWKf1JCSjs3CRkoI02wpophUQr4fhEojZjM79iVrg+OrdrhB1LX8PpdjFzxxKQ4gQVnHNmGmWfyaH71UEO/8Z46tAZF5zCVfdexj8e/yc/vvWnYbc3or7jne6OHgAy8jJiO5BJOGKEEYlhgVU/k/0xUXpWyhOON7RPWGMifuFj15I38DjdzNu4csx84IwlyVR+PZ/+XSPUP9xhWPUsPmcBV/z6y2z69xZu+8J3jJ1sCjLYF5hgpGWFXuTDCEIQlMMajUb7hOMWiZp/1/yKWaH6hCWSffO30pfdxfztK0nrzzzhGGeenaqrChg55KHmvlak58Rzh6NKZ508g2v/dDmH9jVx6Ye+Rv9of8htEwW/38/w4DApGcnmdGgDqXC7yfBXsDVhofOEY0PiGGFdOzrGNEzfT1tpI5UH51HQVnriAXaouqoAYRPUPNCGf9iYWqpcOI07/nkzfR393Paeu+jv1QZ4PIaHRnAmj12lLNLYHAL/GJOryQm0EdoaazSWxGJTDTWfcMYJFZlDw4hPOBQl3FbcSP3MvRQ2TWNa7Wwk8gRfctknc0mfm0zNQ62MtI5fIziUvX2nzS3nO//6Fu4RN9eccyPNTa2Tv5EExuv24Ewy5yNic9qQXokt3HmxLXCf2WX0NkDXaExHAv7E0IgWM8JTh+7cNvYv2EpWVz6z9ywdc2OG7JWpFJ+fTftzfXT/z1gQVMW8ch548W4QcO25N9Ncqw3wZEgpEQoRyyrYkgS+0fAnfX6bH/wCIRPjC0uTQCRIITiLGWE/Kj5h1XxfI3nCE0VHD2T0sHvJBlIG05m37eRjKmKNBt9fcrGLyq8X0H9gmAO/PGyoGEPVggq++8JdAFy96kYa9r61M1AsIpxDUe1jYWbt6CNIk4Kl7Mk2/MNqRthmUDGo1jnXTCXiLE9YkjBGWE+fTWYodYBdy9/A6XGxYMupOHwnFu23JduYd/00pFey/8FDSK+6IZi+qJKHXr4Hv9/PN48zwJqJsTsc+H3mGCh7qh3fiIIRtnux+/RStEZjVSymhEFlxhYLJewZo63bNcqu5a8jJCzYfCoud9KJClfA7MtLSSlxsevOekY7QvNnj7cv8K1/uAG328NlZ13J4ZompfcyHqqK1iq4kpx4RtTiCcLFkWZjuCn8fYF9dh82nwU/xhrNZCTIAo3+9JqEX/jZt3gTXtcoizaeQcrw2Eur0y4oIO+UTGp/2ULfLmNG7rpHvkn5rDKuOOeaiBvgRCA5NZnR4fANowqOdDvegfDX33wOLw6v/hhrphgJtBxtsU+vmk9YNcrZiK/s6H2BJZKauTvpz+lm9o5lpPVnjunjzTozhWkfy+PwC+3UPGOsXvEHv/A+3nPxeay95RdseXnifYE1J+J0OXEluRjqNUftOzIcePrC/9bxOj3YveHvQ6wxgzjzs05InJkCbYQ1kaRlWh1t5Y2U1s4gv3WMXGAgdZ6L8q/l0r2rj32P1hvqb87y2Vzz8JW88fxG1t39mKFzJSqZOYGiKYMmGGFHuh2bQ+DpCf9L2+N0kzoQ36U1NRol9HJ0PKKWJ6xa+cpIRPKRaOeevHbq5uwhp62Q8oOzx4yadpU4qLwuD3ebl9e/uwOPL/z3eMQnnJWbxd1PfZuuti5u+PRNDPonLsZhpfrPxsZaENbR2fmBjRt62/sM9BkarpzAx9Dd7YExUtUmwut043QbKyjijUF8ppF4CzWspEo1iYTFjLC1GEzvo3rxVlIHMpi5a8mYucCObBszbi5E+qF2TTueQfUvC4fDzn1P3E1BaT5fOOvSNzch0IRPfnEeAL1t0d9b2ZUfWE52d3mB0JeW/cKP1+XB6U6K0sg0mhihl6PjFbXa0ao+YSNKeDRpmH3LNuHwOpm75WTsY0Sw2lIEM24qxJFp4+C323C3eg2pvasfuIpTz1nJLZfczs4Nu5XPo4HCskIAOg51Rb2vpKARHm33kBSGEfa4RgFwjWojHD1iUTs6FsThSoE2whpVPA43+5Zvwu/wsWDDqbjcJ24CIFyCGTcWkFTupHZNO8M1xqJwL/zax/n0Ny7kVw88xl9/9XdD59JAcUURAJ2Ho2+Ekwud+D1+3D1ewjGn7uTAhNQ1atImExqNWUi0Tzg+kdgVpkeq/qdw9wQG8Nm8bF/2KiOpg8zdsoKUwfQTFLVwQNW1BaTOSaLuex30bx9WGt8Rlr93Cdf/4Bv86+kXufv6e5XGrYKqao9F5atwKZ9RRtvhNhyjLhyE73NNInTDmFqUjLvNS5IMz5iOJgfum6SRlLDaaTSWIEGUcNQjMoQQ04QQLwkh9gghdgkhrgy+niuEeF4IUR18zIn2WKKNX/jZvWQDfVldzNy5hMzu3BOOEXaourqAzGUpNP60i97XjEXfzlwxnauf/Bp7t+3j2otuwO9PkOljlJk2axqNBw+b0ldSiZORlvCLgmgjrNFYHzOUsBe4Rkq5WQiRAWwSQjwPfA54QUq5RgixGlgN3DDRiWyo+XdVlXA4PmGJZO/CjXTltzJ79zKy2wrHrIZVeWUBWStSaXy0k86XToxcDkddls4qZvXfr6SntYcvv+8yhgbDV9RTveqVKlXzKnnpT+uj35GA5BInfdvC/9uNpA7hcLvGjDfQJCJTyH8d5cAsIcQ04FdAMYGF77VSyu8LIXKB3wNVQB1wgZSyO3ojMUEJSymbpZSbg7/3A3uAMuB8YF3wsHXAh6M9lmghkexbsIX24iam719ISVPlmMeVXZxDzmlpHFrXRcdzxvbxzS/P485/3YIQglvfczcdrZ2Gzqd5i5yCHHILcqjdUxf1vpKKHNhcNkYOh6+ER1IGSR5KjcKoNJo4wG/gZ3KOiMP5wGnA14UQCwiIwReklLOBF4LPo4qp0x8hRBWwHHgdKJJSNkPAUAshCsdpcylwKYCDfKWc3wLFKVUoSlgiOTB/G22ljVQenEt5w8wx9wXOOzuDwg9m0faPXlr/3qM0niNkF2Zx579uIS07jZveeTtN1c2GzmcVzPIlz106G4B926qj3ldKZSAUa7h+NOy2w2mDZHeGl/9sdVRXtYxNeVUxO+I4ztSsEdSKI4ZM0PYcsT/9QoijxeGq4GHrgPVMskJrFNOy9IUQ6cBTwFVSypArIEgp10opV0gpV9iJr8pARwxwS1kD02pnM612zpjHZSxJofLL+fRuHaJxnTHFmlWQyV0v3kp+eS63v/8eDm6pNXQ+zYnMWz4XgP0mGOHUKhfSJxk+FJ4S9jo8eJJGSRmM/yA3jSaemUgcAmOKw0hiytRJCOEkYIAfk1L+MfhyqxCiJKiCS4C2yc5jtk/YO0FfAQO8nbayRsprZzHt4Bxk8B+8tS9wamUSM64uZKhxlL0PNeCbJHBqIh9tdkEWd754C0XTC7jxfbez7b873vw/K1W+incWnrKAhgON9Hb1Ukh5VPtKnZHE8CE30h3ePTqUHtB2umSlJjLEWZ6wHzCWNJIvhNh41PO1Usq1xx90vDgUIryKdZEg6kZYBN7Vz4E9UsoHj/qvvwCXAGuCj09HeyyRQgo/1Qu20V5ymGk1s5lWM2fMalhJBU4W3FiBb9jPnjUN+BQ2bT9CXkku9//rToqrCrnp/bez7eUdkzfSKLH09MVsenmLKX2lzkyid1P4wXGDGYHFpNSBzEgPSaOJPX4wGDPaIaVcMdEBkRKHRjFDCZ8JfAbYIYTYGnztJgLG9wkhxBeBBuATk53IhpqqVd0Naaw9gSWSA0EDXHFgLuV1s958/Zg+M+wsuLkCm8vGjlvrgiUJJ2csRZucmsw9z95GwbR8vvmeG9j6ynaFd6MJhbLppRSWFbLlP1uj3ldSsQNnpp3B/eE7vwYzenGOuo6pljXWRFCjsSQSo0p4QuJJHEbdCEsp/8P4VenPiXb/kUQiqZm3I2iA57xpgI/HliJYcHMFrnwnu79Tz/Ch8INujubGR69lxqLpXP3e1doAR5kV7zgJgM2vbI16X+lzA8U5BhSM8EBmL2n9WdrwaqYmxpejJyNi4tAolgqncyCVIp1Va0Af3U4iqZuzm9byRspqZ1JWN2vM84okQeXqPFIqnGy/7wCd+4xtAPCpqz/Buz51Do/c+DNef27DuMep+oQTxZccamT1GeeeQVdrF4d3tkQ9Gjt9QQrePh8j4QZl2T0Mpw2Q11oSpZFprMkUyhOOMvEkDhPryisSMMB7aKmop6R+OtMOjh0FLVyCquvzSJ3rYuf3DtK51ZgBPvVdK/n6fV/hpaf+za/WPG7oXJrJEUKw8tyT2PCvzab0l7Ewmf49I4Q7RxzI6gEBGb3Z0RiWRhN7jPuELYOljLCqT1hVCfvxI5HUz9lLa0U9xfWVTKs+Ngr6CMIZqAedtjCJxoc7qXtNreThEWVaNbeS7/z+Fg7srOGmS25hOMprMxqYv2IuuUW5vPrM61HvK6nYQVKhk9a/HDtRC2V5uT+7GyRk9ubq5eioo9VlTIj+cnTcoO+SCZBI6ubtpq28keL6Kiqq5475pSdcgunX5ZOxNIWGH3fS/YqxKVxOfjY//PtDeNwervrQNQwrlKPUhM/bPnQGPp+P1559I+p9ZS4NVLrq2xb+vdKX00VaX5YuV6mZumgjHJ/YkUpKWGVXISn8HJy/g87SJkpqp1N+MFBF6XgFbEsSTL+hgPQFSTQ83EHXywElq+prdSW5+N7T36WgNJ8vrbqM5oYWpfNowuedH3s7W17eRl/XW7VkVP3Ctknq4GSfnMZIkwdPq3/SY4/GZ/fSn9VNacMMpXFprICq+rbU1/nEJNBytGkVs6yE3+Zj9+INdJY2UX4gUIhjLAVsT7Mx81uFAQP8o843DbAqNpvgnse/w+LTFnHzxbex841dhs6nCZ0ZC6uoml/J+qdeiXpfthRB+sJk+hTyg3uzO5E2SVZXXhRGptFozMZSUycHanWgw/EJe+0edi19g97cDir2zaOosWLM9o4sO7O+VURSqZPaB9vpfcP4tO1LP/4s53x0Fd/5xj389Y9/D6ut3g3JGO++6Fy8Xh8vPrk+6n1lLk/F5hT0KNwz3flt2Lz2MbfJ1IyPyj7kAD5DX5FWUrRTrmKWZbCUEY42o65hdi5/naG0PubuPInMlvwxj3MVOJh1SxGObDs1a1rp32G80vin7voY531lFT+5+2f8+oePGT6fJnRsNhvnffoc3nhuA93tPVHvL+f0VDzdXgb3h5c/LpF057eT3ZWPTdqjNDqNJg7QRjjxGEzrY+ey1/A63Szcdiq5nUVjRiQ78+zMurUIe5qNA3e0MHTAbbjvD69+Px+96QM895OXePDm7xs+nyY8VpxzEiWVxTx83U+j3pctRZC5PJWOF/rDTk0azOjDnTJMbs3s6AxOo4kXEsgnbCkjrJqi5J1kKaont529izdh89tZtOkM0vuz8OE7IaDLmWdn9m3F2NNtVH+nmaGa8Q1wqIFZH/7m+7nono/z0m9e4Ydf+0lIbTSR5fwvv5/ezl7+/fR/o95X9mlp2FyC7v+EHz/QWdgMfkFOR1EURqbRaGKBpYxwNGgpq+fg3J2kDqUzf+tKkkfG3iTdle9gzm0lAQN8Z8uEBjhUPvzN9/OlBy/hlSf+x0Ofexgp1fKZNeoUlObz9o+cxe8fehKPO7zKVSrkn5POyCE3Q9XhL0V3FDWT3Z2H0+OK0ug0mjghyrWj4wlLGWGBVNqMYbyNGOpm76G5spbsjgLm7FiOw+ccMwgrqcjB7FtKsKdOroCPMFmg1Kdu+DhfWnMJ6594hTs/fT9+X+B9JUoZSbMZL9Xogss+hs0meOaR56NepjK5wknanGQOresKu+1AZg+jqUOU145dr1yjmVpI4i5YLEpYyghHCp/NR/WirXQVtlDcWMn0/QsQcuxsreQyJ7NvKcHmEOz/TjPDtcYV8Be+czEXf+uTvPD4etZc8tCbBlhjLqkZKXzwsvfx3z+/RnNta9T7K3hPJv5RP13rB8Ju21Z6CJvPTn6brhetSRS0EY5LVEpQHt3G7Rpl39KNDGT2UrVvPsWNVQQ09onndc6EOTeWIH2SHd82vhsSwNcf+jIfu+p8/v6zf/LQVx/G748fA3zdddexYcMG1q9fP+4xq1atYuXKldx///3mDSxIpJXqBdd+lKy8TB6/58mInncs7Ok2cs9Ko+s/g/gGx/+bj5WP7rP56ChuIq+tGIfPOW7bcIp+HE1/DEpfDupym1FgKpXYTBwlnFDFOgbT+9hxyn8ZSh9gzvaTKGmcPm7t3fQlSSy8rQrfsJ+dt9ZFxAB/+sZP8LGrzucPD/2ZBy79YVwZYIANGzbwxBNPsGrVqjH/f9WqVTzxxBNs2DD+bk5WoaA8nwuu/Qgv/vZl9m2sjn5/787AlmSj/W99kx98HB1FTfgcXooOV0RhZBqNJpbE2/RnQiQSn0LSvR8/XQUtHFy4A7vXwYKNp5LWnzluOcvss1KZdlkeg4eG2Xr3ftw94QfsHO/bPe+TZ/Oluy/hn4/9iweujnwaUiR8yevXr+eCCy7giSee4IILLjhGER8xwMe/blUuvfdz2Gw2frb6lxMel8rYgXqTcfTkzpYkKHhvJr2bhhg97A1r0wWJpKW8jpSBdDJ7dIEOTaKglfCUQSI5NKOa6qVbSRlMZ9Ebp5PWnznu8YUfyaTyinwG946y+dt7lQzw8Sx/x1JuWbeazeu3ctcX7jN8vmhytCE+ooinmgFe+o7FnPPpVfzuvqdobWiPen9556bjyLTT+qfwVXB/Vg8DWb2UHKrSOyZpEogjRlj1xzpYSglDeD5hj8PNvoWb6SpoJb+pjKq987H57WOeQ9ih/Mt55J2dTte/B2h8pJM+X7/h8c5cPIN7/3wnjdWHWf2RWyZNg4mH6OijDfEjjzzCZZddFlEDrKouI4HD6eDKh79Kc22LKb5gW5Kg8Pws+ncMMxRmhSyApooa7B4nhU3lURidMbRfdyoRb4YrcZSw5YxwqAyk97J7yRuMJg9TuXc+hYemjask7Bk2pl9TQPqCZFr+0EPLE71jHhcuJVXFPPTsvQwPDHPN+1bT3xN+VGysWL9+PY888gi33nord9xxx5RQwACfWv1xqhZWcuP7v417xHik+2TkvycDZ7adugfCv6eGUgfoLGymvG4Wdv+U/ahqNOOgjXBcMpkSlkhaShuombsDh8fFkk1nktQ7flRt8jQn068rwJnroO777fT89638XiOqNK04me89dy+uZCdfOOsr1DbUKJ8rFqxatYrLLruMO+64g8suu4yXXnrJ0oY4lVQqFpRz8c0Xsv7x/7Djmd1RV+T2NBuF52fSu3mIwX3hq+DDVQex+W2UNEyPwug08YuVNn7QGGVK/dW8dg8H5m+nvfgw2Z0FzN15Ei5PEqOM/QWYuSKFym/k4x/yc+D2FoaqI6OM0nPTuOf5H5Ffks9Xz7uCg7utZ4CP9gG/9NJLlvcJ2x12rv7V5Qz1DfOzq35pSp/FH8/Cnmqj+bGesNsOpwzSVnyIkkOVuDxJkR+cRhPX6OXouETCuBHN/Zk97F+0mZGUISoOzKW8bhYCgR//mOq58AOZlF6cw1CNm5r72/B2q211djypWSl865/XUD6rlMvf9022v7YjIuc1i7GCsCaKmjYTI8r1wm99lFknz+Cuj36X3vbQA6RUc5OTS5zkvzuDrhcGGG0MPSL6SK7voRnV2KSNiro5Ief/6sAtzdQhcYzwlIiO7s5tY8eK/+K3+Vm06XSm1c0e/wvJDtO+kkfZZ3PpeWOI6m+3RM4AZ6Zwy3PXUrlkGtd+bDUbXtoYkfOaxURR0GNFTVuF5W9bxoXf+hgv/upl/venN0zps+yLufhGJC2/D98XPJDeS1vxIUobp+NyJ0dhdBpNvKOjo+MUOaaqzejJobiximm1s3F4T6z/fEQ92zNszPhmEZmLUmh+qpumJ7on3E4uHJ9wWlYqNz57FZXLKrj7Yw/wz2eeD7ltvLBy5coJle4RQ7xy5UrLLEtn5WZx92N30FrbxiNf/7kpfeafmUnGkmQOPdqFty+8giwSSe3s3Tg8TqbVmbNloY5y1mhih8WM8NjY/Q6mVy+Y8JiUChczry/CmW2n9odtdL0SuUjljNx0vvPct6hcXMGaTzzIG3/bFLFzm0kopSjXr19vGQMshOCuX3+bvKJcrj/zVoYHRqLepyPDzvTPFzN0YJTO58O/xzoLWujJa2fGvkU4vOOXqNRopjaJsxxtOSOsUjs6/XQXsy4rxTvkY+dtdQwcjNyXcU5xNt957luUzi7mzvPvY9OzWw2db7Ldl+IJVX9ptHcrOsKlt3yRt73vTO786hoOb2o1pd/pny/Gnmqn5pH2CVdZxsJn81E7exepAxmUHqqKyvg0ZqKjnI2hjbD1sUPJRVnkfyCDvr1D7HvwEJ6eyP1hi6YXcufz3yK7KJvb37+GbS/ujNi5NcZ454dXcdntl/LXXz3DH376R2azMOp95p2eScHbsmj4fRsjjeFXWjs84wAjqUMs3njGuLt6aTSJgcSuUKL4CJGJ8jEHSxnhiaKjj8eRY6fyyjzSFyTT+I9Wqn/ViPSFJ00mUqUzFldx7z/vwOFycM3ZN7F3w/6wzq0Jn1CV7KwlM7nr17ez+409PPCV7xtSwEmEFhjlzLMz88slDOwfof1PA7hwhdXPYEYvzRV1FB+uJKenQGWomqhgpZ2JrDTWibEBqQp7xx/BeK1D84ivKx8hMpYmU3FFHsIlqP9BBwf+0xDR889dMZv7n/8OwwMjXHXWDdTvaYzo+TXqFJTl88Df72GgZ4DVH7mVUROqYmGHGVcWIeyCuh+1Ee53h1/4qF64DafHNWlsg0ajmVpYzghP5BMWdii+MJuiD2cx3OCm7sF2Rpu8yn7WsaKjK+dVcM8/vk1vVx+Xv/NqWhrG3gw+HmpAJxrpWWk88Pc1pGWm8pW3fYP2pg5T+i37ZC7p85Kp+V4roy3hq5HGmdUMpw8wb8sKnN7wFPTRmJ0nPDQ1Mhw1cYgNSZpC/M8RtBKOAa5iB1XfyCd1VhId/+rn8C+7kW71P+JYFFcU8YPnv4vX4+XK864b1wBrzCcpJYn7/3o3VfMruOb9N3JwhzlVynJOT6P4/Gza/tlL96vhT7y689ppqqqh6FAFOZ2FURihRmM97GDICFsJixnhsfOE0xckM2N1IdIrqXmgjd7XIx9hnFOYw/eeu4/ktGS+9o6rOFzTFPE+rIZqBatIR1U7k5zc8dTNLD5jId+58D52/mvvCceq9jlRtaqUKheVlxUwsHeEpl92h1zZ6gijSSMcWLiNlIF0KvfPVxpfJLCSotU5zdEiviKRjfqErUTUjbAQ4hfAB4A2KeWi4Gu5wO+BKqAOuEBK2a3ax3C9m94NQzQ93o2nM/JxcRnZ6Xz/ufspLC/gyvOuM01laSYnMzeDO/78LZactZDvfukH/Pup/5rSrzPXzozVhfgG/NQ+0I4M87bzCx/7l2zGb/MxZ/tJ2P326AxUo7EgNqKvhM2wTaFghhL+JfAj4FdHvbYaeEFKuUYIsTr4/IZQTjZWdLR/0E/tD9vGbaPqnx1mkNT0VL77j7upnDeNb3zgat74nzllD83CrJzdaFA+u5S7/norRZWF3HHhvax/4hVT+rWn25h5UxH2ZBv7b2nG23uiBZ7IPyuR1M7bzUBWD3O3nUTaUEZI7aYSVlLfminLL4mgbVIl6kZYSvlvIUTVcS+fD6wK/r4OWE+U36gKySlJ/OCvD7BgxXyu+8RNvPavqWWArcyK85Zzy+9vwOf1cc05N7Pr1T2m9GtPtTHr5iKSip0cXNOqlg9cWUNbWSPltbPIay+Jwig1GmtjNDArFOLFNsXKJ1wkpWwGkFI2CyHGjUgRQlwKXAqQQRajhF/tSkUJO1wOHvjjfZz09uXcdNGtvPTn9VHvM1FQ9SUfaXf+Ve/nC9/9DA27GvnOh+6jrb590nOq9nm0MrWlCGbeVERypYu6+9sZ3DkatnLtKGqiYfZe8ltKmXZwjtKYxsOr1aVmihBDn3DItilSxH1glpRyLbAWoEiUmhIu53DaufHJqzn1PSv49hfv5NnfPWdGt5pJSElP5opHv8rbLzyDV//4Og9+9keMDI69V3SksafZmHlzISlVLuoeaqdvy3DY5+jOa6N64VYyu3OZtXtJwiw9azThEgGfcL4Q4uht7NYGbUncESsj3CqEKAnONEqA8R26JmOz27j2sW9w6odWcPfX7uXPv/hLrIcUErHw7ZrZ58wFM3jgybspm1vK/13/G56637y/iyPLxsybi0gqdVJ7f2gG+HgD25vdyb4lm0gdyGD+1pXY/WN/9MKNsD6axIgl1WhCokNKuUKhnem2KVbrV38BLgn+fgnwdIzGcQw2m+DqdV/nbZ84nUevXscTjzwV6yFpgI99+cM8tvGXZOSlc8t53zHVACeVOJh9VwmuYgc197YpKeDe7E52L3uDpJFUFmw5FYdP746k0UyEPegTVv0xgOm2yYwUpd8ScHTnCyEOAbcBa4AnhBBfBBqAT4RyLj9+JV9rKBWzhBBc+7NvsOqis/jZ6l/y24f+YKkdjVRR9ZUaIVQFnZmTwfU/uZpzLljF689t4NHP/oae1j4lBa7SJmNOCrOvL0ZKOPjtNoZr3GEvIffkdLBn2QaShlNYtPk0XJ6ksMcRKqo5tGa302gmwwyfcCRtkxHMiI7+1Dj/dU60+w6Hb/zoq7z3C+ex7tuP89t7/xDr4SQ8p717JTf94npyCrL58eq1/Oa+31EmK03rv+AdWcy8tARPh4+au9txt4ZfzKCjsJn9i7aQMpTGws2n4XJHzwBrNFMJM/KE48U2xX1g1tH48UdFnX71u1/k/K+9n9/e+wfW3f54xM8/VYmGTzg1I5XL7/sKH/nqhzi4s5ZrP3AT+7dUGz5vqLshYYOyi3Ip/mA2fduHaHioC99g+DPy1rIGDszbTmZvDgu2nhpyTWgdrKXR+xCbY4TjhanzV1Pk4psv5IJrPsKffvhXfrb6l7EeTkJz9sffwZUPfZ380jwe/+7v+em3fo57NPw8XFUc2XZmXFVIxoIU2p7tpXFdJw5feB8RiaRhzl5aKurJbS9i3o6Txw3C0mg0moT+dnj/l97NF+78DM//+kV+dGV8RK+b7aONh6jqFect5wt3fYZ5K+dwcFst3/7YPex9Yz92XKQctS9vNMeasSiZ6d8oxJZso/YHbXT9ZwAIL1rZ6/BQvWgrPfntlDXMYEb1QoSM/9zdRPDt+tBlQa2EQOJIkHj/hDXC0xdVcsUPv8Ibz27ivi98HykTY+kjnlh+9lI+d/unWfy2hbTWt7Hmkgf512/W4/eb9+ETTkHZp3Io+kA2w4fc1NzezMjh8NX3YEYv+xZvwZ08zIw9iyg/PDMKo9VoEoeJtq2dSljKCEvF6Ojj2zhdTm587GoGegf49mfvZMDbF6khJhSqqv30s0/hU7d9nEVvX0B7YwcPf/VnPP9/L+F1e0mewHdrZJVgLEWbOtNFxWX5pFS4aH+2j8O/CWx/efSxk/loJZLWskZq5+zC6XGxaNPpZPTmKI/TiE/YbEVrpL9EUN8adSTaCE9pLrruQmYtmck177+R7vaeWA8nIRBCsOpDb+cLN17C4lMX0Xm4i0e+/nP++egLeN3mbqNmT7NR8sls8s/LwNPj4+A9rUr5v27nKAcX7KC7oJWsznzm7FyGM4opSBpN4iDH3KxnKpJwRri4sojP3XwxLz65nlefeS1q/aj6L1VrTkd6j95ItXU4HZz3ybO5+IZPMWNhFYcOHuanl/6S9ev+i9ftxUkSTkIzXEbGKhAgIOftaZRenIMjw0b7P/ppeaIH/7AMW4F25bdycP52vA4vVfvnU9Iw/Zhz6Cjn6KB3X9JMNRLOCH/mhk+BEHz/6h/HeihTmsycDD78lQ/y8Ss+QkFpPgd31nLbp+/khSdeosRXYfp40uYnUXZJDqkzkhjcP0rNnZ0M14fv+x1NGqF27i66CltI7c9gwebTSBvMmLyhRqMJC70cHYeoVsw6QkZ2Ou/97Lt47rF/0XaoPYIjixxmRysb8bMeP1abzcaity3gnE+/g3MvficpaclsfG4z3/3CD9j43BaklLhIMTrk8MZY6qLi04XknZKJu8NL/Q866PnvEMjQ1OqRY6Tw01xeT8PMfUjhp/LAPErrZ2AbJ/pZtQa0kZ2QVFWiVpeaeEQb4SnI6e87jZS0FP70E2tsymAFbDYbi89ayDs+cSZv/9iZ5BbnMDw4wstPvMKTD/6Z2p31MRlXcpGT8o8XUHBWFv5RP82/7aH9b/1IT/gf7O68Nmpn72Y4fYDsznxm7F1MyrD5qV0aTaIgQfuEpyKzlszAPepm/9YDsR5KxDHTJ2x32Fl21hLe+/F3c8ZHTyWnOJuRoVE2/n0z/3nyNTY+s/nNLQbHUtqqYw2l6pWryEHJR7LJe0cG0itp/VsvrU/3YOt3IIL/QmUorZ+G2fvozm8jeSiV+dtWktdeHFV/r5UijhNDQSfUV6QmBiTUHZZfkkdvZx8+ry/WQ7EcyanJnPquFbzjI2dx5gdOIzM3k5HBETb8fQv/efJ/bHxmC6ND5uztOxapM1wUfSibnNPSkF5J27N9tPy5B29v4G8dWtHIAEOpAxyaUU1nUTN2r4Pp+xdQ2jgdm9QFHzQac5B6OXoqcmD7Qd772XeRU5AdUmpSLHYYUkV1rBOp0pzCbE5578m87SOns+Jdy0lKSaKvq5///fUN/vOn/1HzXCPuYTcANhykhHg7qY71eD+rcEDWKWkUvCeD9HnJeAf9tD7dS/s/+vH2+I5pE4p6HU4d4ND0A7QXH8bms1NWN5OKhjk4PeGYcELuL16IRc6uzhOOBkZS/eLPFGgjPAXZ+u/tAHz+ls/w4Dd+GOPRxBd2h50ZS6pYcNo8Fpw2l/mnzaN8dikArQ1t/P1n/+S/f36N7a/senMlIY+CmIw1eZqTvLMzyD0rDUemndEWD4fWddH5Yj/+4fA+uBJJf3Y3hytr6C5oxeazUVo/g7L6GTg9STgS6yOi0cQF2ic8Rdm9YS+Pf/f3fPraC6nf28hTP/5z1PqK9xrQRdMKmb9yHitOO5l5p81m1oqZJKUEFF9Xczd7/7ef5x99ia3Pb+fglto32yUdldNrph/akWYj//QMcs9OJ21WEn6PpHfDEF0vDdC/fYQjk+bxFOjxr/uFj87CFpoqahnI6sHhdjKtZjbFjVXH7PmrqmhV21nJJ6zRRBOthKcoD9+wlsp5FVz78JWc8q4VPHz9T2nY3xjrYUWV7Pws5q+cx/yVc1kQfMwtygXAM+rhwKYa/vHIc+x9rZp9r+2nvbEzxiMOYEu2kbsyg/wzMslemo7NIRhucHPo/7ro/s8gvv7wZ8qD6X20ljbQXnwYr8tD8mAaM/YuorBpGna/9vlqNBpzSQgjfIzy8sMdH7uXj3/zw1x08wX8ft+v2Lexmn8/9Sqb/7WV+j2NjAyOxG6wCgghKCsup3h6EaUzSyidWUzZrNLgYwlZ+VkA+P1+GvY08sY/NrNvQzX7NlTTva3/hLKRoSpVdT/0+O1sSYLMZSlkn5lG1kkp2Fw23B1e2p/po++/I4zUHimwIbCHuDOO1+GhtaiBltIGBrJ6EH4beW3FFDdVkN1VMKFqNTvf10jEsWpbVQWtlbcmeujArCmNx+3lt/f+gWd/+S/e9dlzOOujp/Pley6Bey4BoLW+jbrdDRzcXUPd7npqd9dRt6eegV71QiFGsDvsZGSnk1+aT+n0YkpnlFI2o4SS6SWBx6piklLeWkL1+Xy0NbTTdLCFl//wKoerm9i3sZrqzQcZHji2RnKs/LpHkzzNScaSFDKXppC+IBmbS+Dp8dH5wgDdrw4yuH8UJGH5Z712D92FrXQUNtOb14G0SVIHMpixbyGFLeW6xrNGE+don3AcYsOm5E8cYmjM17tbe/j9/U/x+/ufIr8sj3kr51C5YBqVCyqoXDCNZavOP8a4dRzupH5PI/W7Gzi0v4mRoVF8Hi9ejw+vx4vP48Xn9eMNvubz+vD7/G89+vz4fX6S05LIyEknIyed9OBjRm7GCa+l56SRkZNOasaJynGgZ4DmmhYadh/i9b9tpLO2m5aaVlpq2mira8PrOTENSyBOUK/Rrh09Fs4sO+mLkslcmkLGkmScuYHbcOSQm47n++nbOMTAnlGOfAaPKNXJ/Kxeu4fu/DY6i5rpzm9H2vy4RpIpaayiqHUa6X3ZYftqzfbtxkJdxiLfV7VPvS9wYuDHzyjWWpFUxVJGOJp0HO7kP4f/x3/+/L83X7PZbBRVFgSNcsAwVy2o4D2fP3dMw2iE4cERBroH6O8eYKB7gJa6Vvq3DNDf9dZr3a3dNNW00FzTwkDPsao8ntOpkoqcZM1PJWN+KpnzUkkpCUxsvP0++reP0L+9h/7tI3g6w8vflkiG0vvoyeugO6+d/uwupE3iGkmm+FAl+a0lpPcGDK+OctZoNPGI/maaAL/fT3NtK821rbz29w3H/F9OYTbOJCcOpx2704HDacfhdGAPPjqcduwOOza7DZvdht1hxx78fWRwlP6jDO5AzwAek7fzixbCIUgtTyJjbgqZQaPrynUC4On30r93mK7nBxnYM8pwjZtw3T5u1wi9uR305HXQk9uBJylYmas/g9L6GeR0FJLRm2OpPF2NRnMsRvcJsBIJYYSjoRJH29yM4o7IuZy4cAZrOpm9JaERn3BmTgYplS5SKl0kVzpJqXCRXOZE2AMG0N3hZWD3CAN7+hjcM8LIYc8xvt3JAqskkuG0Afqzu+nL7qI/q4fR1IBrwel2kdWZT05XAVmdBSS5Jy5paZVUIyulKFkrMCshvuqmDH7847oRpxr6ztRMii3ZRmqZi9SKZFIrk0gLPjoz3rp93O1ehhvc9G4cYrjezdABN+728NS9z+5lMKOP/uwu+rK76c/qxucMnMM5mkRmTw4lh6rI6s4jrT9Tq12NZooitRKOT4RiYJaVMFvRppAKApzZdlyFjjd/kgqCj6VOXLlv3Sa+ET8jDW56XxtmpMHDSL2H4QY3/qET15XHTe8RkuGUIYbS+xhK72covZ/B9P43VS5AykA6+a0lZPTmktmTQ9JwKk6cSu8RJlfd49Gr2C4W2wpaKRhMo5kIvRytmVoIcGbaSSpwkVTgJLnQSVKhk6QCFymFSbgKHNicx34Re7q9uNu99O8YYbTJw8hhD8P1btxt3jf9uJPl0EokHtcoQ+kDDAeN7VD6AMNpA0h7MPTZL0gZSiO9L4uipnJSBzLJ6M1Rqtes0Wg0VsNSRtiBQ0nxDcdx5PDxhKSEBThSbTgyHTgz7TgzHaRmpeDItAd/bDiP+t2RaUc4jjOyfT7cbV5G6t30bRjG3RYwuoFH37j77h5teO3Y8dm9jCQPMZoyzGjKECPJw4ymDDOSMsRIyiB+x1sRz87RJNIGMshpnE7aYCZpA5mkDqZjC7FSlaqaBfN9u20mF84A84t1GFHtOtUoGljq63xCtE9YE3VsLhuOFBv2VDuOVAfOVDuOVDspqak40gKv21NtONKCjyl2HBn2gHHNcGBzjP3F6Rvy4+n14e3zMdruZfBg4HdPlw93h/dNY+sfCRjZidSsROJzeHEnjzCaHDCyo8kjAYObPIw7eRivy3NMG5vPRtJICslDaWR155E8lErqYDqpA5lvqlsjxlSj0Ux9/Pj0cnQ8YsM25sbu9jQbFV/MR/ok0g/SL8HHic/9Ehl8xB/4fyQgAj/CJgKPAjjm90BpSN78PfD/NodAOAU2Z+BROAj87hDYnLbA/zkI/J/zqOOTRKCvCZB+iX9Y4hv04xsK/Hjb/Iwc8ODt8+Pt9QUe+3zB3yW+fh9yglgoicRn9+J1ufFkufE6PfhcHjxONx6XO/B68HeP043XNYrffmzVGpvPRtJwKkkjKWT15gYM7nAqySOpAb+txzWp6lTN2VUtIQnqvt12k33CVoqOjg2W+srSKKJ9whZDOAWpM5MQ9qAhtYGwgbCLQLrM8c9DRPoDRlpK4JjfQUqJ9IL0SPxeGXj0BB6lJ2A8/R6JPPr/jvw+KvEN+/GPyDcN7BFj6x8KvOYfkePm0PqFH5/Dg9fpDT568Lm8+Iq9eI88d3iDj4HnAaPqRtrGWWb22XG6XTjcLpxuF6kD6Tg9LpyjySSNJJM0kkrScMoxRlYrWo1GEw30cnTcMk7lox7Yd2VzOKdB2AGbQNhA+glYV3/QyB75CSKRIALlxKUI/i4kUvjx2yTS5kcKP9Im8Qs/0uYHG8HX/PhF8Jjg7367D7/Nh9/mD/4efMz04c/xB//Ph9/uP/Y4uw+f3TuuIT2CzWfD4XFi9zpxeJ0kjaSS3pcdMKpuF05P0psGN9mTgtPtwu4P/1ZQNcKqirbJwO2qqmhV26n6hFX7A/N9u8b8uhb76jEVfW0SiSnx13a7Rth1ymtKbd80sEIiBUFjK48xttFc5RN+GzafDZvfjs1nx+a3YfPZsfvtOLwObL6k4P8FjrF7Hdh9Duw+e8DQegKG1ulx4fA6cXic2GToX45azWo0mnhDL0ebhBDiPcD3ATvwqJRyzUTH2xC4ODF1Rfghp7NQfRxSIAIWOPj7kddsx70ePE4KBAKb3xYwolIEHv22QBu/DduRR3/gPEf/nz1oaG1++wn+U9UoXlVjaqTghdl+VlV1aaRP1XaxyNk1X9FOiTl8FNHXRxUzjHC49idaxOwuEULYgYeB84BDwAYhxF+klLvDPZfT62LOnmURHqExA6XRaDQaNaLtE46k/TFKLKdqpwAHpJQ1AEKI3wHnA+NeBC82OiyS82uVDdaN5HqqKlMrRQ5ba8u9RFBeifAeNSYQtv2JFrG8o8uAxqOeHwJOjdFYNBqNRhMnmFA7Om7sTyyN8FjS5ISwXyHEpcClwaejP+RLO6M6KuuSD3TEehBxjL4+46OvzfjoazMxc6N03n8SuPaqJAshNh71fK2Ucu1Rz0OyP2YQSyN8CJh21PNyoOn4g4IXbi2AEGKjlHKFOcOzFvraTIy+PuOjr8346GszMccZuoghpXxPNM57FCHZHzNQdwgaZwMwWwgxXQjhAj4J/CWG49FoNBpNYhA39idmSlhK6RVCXE5g2cEO/EJKuStW49FoNBpNYhBP9iemoYZSymeAZ8JosnbyQxIWfW0mRl+f8dHXZnz0tZkYy14fBfsTFYSUMfFFazQajUaT8MTSJ6zRaDQaTUITd0ZYCPELIUSbEGLMVCQhxCohRK8QYmvw51azxxgrhBDThBAvCSH2CCF2CSGuHOMYIYT4gRDigBBiuxDipFiM1WxCvDaJfO8kCyHeEEJsC16f28c4JlHvnVCuTcLeOxCoMCWE2CKE+NsY/5eQ902kiMfyM78EfgT8aoJjXpFSfsCc4cQVXuAaKeVmIUQGsEkI8fxxpdbeC8wO/pwKPEJiFEEJ5dpA4t47o8DZUsoBIYQT+I8Q4h9SyqN3PknUeyeUawOJe+8AXAnsATLH+L9EvW8iQtwpYSnlv4GuWI8jHpFSNkspNwd/7yfwoSg77rDzgV/JAK8B2UKIEpOHajohXpuEJXg/DASfOoM/xweEJOq9E8q1SViEEOXA+4FHxzkkIe+bSBF3RjhETg8uHf1DCLEw1oOJBUKIKmA58Ppx/zVWObaEMkYTXBtI4HsnuKS4FWgDnpdS6nsnSAjXBhL33vkecD3gH+f/E/a+iQRWNMKbgUop5VLgh8CfYzsc8xFCpANPAVdJKfuO/+8xmiTMrH6Sa5PQ946U0ielXEagOtApQohFxx2SsPdOCNcmIe8dIcQHgDYp5aaJDhvjtYS4byKB5YywlLLvyNJRMM/LKYQwUmPUUgR9Vk8Bj0kp/zjGIXFTjs1sJrs2iX7vHEFK2QOsB44vDZiw984Rxrs2CXzvnAl8SAhRB/wOOFsI8Zvjjkn4+8YIljPCQohiIYQI/n4KgffQGdtRmUPwff8c2COlfHCcw/4CfDYYsXga0CulbDZtkDEilGuT4PdOgRAiO/h7CnAusPe4wxL13pn02iTqvSOlvFFKWS6lrCJQ2vFFKeXFxx2WkPdNpIi76GghxG+BVUC+EOIQcBuBQAmklD8BPg5cJoTwAsPAJ2XiVBw5E/gMsCPovwK4CaiAN6/PM8D7gAPAEPB584cZE0K5Nol875QA60RgM3Mb8ISU8m9CiK9Cwt87oVybRL53TkDfN5FDV8zSaDQajSZGWG45WqPRaDSaqYI2whqNRqPRxAhthDUajUajiRHaCGs0Go1GEyO0EdZoNBqNJkZoI6zRaDQaTYzQRlij0Wg0mhihjbBGYwJCiNlCiDohxKzgc2dwM4DyWI9No9HEDm2ENRoTkFJWA2uBdwdfuhx4Wkp5KHaj0mg0sSbuylZqNFOYncC5Qohc4Ivojc81moRHK2GNxjz2A3OBbwPflVIOxnY4Go0m1uja0RqNSQS3WmwCDgJnSCnH2yRdo9EkCFoJazQmIaX0AH3Aam2ANRoNaCOs0ZiNE3g51oPQaDTxgTbCGo1JCCGqgPpE3odWo9Eci/YJazQajUYTI7QS1mg0Go0mRmgjrNFoNBpNjNBGWKPRaDSaGKGNsEaj0Wg0MUIbYY1Go9FoYoQ2whqNRqPRxAhthDUajUajiRHaCGs0Go1GEyP+H9CfO/arAd72AAAAAElFTkSuQmCC",
      "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, delta_ts, cmap='nipy_spectral')\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label(r'$\\Delta$TS')\n",
    "plt.contour(gamma_vals, ns_vals, delta_ts, [chi2_68_quantile], colors='#FFFFFF')\n",
    "plt.contour(gamma_vals, ns_vals, delta_ts, [chi2_90_quantile], colors='#AAAAAA')\n",
    "plt.contour(gamma_vals, ns_vals, delta_ts, [chi2_95_quantile], colors='#444444')\n",
    "plt.plot(gamma_best, ns_best, marker='x', color='white', 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)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the significance (local p-value)\n",
    "---"
   ]
  },
  {
   "attachments": {},
   "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": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.utils.analysis import create_trial_data_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function create_trial_data_file in module skyllh.core.utils.analysis:\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 : instance of 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 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)"
   ]
  },
  {
   "attachments": {},
   "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": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.timing import TimeLord\n",
    "tl = TimeLord()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10001/10001 [08:52<00:00, 18.78it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "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.006 sec/iter (10000)\n",
      "[Generating background events for data set 4.]   0.019 sec/iter (10000)\n",
      "[Generating pseudo data.                     ]   0.027 sec/iter (10000)\n",
      "[Initializing trial.                         ]   0.030 sec/iter (10000)\n",
      "[Get sig probability densities and grads.    ] 4.4e-06 sec/iter (1950580)\n",
      "[Get bkg probability densities and grads.    ] 3.3e-06 sec/iter (1950580)\n",
      "[Calculate PDF ratios.                       ] 9.9e-05 sec/iter (1950580)\n",
      "[Calc pdfratio value Ri                      ] 5.4e-04 sec/iter (975290)\n",
      "[Calc logLamds and grads                     ] 2.7e-04 sec/iter (975290)\n",
      "[Evaluate llh-ratio function.                ]   0.003 sec/iter (195058)\n",
      "[Minimize -llhratio function.                ]   0.058 sec/iter (10000)\n",
      "[Maximizing LLH ratio function.              ]   0.058 sec/iter (10000)\n",
      "[Calculating test statistic.                 ] 5.1e-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)"
   ]
  },
  {
   "attachments": {},
   "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": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGwCAYAAACgi8/jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGUElEQVR4nO3de1xUdf7H8fdwFVQQNG4qaLVSmOkGYqy2aVGolavdNPsVWdlWaLa01Vqmbqnd1nIzdt0uanupdd1da7OLm2hZLV7QrFzMsiUtE7wiCepwOb8/aCZmuM3ADANnXs/HYx5yLvM9n2E8Mx++V4thGIYAAABMKMDXAQAAAHgLiQ4AADAtEh0AAGBaJDoAAMC0SHQAAIBpkegAAADTItEBAACmFeTrAHyttrZW3377rbp37y6LxeLrcAAAgAsMw9B3332nhIQEBQQ0XW/j94nOt99+q759+/o6DAAA0Apff/21+vTp0+Rxv0108vLylJeXp+rqakl1v6iIiAgfRwUAAFxRXl6uvn37qnv37s2eZ/H3JSDKy8sVGRmpY8eOkegAANBJuPr9TWdkAABgWiQ6AADAtEh0AACAafltZ2QA/qGmpkZVVVW+DgOAm4KDgxUYGNjmcvw20bGNuqqpqfF1KAC8wDAMlZSUqKyszNehAGilHj16KC4urk3z3DHqilFXgCnt379fZWVliomJUXh4OBOCAp2IYRiqrKzUgQMH1KNHD8XHxzc4x9Xvb7+t0QFgXjU1NfYkp2fPnr4OB0ArhIWFSZIOHDigmJiYVjdj0RkZgOnY+uSEh4f7OBIAbWG7h9vSz45EB4Bp0VwFdG6euIdJdAAAgGmR6AAAANPy20QnLy9PKSkpGjp0qK9DAQC7kSNH6u677/Za+TfddJPGjx/vtfJ94auvvpLFYtH27dt9HQo6IL9NdHJyclRUVKQtW7Z4vGzDMFRprValtVp+PnofAACfYni5F5yoqlHK7DWSpLSkKK28PYNOkQBMy2q1KiQkxNdhAI3y2xqd9lK456hOVDH7MuBr9Wta2/vhbs1udXW1pk2bpsjISPXq1UsPPfSQvYw//elPSktLU/fu3RUXF6fJkyfrwIEDkqSaWkOffFOmf+YX6LLLLldERIS6d++uCy64QF9++WWj19qyZYtOO+00Pf744/Z98+bNU0xMjLp3765bb71Vv/rVrzRkyBD7cVvz1/z585WQkKDk5GRJ0qeffqqLLrpIYWFh6tmzp2677TYdP37c/rzGmuXGjx+vm266yb7dr18/LViwQDfffLO6d++uxMREPffccw7P2bx5s3784x+rS5cuSktL00cffeTW7xf+hRodLwgLDlThrEylzVvr61AAfK9+TWt7K3o4S+Ehrn/cvvTSS7rlllu0efNmFRYW6rbbblNiYqKmTp2qqqoqPfLII0pOTtaBAweUm5urm266SW+++aYkqXT/t7r56st00ahRWrdunSIiIvThhx+qurq6wXXWrVunK6+8Uk888YRuu+02SdJf/vIXzZ8/X7/73e80fPhw/fWvf9XChQvVv39/h+fm5+crIiJC77zzjiSpoqJCWVlZysjI0JYtW3TgwAHdeuutmjZtmpYvX+7W72vhwoV65JFH9MADD+jvf/+77rjjDl144YVKTk7W8ePHdfnll+uSSy7Rn//8ZxUXF2vGjBlulQ//QqLjBRaLReEhbV+IDIB/6tu3r55++mlZLBYlJyfr008/1dNPP62pU6fq5ptvtp93+umn65lnntHQoUN1/PhxhYV31YqXXlC3iAi9/Mor6hJa15w0YMCABtdYtWqVbrzxRr3wwguaOHGiff/ixYt1yy23aMqUKZKk2bNn69///rdDzYwkde3aVS+88IK9yer555/XyZMn9cc//lFdu3aVJD377LO64oor9Pjjjys2Ntbl1z927FjdeeedkqT7779fTz/9tNavX6/k5GS9/PLLqq2t1YsvvqguXbpo4MCB+uabb3THHXe4XD78C4kOAL8QFhyoooezfHZtd5x//vkO/foyMjK0cOFC1dTUaPv27Zo7d64+/vhjHT16VLW1tZKkvXv3Kvmss7Wr6FOdl56h4ODgJsvftGmTVq9erb///e8NRmDt2rXLnmTYpKena926dQ77Bg0a5NAvZ+fOnRo8eLA9yZGk4cOHq7a2Vrt27XIr0Tn33HPtP1ssFsXFxdmb53bu3Klzzz1XXbp0sZ+TkZHhctnwPyQ6APxCXU1r5/7IO3nypLKyspSVlaW//OUvOu2007R3715lZWXJarVKkkK7hLVYzhlnnKGePXtq6dKluuyyy5pNippSP6FxVUBAQIP+So1N7e8cj8VisSd0gLvojAwAHcymTZsctjdu3Kgf/ehH+uyzz3T48GE99thjuuCCC3TWWWfZazpsBpw9UNs2FzS7NlCvXr20bt067d69W9dee63DucnJyQ2m3XBlGo6zzz5bH3/8sSoqKuz7PvzwQwUEBNg7K5922mnav3+//XhNTY127NjRYtnO1/nkk0908uRJ+76NGze6VQb8i98mOkwYCKCj2rt3r3Jzc7Vr1y698sorWrx4sWbMmKHExESFhIRo8eLF+t///qd//etfeuSRRxyeO+mmqar47jtNvu46FRYW6osvvtCf/vQn7dq1y+G8mJgYrVu3Tp999pmuu+46e2fl6dOn68UXX9RLL72kL774QvPmzdMnn3zS4hQZ119/vbp06aLs7Gzt2LFD69ev1/Tp03XDDTfYm60uuugivfHGG3rjjTf02Wef6Y477lBZWZlbv5vJkyfLYrFo6tSpKioq0ptvvqnf/OY3bpUB/+K3iY43JwwEgLa48cYbdeLECaWnpysnJ0czZszQbbfdptNOO03Lly/XypUrlZKSoscee6zBl3yPqGg9v+I1HT9+XBdeeKFSU1P1/PPPN9o8FRcXp3Xr1unTTz/V9ddfr5qaGl1//fWaOXOmfvnLX+q8885TcXGxbrrpJoc+MY0JDw/XmjVrdOTIEQ0dOlRXX321Lr74Yj377LP2c26++WZlZ2frxhtv1IUXXqjTTz9do0aNcut3061bN73++uv69NNP9eMf/1gPPvigw9B4wJnF8POpe8vLyxUZGaljx44pIiLCY+VWWqvtQ1ndHVoKoG1Onjyp4uJi9e/fv8UvaDOpqTX032+PSZIGJkQqMMAzE5VecskliouL05/+9CePlAe4qrl72dXvb759AQB2lZWVWrJkibKyshQYGKhXXnlFa9eutc+XA3Q2JDoAADuLxaI333xT8+fP18mTJ5WcnKx//OMfyszM9HVoQKuQ6AAA7MLCwrR2LbO6wzz8tjMyAAAwPxIdAABgWiQ6AADAtEh0AACAaZHoAAAA0yLRAQAApuW3iQ5rXQEwu5/+9Kd6+eWXfR0GXDRp0iQtXLjQ12GYjt8mOqx1BaCjsVgszT7mzp0rSVq1apXOP/98RUZGqnv37ho4cKDuvvtuh7Je/9e/VFpaqkmTJundd99tsex3331X999/v/r166fvvvvOoawrrrhCP/3pT1VbWytJ+vjjjzVu3DjFxMSoS5cu6tevnyZOnNhgJfX6DMPQ7NmzFR8fr7CwMGVmZuqLL75wOKdfv34N4nrssccczvnkk090wQUXqEuXLurbt6+eeOKJBtcqKytTTk6O4uPjFRoaqgEDBujNN9909W1o4MiRI7r++usVERGhHj166JZbbtHx48ftx7/66qtGf6fOq6q3FNesWbM0f/58HTt2rNWxbtiwQVdccYUSEhJksVj06quvNjhn7ty5Ouuss9S1a1dFRUUpMzNTmzZtanO59d1+++2yWCxatGiRw35X3mNP89tEBwA6mv3799sfixYtUkREhMO+X/7yl8rPz9fEiRN11VVXafPmzdq6davmz5+vqqoqh7IWP7tYU6ZMUUBAgH7yk584lHPttddq9OjRDvt+8pOf6OGHH1a3bt2Um5trL2fp0qVav369li1bpoCAAB08eFAXX3yxoqOjtWbNGu3cuVPLli1TQkKCKioqmnxtTzzxhJ555hktWbJEmzZtUteuXZWVlaWTJ086nPfwww87xDV9+nT7sfLycl166aVKSkrS1q1b9eSTT2ru3Ll67rnn7OdYrVZdcskl+uqrr/T3v/9du3bt0vPPP6/evXs3GdvIkSO1fPnyJo9ff/31+u9//6t33nlHq1ev1oYNG3Tbbbc1OG/t2rUOsaemproV1znnnKMzzjhDf/7zn5uMpSUVFRUaPHiw8vLymjxnwIABevbZZ/Xpp5/qgw8+UL9+/XTppZfq4MGDbSrXZtWqVdq4caMSEhIaPd7ce+wVhp87duyYIck4duyYR8utOFVlJN2/2ki6f7VRcarKo2UDaN6JEyeMoqIi48SJEw2OVZyqatdHay1btsyIjIxssH/GjBnGyJEjG31OdU2t8fHXR431278wLBaLsWPHjkbPy87ONn72s581eqywsNAIDg423nrrLWPPnj1GRESEkZeXZz++atUqIygoyKiqcv211dbWGnFxccaTTz5p31dWVmaEhoYar7zyin1fUlKS8fTTTzdZzu9+9zsjKirKOHXqlH3f/fffbyQnJ9u3f//73xunn366YbVaXY7vwgsvNJYtW9bosaKiIkOSsWXLFvu+t956y7BYLMa+ffsMwzCM4uJiQ5Lx0UcfNXkNV+P69a9/bYwYMcLl2JsjyVi1alWL59m+B9euXdvmcr/55hujd+/exo4dOxp9P1t6j501dy+7+v3NEhAA/ErK7DXter2vHrvMo+XFxcXp5Zdf1o4dO3TOOec0es5HmzcqPDxcZ599ttvlp6amaubMmbr11lt1xhlnKD09XXfccYfD9aurq7Vq1SpdffXVslhaXiG9uLhYJSUlDutlRUZGatiwYSooKNCkSZPs+x977DE98sgjSkxM1OTJk/WLX/xCQUF1X1UFBQX66U9/qpCQEPv5WVlZevzxx3X06FFFRUXpX//6lzIyMpSTk6PXXntNp512miZPnqz7779fgYGBbv8+CgoK1KNHD6Wlpdn3ZWZmKiAgQJs2bdKECRPs+8eNG6eTJ09qwIABuu+++zRu3Dj7MVfjSk9P1/z583Xq1CmFhoZq7969SklJaTbGBx54QA888IDbr02qq2l67rnnFBkZqcGDB7eqDJva2lrdcMMNuvfeezVw4MAmz2vuPfYGEh0A6ESmT5+u999/X4MGDVJSUpLOP/98XXrppbr++usVFFyXAOzf97ViY2MVENC63gmzZs3SsmXLtGnTJn3++ecOycz555+vBx54QJMnT9btt9+u9PR0XXTRRbrxxhsVGxvbaHklJSWS1OB4bGys/Zgk3XXXXTrvvPMUHR2t//znP5o5c6b279+vp556yl5O//79G5RhOxYVFaX//e9/Wrduna6//nq9+eab2r17t+68805VVVVpzpw5bv8uSkpKFBMT47AvKChI0dHR9ti7deumhQsXavjw4QoICNA//vEPjR8/Xq+++qo92XE1roSEBFmtVpWUlCgpKUkJCQnavn17szFGR0e7/bpWr16tSZMmqbKyUvHx8XrnnXfUq1cvt8up7/HHH1dQUJDuuuuuJs9p6T32BhIdAH6l6OEsX4fQJl27dtUbb7yhL7/8UuvXr9fGjRt1zz336Le//a0++PA/kqRTJ0+oS5curb7GO++8Y/8S37JlixITEx2Oz58/X7m5uVq3bp02bdqkJUuWaMGCBdqwYYMGDRrU6uvW7xt07rnnKiQkRD//+c/16KOPKjQ01KUyamtrFRMTo+eee06BgYFKTU3Vvn379OSTT9oTigULFmjBggX255w4cUIbN27UtGnT7PuKiooavO6m9OrVyyH2oUOH6ttvv9WTTz5pT3RciUuqW1RVkiorKyXVJVVnnnmmS3G4Y9SoUdq+fbsOHTqk559/Xtdee602bdrUIKlz1datW/Xb3/5W27Zta7aWzxPvsbvojAzAr4SHBLXrw1vOOOMM3XrrrXrhhRe0bds2FRUV6W8rVkiSekT31NGjR1tV7tGjRzV16lTNmjVLDz74oO68804dOnSowXk9e/bUNddco9/85jfauXOnEhIS9Jvf/KbRMuPi4iRJpaWlDvtLS0vtxxozbNgwVVdX66uvvrKX01gZ9a8RHx+vAQMGODQHnX322SopKZHVapVUNyJo+/bt9kdaWpoefvhhh322jrRxcXENRpNVV1fryJEjLca+e/du+7YrcUl1I7wk6bTTTpMk7d27V926dWv2UT9pc1XXrl115pln6vzzz9eLL76ooKAgvfjii26XY/P+++/rwIEDSkxMVFBQkIKCgrRnzx7dc8896tevX5PPc36PvYEaHQDo5Pr166fw8HBVVNaNejpr4LkqKSmx91txx/Tp0xUXF2fv8/Haa68pJydHK75PohoTEhKiM844o8lRV/3791dcXJzy8/M1ZMgQSXUjqDZt2uTQ/8fZ9u3bFRAQYK9lyMjI0IMPPqiqqioFBwdLqqt9Sk5Otr/O4cOH6+WXX1Ztba296e7zzz9XfHy8vW9PdHS0Q3NPWFiYYmJiGq05ycjIUFlZmbZu3WofRbVu3TrV1tZq2LBhzcYeHx9v33YlLknasWOH+vTpY29G8lbTlbPa2lqdOnWq1c+/4YYbHPpgSXX9p2644QZNmTKlyec5v8feQKIDAJ3I3LlzVVlZqbFjxyopKUllZWV65plnVFVVpczMS1Qt6axzzlWvXr304Ycf6vLLL3e57FWrVmnlypXaunWrvXPoSy+9pLS0NP3jH//QVVddpdWrV+uvf/2rJk2apAEDBsgwDL3++ut68803tWzZskbLtVgsuvvuuzVv3jz96Ec/Uv/+/fXQQw8pISFB48ePl1TX6XfTpk0aNWqUunfvroKCAv3iF7/Q//3f/9mTmMmTJ+vXv/61brnlFt1///3asWOHfvvb3+rpp5+2X+uOO+7Qs88+qxkzZmj69On64osvtGDBgmb7jTTn7LPP1ujRozV16lQtWbJEVVVVmjZtmiZNmmSv9XnppZcUEhKiH//4x5Kkf/7zn1q6dKleeOEFt+N6//33demll9q33W26On78uENNUnFxsbZv367o6GglJiaqoqJC8+fP17hx4xQfH69Dhw4pLy9P+/bt0zXXXGN/3sUXX6wJEybYm/NaKrdnz57q2bOnQyzBwcGKi4tTcnKyJNfeY69weYyXSTG8HDCf5oakdhZNDS9ft26dcdVVVxl9+/Y1QkJCjNjYWGP06NHG+++/bx9e/vHXR41f3nuvMWnSpEbLbmx4+cGDB42YmBhj/vz5Dc6fP3++ERMTYxw8eND48ssvjalTpxoDBgwwwsLCjB49ehhDhw5tcni2TW1trfHQQw8ZsbGxRmhoqHHxxRcbu3btsh/funWrMWzYMCMyMtLo0qWLcfbZZxsLFiwwTp486VDOxx9/bIwYMcIIDQ01evfubTz22GMNrvWf//zHGDZsmBEaGmqcfvrpxvz5843q6uomY2tueLlhGMbhw4eN6667zujWrZsRERFhTJkyxfjuu+/sx5cvX26cffbZRnh4uBEREWGkp6cbK1eudDuuEydOGJGRkUZBQUGTsbRk/fr1hqQGj+zsbPs1JkyYYCQkJBghISFGfHy8MW7cOGPz5s0O5SQlJRlz5sxxudzGOA8ld/U9rs8Tw8sthmEY3kujOr7y8nJFRkbq2LFjioiI8Fi5ldZq+zDWooezvNpWD8DRyZMnVVxcrP79+7epU25nU1Nr6L/f1s2q2zPghM4ddI62bdumpKQkH0cGV/z+97/XqlWr9O9//9vXoXQYzd3Lrn5/0xkZAEwoLi5OL774ovbu3evrUOCi4OBgLV682NdhmA7VDABgUrb+L+gcbr31Vl+HYErU6AAAANMi0QFgWn7eBRHo9DxxD/ttopOXl6eUlBQNHTrU16EA8DDbHCu22WUBdE62e9h2T7eG3/bRycnJUU5Ojr3XNgDzCAwMVI8ePewz2oaHh7u0+GRnV1NryKium2X35MmTCgww/2uGORmGocrKSh04cEA9evRo1YKsNn6b6AAwN9v0/M7T95tZrWHoQNlJSVJQZRcF+EFyB3Pr0aNHs0ttuIJEB4ApWSwWxcfHKyYmRlVVVb4Op12csFbrtlUfSJJWTx+hMObvQicWHBzcppocG+4CAKYWGBjokQ/LzqA2oFr7vquRJIV26aIuJDqA/3ZGBgAA5keiAwAATItEBwAAmBaJDgAAMC0SHQAAYFokOgAAwLRIdAAAgGmR6AAAANMi0QEAAKZFogMAAEyLRAcAAJgWiQ4AADAtEh0AAGBaJDoAAMC0SHQAAIBpkegAAADTItEBAACm1ekTnbKyMqWlpWnIkCE655xz9Pzzz/s6JAAA0EEE+TqAturevbs2bNig8PBwVVRU6JxzztGVV16pnj17+jo0AADgY52+RicwMFDh4eGSpFOnTskwDBmG4eOoAABAR+DzRGfDhg264oorlJCQIIvFoldffbXBOXl5eerXr5+6dOmiYcOGafPmzQ7Hy8rKNHjwYPXp00f33nuvevXq1U7RAwCAjszniU5FRYUGDx6svLy8Ro+vWLFCubm5mjNnjrZt26bBgwcrKytLBw4csJ/To0cPffzxxyouLtbLL7+s0tLSJq936tQplZeXOzwAAIA5+TzRGTNmjObNm6cJEyY0evypp57S1KlTNWXKFKWkpGjJkiUKDw/X0qVLG5wbGxurwYMH6/3332/yeo8++qgiIyPtj759+3rstQAAgI7F54lOc6xWq7Zu3arMzEz7voCAAGVmZqqgoECSVFpaqu+++06SdOzYMW3YsEHJyclNljlz5kwdO3bM/vj666+9+yIAAIDPdOhRV4cOHVJNTY1iY2Md9sfGxuqzzz6TJO3Zs0e33XabvRPy9OnTNWjQoCbLDA0NVWhoqFfjBgAAHUOHTnRckZ6eru3bt/s6DAAA0AF16KarXr16KTAwsEHn4tLSUsXFxfkoKgAA0Fl06EQnJCREqampys/Pt++rra1Vfn6+MjIy2lR2Xl6eUlJSNHTo0LaGCQAAOiifN10dP35cu3fvtm8XFxdr+/btio6OVmJionJzc5Wdna20tDSlp6dr0aJFqqio0JQpU9p03ZycHOXk5Ki8vFyRkZFtfRkAAKAD8nmiU1hYqFGjRtm3c3NzJUnZ2dlavny5Jk6cqIMHD2r27NkqKSnRkCFD9PbbbzfooAwAAODM54nOyJEjW1yyYdq0aZo2bVo7RQQAAMyiQ/fR8Sb66AAAYH5+m+jk5OSoqKhIW7Zs8XUoAADAS/w20QEAAOZHogMAAEyLRAcAAJgWiQ4AADAtv010GHUFAID5+W2iw6grAADMz28THQAAYH4kOgAAwLRIdAAAgGmR6AAAANPy20SHUVcAAJif3yY6jLoCAMD8/DbRAQAA5keiAwAATItEBwAAmBaJDgAAMC0SHQAAYFp+m+gwvBwAAPPz20SH4eUAAJif3yY6AADA/Eh0AACAaZHoAAAA0yLRAQAApkWiAwAATItEBwAAmBaJDgAAMC2/TXSYMBAAAPPz20SHCQMBADA/v010AACA+ZHoAAAA0yLRAQAApkWiAwAATItEBwAAmBaJDgAAMC0SHQAAYFokOgAAwLRIdAAAgGn5baLDEhAAAJif3yY67bkERKW1RoZheP06AADAkd8mOu0pbd5aXbOkgGQHAIB2RqLjJWHBgUpLirJvF+45qhNVNT6MCAAA/0Oi4yUWi0Urb89Q4axMX4cCAIDfItHxIovFovCQQF+HAQCA3yLRAQAApkWiAwAATCvI1wH4k0prXWfksOBAWSwWH0cDAID5kei0o7R5a+v+TYrSytszSHYAAPAymq68zHmYucRQcwAA2gs1Ol5mG2Z+oqpGldYae60OAADwPrcTnZqaGi1fvlz5+fk6cOCAamtrHY6vW7fOY8GZRd0wc3JKAADam9vfvjNmzNDy5ct12WWX6ZxzzqGfSStVWmvolAwAgJe5nej89a9/1d/+9jeNHTvWG/H4jbR5a+mUDACAl7ndGTkkJERnnnmmN2IxvcbWvzpcYVWltVqV1moW/QQAwMPcTnTuuece/fa3v+30X8p5eXlKSUnR0KFD2+2aja1/lTZvrVJmr1HK7DWscA4AgIe53XT1wQcfaP369Xrrrbc0cOBABQcHOxz/5z//6bHgvCknJ0c5OTkqLy9XZGRku13XYrGoZ9cQpSVFqXDPUYdjtmHndFwGAMAz3P5G7dGjhyZMmOCNWPxG/SHnkhh2DgCAl7id6CxbtswbcfgdhpwDAOB9zIzcwVRaa+z9dAzDoKMyAABt4FKVwnnnnaf8/HxFRUXpxz/+cbPDobdt2+ax4PyRbdj5336eoWv+UKCt3/fjse07Wc3CoAAAuMqlROdnP/uZQkNDJUnjx4/3Zjx+yTbs3NY5uXDPUe0rO2FPcmz7Llv8gXbuL5fEwqAAALjCYvh5m4ht1NWxY8cUERHhszgMw9DhCqtbnZKLHs6inw8Au0prtVJmr5HE5wPMz9Xv71bfBYWFhdq5c6ckKSUlRampqa0tCmp62HlKfISKvq/FAQAA7nE70fnmm2903XXX6cMPP1SPHj0kSWVlZfrJT36iv/71r+rTp4+nY/QbzsPOJckwpIFz1jR6PutlAQDQPLdHXd16662qqqrSzp07deTIER05ckQ7d+5UbW2tbr31Vm/E6Fdsw85tj+ZymLR5a5lNGQCAZrhdo/Pee+/pP//5j5KTk+37kpOTtXjxYl1wwQUeDQ6Nq9+cxWzKAAA0ze0anb59+6qqqqrB/pqaGiUkJHgkKPzAeSHQtKQorZ4+wmG9LMOQfb4d5twBAOAHblcDPPnkk5o+fbry8vKUlpYmqa5j8owZM/Sb3/zG4wH6O+d+O7Y+OeEhgfZzrl5SYB92LtXV+Ky8PUPhIfTfAQD4N5cSnaioKIcvzIqKCg0bNkxBQXVPr66uVlBQkG6++Wbm2fGClpaL2Ok0Kqtof7kGzlnDXDsAAL/nUqKzaNEiL4cBT9jyYKayl26m/w4AAN9z6RswOzvb23HAA7qGBuqNu0a4PfEgAABmxaKenVRjnZQb67/jrP5CoXRcBgCYHW0anVRTnZSd2UZkSVKXoECHhUIl1swCAJgbiU4n1lInZUm6fPEHKj5UIanx5SToxwMAMDOarkzOluRIckhy3r9vlP3n+vPw0JQFADATtxKdqqoqBQUFaceOHd6KBx7g3H+nf6+uDc9xmocnZfYapcxew5ISAABTcau9Ijg4WImJiaqpqWn5ZPiMc/+dLkGBuvYPBfZV0W0dl23qz8NDUxYAwEzc/jZ78MEH9cADD+hPf/qToqOjvRGTW77++mvdcMMNOnDggIKCgvTQQw/pmmuu8XVYPufcf8e543L9FdIBADArtxOdZ599Vrt371ZCQoKSkpLUtatjs8i2bds8FpwrgoKCtGjRIg0ZMkQlJSVKTU3V2LFjG8Tl71zpuGxTaa1pchQXAACdiduJTkdb4iE+Pl7x8fGSpLi4OPXq1UtHjhwh0WmDtHlrGXYOADAFtxOdOXPmeDSADRs26Mknn9TWrVu1f/9+rVq1qkEylZeXpyeffFIlJSUaPHiwFi9erPT09AZlbd26VTU1Nerbt69HY/QX9Yef01cHAGAGrRpeXlZWphdeeEEzZ87UkSNHJNU1We3bt8/tsioqKjR48GDl5eU1enzFihXKzc3VnDlztG3bNg0ePFhZWVk6cOCAw3lHjhzRjTfeqOeee67Z6506dUrl5eUOD3/U2MzKq6ePUOGsTB9GBQCAZ7n95/onn3yizMxMRUZG6quvvtLUqVMVHR2tf/7zn9q7d6/++Mc/ulXemDFjNGbMmCaPP/XUU5o6daqmTJkiSVqyZIneeOMNLV26VL/61a8k1SUv48eP169+9Sv95Cc/afZ6jz76qH7961+7FaMZNTWzcv3lI+irAwDo7Nyu0cnNzdVNN92kL774Ql26dLHvHzt2rDZs2ODR4KxWq7Zu3arMzB9qGQICApSZmamCggJJdWs33XTTTbrooot0ww03tFjmzJkzdezYMfvj66+/9mjMnYmtg3J4SFCjyUzavLXMqwMA6NTcTnS2bNmin//85w329+7dWyUlJR4JyubQoUOqqalRbGysw/7Y2Fj7tT788EOtWLFCr776qoYMGaIhQ4bo008/bbLM0NBQRUREODzwA+cmLVtfHQAAOiO3m65CQ0Mb7dfy+eef67TTTvNIUO4YMWKEamtr2/26ZmVr0jpcYVXavLVNnmcYRosLigIA4Gtu1+iMGzdODz/8sKqqqiTVfTHu3btX999/v6666iqPBterVy8FBgaqtLTUYX9paani4uLaVHZeXp5SUlI0dOjQNpVjRo311am/DpZhGCwbAQDoFNxOdBYuXKjjx48rJiZGJ06c0IUXXqgzzzxT3bt31/z58z0aXEhIiFJTU5Wfn2/fV1tbq/z8fGVkZLSp7JycHBUVFWnLli1tDdP00uatdUhoTlTVaOv3y0lING8BADout5uuIiMj9c477+iDDz7QJ598ouPHj+u8885z6DDsjuPHj2v37t327eLiYm3fvl3R0dFKTExUbm6usrOzlZaWpvT0dC1atEgVFRX2UVjwDltfnUISGgBAJ9bq2eBGjBihESNGtDmAwsJCjRo1yr6dm5srScrOztby5cs1ceJEHTx4ULNnz1ZJSYmGDBmit99+u0EHZXhW/eHnldaaZvvrAADQUbUq0cnPz9fTTz+tnTt3SpLOPvts3X333a2q1Rk5cmSL/TumTZumadOmtSZUtIE762MBANARud1H53e/+51Gjx6t7t27a8aMGZoxY4YiIiI0duzYJmc37ojojAwAgPm5/ef6ggUL9PTTTzvUsNx1110aPny4FixYoJycHI8G6C05OTnKyclReXm5IiMjfR0OAADwArdrdMrKyjR69OgG+y+99FIdO3bMI0Gh4zKMuuHmAAB0Bq2aR2fVqlUN9r/22mu6/PLLPRIUOq6rlxQ02jHZNteO7cG8OgCAjsDtpquUlBTNnz9f7777rn0um40bN+rDDz/UPffco2eeecZ+7l133eW5SNEh7Nz/w6zYKfERKvp+2zn5SYmP0MrbM2SbMJnZkwEAvmAx3PzTu3///q4VbLHof//7X6uCag95eXnKy8tTTU2NPv/8cx07dox1r5pQaa1Wyuw1DvsKZ2UqOjxE1/6hwGGunaakJUV9n/iQ7ADeUv9eLXo4i1GTMDVbH9uWvr/dvguKi4vbFFhHQWdk1zlPHpiWFKWeXUMc5tqxMQzpmiUF9poeG9tkg3zwAgDaE986aJFzQlO/GaqxuXbeuGuE/VwmGwQA+BKJDlzizuSBTDQIAOgo3B51BQAA0FnwZzc6FNvq6BIjtQAAbee3iU79UVfoGAzD0NVLCrS1XqdnRmoBANrC7aart99+Wx988IF9Oy8vT0OGDNHkyZN19GjLw4w7ipycHBUVFWnLli2+DgXfO1FVY09ypB9GagEA0FpuJzr33nuvysvrhg5/+umnuueeezR27FgVFxcrNzfX4wECAAC0Vqvm0UlJSZEk/eMf/9Dll1+uBQsWaNu2bRo7dqzHAwQAAGgtt2t0QkJCVFlZKUlau3atLr30UklSdHS0vaYHAACgI3C7RmfEiBHKzc3V8OHDtXnzZq1YsUKS9Pnnn6tPnz4eDxDmVX+ElcSq6AAAz3M70Xn22Wd155136u9//7t+//vfq3fv3pKkt956S6NHj/Z4gDCP+olMl6BAXfOHAofOxwAAeJrbiU5iYqJWr17dYP/TTz/tkYDaC8PL21/9pSDqr3zurLljAAC4w6VEx52+N51lBXAW9WwfzguC2tRPZApnZSo8JNC+bRjSwDmOq6UDANAaLiU6PXr0aHHSNsMwZLFYqCGBA+cFQRtb5DM8JNBhbaxKa3W7xggAMC+XEp3169d7Ow6YWHst8snyEQAAZy59+1x44YXejgN+wrkpKy0pSmHBgU2ebxg/1PA0l7ywfAQAoDGt/jO7srJSe/fuldVqddh/7rnntjkomJdzU1ZLNS9XLynQzu/78zSXvDS1fER71CQBADout78FDh48qClTpuitt95q9Dh9dNASd5qydtbvtLznqA5XWBUeEmivBbIlTIbh+TgBAJ2f24nO3XffrbKyMm3atEkjR47UqlWrVFpaqnnz5mnhwoXeiBGws3VkTk2KkiR7LU7/Xl19FhMAoONyO9FZt26dXnvtNaWlpSkgIEBJSUm65JJLFBERoUcffVSXXXaZN+KEn3OeW8d5osHiQxXtHRIAoBNwe62riooKxcTESJKioqJ08OBBSdKgQYO0bds2z0bnRXl5eUpJSdHQoUN9HQqc2Dos26QlRWn19BEqejhLhbMyG5xfvzYnJb5zzOMEAGgfbtfoJCcna9euXerXr58GDx6sP/zhD+rXr5+WLFmi+Ph4b8ToFUwY2HE11WG5qX49q6ePkK1/MpMNAgDqczvRmTFjhvbv3y9JmjNnjkaPHq2//OUvCgkJ0fLlyz0dH/yUOx2WLRbZz60/2aBtba3GRnYx5w4A+Ae3E53/+7//s/+cmpqqPXv26LPPPlNiYqJ69erl0eCAtrB1XHYels6cOwDgP9zuo+MsPDxc5513HkkOOgTn/j3SD3Pq2DQ15w4AwHxcqtHJzc3VI488oq5duyo3N7fZc5966imPBAa0Rv3+PY2tqwUA8C8uJTofffSRqqqqJEnbtm1rsoqfqn90BO21thYAoONze1HPd99911uxAC1yd60sAIB/c+vP3qqqKoWFhWn79u0655xzvBUT0CR318qyqbTWMLoKAPyQW52Rg4ODlZiYyHpW8Clb01R4SJDLiUvavLW6ZkmBDBbFAgC/4vaoqwcffFAPPPCAjhw54o14AI9xHoHl7ugqwzBUaa1WpbWaBAkAOim3e2w+++yz2r17txISEpSUlKSuXR0XU+xMy0DA3GzNXIcrrG6PvmKuHQAwB7cTnZ/97Gem+LDPy8tTXl4ezXAmV9fM5X5n5abm2mE0FwB0Lm5/as+dO9cLYbQ/1roCAMD83O6jc/rpp+vw4cMN9peVlen000/3SFAAAACe4Hai89VXXzXa3HPq1Cl98803HgkK8KW6Tsg0aQKAGbjcdPWvf/3L/vOaNWscmntqamqUn5+v/v37ezY6oJ05d0IGAHRuLic648ePl1TXuTM7O9vhWHBwsPr166eFCxd6NDigvdgmFHTuhJwSH6Gi/eU+jAwA0BYuJzq1tbWSpP79+2vLli2sVg5TSZu3VmlJUXrp5nT7vsJZmQoLDtTAOWt8GBkAoC3c7qNTXFxMkgPTSImPsP/sPKFgeEigTDCTAgD4NbcTHUmaNm0aMyOjU6q01qj+JMcrb89Q4axM3wUEAPAqlxOd+iOqXn75ZR0/flySNGjQIH399deejwzwAtuaVzYWi1o1oSAAoHNwOdE566yzlJSUpMmTJ+vkyZP25Oarr75SVVWV1wIE2sp5zavmOhefYFg5AJiKy4lOWVmZVq5cqdTUVNXW1mrs2LEaMGCATp06pTVr1qi0tNSbcQKtZlvzypUmqgueWN8OEQEA2ovLiU5VVZXS09N1zz33KCwsTB999JGWLVumwMBALV26VP3791dycrI3YwVarbk1r5xrfKS6RTzDgmnSAoDOzuXh5T169NCQIUM0fPhwWa1WnThxQsOHD1dQUJBWrFih3r17a8uWLd6MFfAKW41P/RFXYcGBpli8FgD8ncs1Ovv27dOsWbMUGhqq6upqpaam6oILLpDVatW2bdtksVg0YsQIb8YKtIlzzU39Wpu6Gp8g+6OxJKduxJbRYD8AoOOyGK345I6KitKGDRu0c+dO3XjjjYqLi1NpaanS09P13nvveSNOr7GtXn7s2DFFRES0/AR0aoZh2GtuXKm1qbRWK2X2DxMGpiVFaeXtGdT2oEOq//+16OEshYe4XGkPdDqufn+3ah4dSYqMjNS1116r4OBgrVu3TsXFxbrzzjtbW1y7y8vLU0pKioYOHerrUNCO6tfcuJKsONcCOU8qCADo2FqV6HzyySfq06ePJCkpKUnBwcGKi4vTxIkTPRqcN+Xk5KioqIh+RWiWOyO2AAAdT6vqNfv27Wv/eceOHR4LBuiInEdsVTrNtdNYE5i7TWQAAO+gARdwU9q8tY7bTv12DMPQ1UsK7Kug068HAHyn1X10AH/S2Fw7NoV7jupwhdU+IutEVY09ybEdP1FVN2Kr0lqtSms1o7cAoJ1QowO4oLG5diqtNfbanbR5a+01N40xDFHLAwA+QI0O4CLnuXZ6dg1xeURWU7U8AADvItEBWokRWQDQ8ZHoAG3Q3BpaAADfI9EBPKhumYiG+09YaaYCAF+gMzLgQWnz1iolvuFU5Bc8sd4H0QAAqNEB2sh56HnR/nL7z85JT2NJEADAe6jRAdrI1in5cIW1wWSCdUPIf9g2DGngnDUCALQPEh3AA5rqlGyxyGEF6UprdXuGBQB+j6YrwEOcm7DSkqIUFsyILADwJWp0AA9xnj25My3mWX8RUpvOFD8ANIVEB/Ag2+zJnYnzIqQ2LFMBwAxougL8nPPyFDYsUwHADEh0ANgVzspkSQsAptK56tgBeBXLWQAwG2p0AACAaZHoAAAA0zJFojNhwgRFRUXp6quv9nUogMsqrTWqtFar0lotwzBkGIbDdnPqn1v/0Vg57pQLAGZjij46M2bM0M0336yXXnrJ16EALqu/XETq9xMN2kY/NTe0u6nh4I2V4065AGBGpkh0Ro4cqXfffdfXYQAtss2eXOiUpDgnLbah3Y3NydPUcPDGynGnXAAwI583XW3YsEFXXHGFEhISZLFY9OqrrzY4Jy8vT/369VOXLl00bNgwbd68uf0DBTzANnty0cNZKno4q8Wh3I03Q/0wt03hrEyXygEAf+XzP+sqKio0ePBg3XzzzbryyisbHF+xYoVyc3O1ZMkSDRs2TIsWLVJWVpZ27dqlmJgYt6936tQpnTp1yr5dXl7epvgBd7k6e7JhyKGJyrkZSqobDk7tDAA0zec1OmPGjNG8efM0YcKERo8/9dRTmjp1qqZMmaKUlBQtWbJE4eHhWrp0aauu9+ijjyoyMtL+6Nu3b1vCB9rEeSHQlPgI+8/OTVRb9xx12GbRUABoWYf+U9BqtWrr1q2aOXOmfV9AQIAyMzNVUFDQqjJnzpyp3Nxc+3Z5eTnJDnzGeSFQw5AGzlnT4vMKZ2WqZ9cQOhUDQAs6dKJz6NAh1dTUKDY21mF/bGysPvvsM/t2ZmamPv74Y1VUVKhPnz5auXKlMjIyGi0zNDRUoaGhXo0bcEf9pqxKa7VLzwkPYWVxAHBFh050XLV27dqWTwLgNlvH57BgEisAnVOHTnR69eqlwMBAlZaWOuwvLS1VXFycj6IC/Idtrh/m3wHQWfm8M3JzQkJClJqaqvz8fPu+2tpa5efnN9k05aq8vDylpKRo6NChbQ0T6JDqd2yu/3NLnDtISz/MvwMAnY3Pa3SOHz+u3bt327eLi4u1fft2RUdHKzExUbm5ucrOzlZaWprS09O1aNEiVVRUaMqUKW26bk5OjnJyclReXq7IyMi2vgygw6mrgan72dVOzpJjB+lKa43DDM4A0Nn4PNEpLCzUqFGj7Nu2EVHZ2dlavny5Jk6cqIMHD2r27NkqKSnRkCFD9PbbbzfooAyYzQlr22pQLBY12cnZMAyHkV4Nn+vaXD8A0NH5/JNs5MiRLS40OG3aNE2bNq2dIgI6hgueWO+Vcp0nInSnWQsAOpsO3UfHm+ijg46osf4xriQizs9rbjJB54kIi/YzOzgA8/J5jY6v0EcHHZHzBIKSa/1rnJ/HcHAAqOO3iQ7QUTn3j3F1EkH61QBAQ37bdAUAAMyPP/8AtEr9kVveaCrzdvkA/AOJDgC3GYbhMHLL0zMne7t8AP7Db5uuGHWFzsKdEVXuqj+iy51ynUdueXrmZG+XD8B/+G2NDqOu0Fl4c0RV/dmTaR4CYEZ+m+gAnYm3RlTVnz0ZAMzIb5uuAACA+fGnHOBn2rqGlmEYqnSxDFdGTjG6CoA3kegAfqYta2g5j4Zy59zGRk41dQ4AeIrfNl0x6gr+pLE1tFozest5NFRz63C5MnKK0VUAvM1va3QYdQV/0tgaWm1tJiqclamw4MAW1+ECAF/y20QH8DeeHrkVHuKZuXwAwJv8tukKAACYH4kOAAAwLZquALjEMORyR2F3hqC35ToA0BISHQAuuXpJgXbuL2/xPHeGoLflOgDgCr9tumJ4OeCe+slHc0PT3RmC3tJ13H0uADjz2xodhpcDrVM4K1M9u4a4NDS9LUPQGb4OwBP8NtEB0DrhIa7Pv9OWIegMXwfgCX7bdAUAAMyPRAcAAJgWiQ4AADAtEh0AAGBadEYG4BGV1poWV0M3DKnSWu3wnLYwDMM+uWBbFykFYE4kOgA8Im3eWqUlRemlm9ObPMeTkwE6T0yYlhSllbdnkOwAcOC3TVdMGAi0XVhwoNKSouzbhXuONrt8Q1NJTmpSlFLrldPchIQ2zhMTtnRtAP7Jb2t0mDAQaDuLxaKVt2focIVVafPWuvy8wlmZDvPk2JIamqEAeJrfJjoAPMNisbg9uV94SKDCQxp+/DS2DwDawm+brgAAgPmR6AAAANMi0QEAAKZFogMAAEyLRAcAAJgWiQ4AADAtxnIC6DTqL/lgGG0vw3n+Hts+5vABzINEB0CnYBhyWPIhJT6iFWU4Lhthm425/gzLLCUBmIvfNl2xBATQMuclHlxZmsFbnJd8KGrFmlnOZWzdc9RhW2IpCcBs/LZGhyUggJbZlnjwh6UZ3r9vlC54Yr2vwwDgYX6b6ABwTd0SD+b/qAhzcxkLAJ2D3zZdAQAA8yPRAQAApkWiAwAATItEBwAAmBaJDgAAMC0SHQAAYFokOgAAwLRIdAAAgGmR6AAAANMi0QEAAKZFogMAAEzL/AvYAGhXJ6zeWfnblXINQ6q0Vkv6YQFSwzDsi5IahmvXqrTWNLqAaf2ybNeQ1OZFT+uXa+aFUwFfINEB4FHeWgHclXKvXlKgnfvLJUlpSVH6288zdM0fCrR1z1FJUkp8hEvXSpu3VmlJUVp5e4Y96TAMQ1cv+aEsSUpNipIk+z7n57jCudzWlAGgaX7bdJWXl6eUlBQNHTrU16EAnV5YcKDSvv/St0lLirLXeHiy3NSkKHuCITkmL7YkR5IK9xzVkUqrQ2JSVO94S9cq3HPUofbmRFWNQ1lSXYJTf5/zc1zhXG5rygDQNL+t0cnJyVFOTo7Ky8sVGRnp63CATs1isWjl7RkNmnXaWivRVLmSHJqjBs5Z06br1F1LWnl7hg5XWJU2b22z575/3yiv1VwB8Cy/TXQAeJbFYlF4iOc/Upoq17bP1ifHc9dquRYqzIVzAHQMftt0BQAAzI9EBwAAmBaJDgAAMC0SHQAAYFokOgAAwLRIdAAAgGmR6AAAANMi0QEAAKZFogMAAEyLRAcAAJgWiQ4AADAtEh0AAGBaJDoAAMC0SHQAAIBpkegAAADTItEBAACmRaIDAABMi0QHAACYFokOAAAwLRIdAABgWqZIdFavXq3k5GT96Ec/0gsvvODrcAAAQAcR5OsA2qq6ulq5ublav369IiMjlZqaqgkTJqhnz56+Dg0AAPhYp6/R2bx5swYOHKjevXurW7duGjNmjP7973/7OiwAANAB+DzR2bBhg6644golJCTIYrHo1VdfbXBOXl6e+vXrpy5dumjYsGHavHmz/di3336r3r1727d79+6tffv2tUfoADq4E9aaNj2/0lojwzDa/BzDMFRprValtdrhWN1+12KsX4ZzWfWPuco5pua2mzqnteU39hrayp3Y2lq+J67h7Xg7go7yGn3edFVRUaHBgwfr5ptv1pVXXtng+IoVK5Sbm6slS5Zo2LBhWrRokbKysrRr1y7FxMS4fb1Tp07p1KlT9u3y8vI2xQ+g47rgifVten7avLVKS4rSytszWvUci8UiwzB09ZICbd1ztO54vfLq72+Ocxn2ayVF6W8/z9A1f3CtnKbKS02KkqQmtxvbV/81ulu+82toqpzWvh5PlNlc+W29hrfj7QicX2PRw1kKD/FNyuHzGp0xY8Zo3rx5mjBhQqPHn3rqKU2dOlVTpkxRSkqKlixZovDwcC1dulSSlJCQ4FCDs2/fPiUkJDR5vUcffVSRkZH2R9++fT37ggA4CAsOVNr3X3RS3Yd6WHCgV6+ZEh/R7LbzvvoxOcdbuOeoTlQ51ro4n5OaFGX/Mnd+zomqGocvSNsx5/2NxWjjfG79so5UWhsca+lLxbm8rXuONrvd2L7Gfi+ulu/8Gpoqx1VN/Y49pbHff1uu4e14O4Km/s/6gs9rdJpjtVq1detWzZw5074vICBAmZmZKigokCSlp6drx44d2rdvnyIjI/XWW2/poYcearLMmTNnKjc3175dXl5OsgN4kcVi0crbM+wf5GHBgV7/y7Xur+Mftg1DGjhnTZPn1I/JFu/hCqvS5q1ttHyLRQ1ek6Rmn9OcwlmZCgsObBBjU+dKatV1WuP9+0a1uWasvsJZmQoPCVSltabdXoMnefr34Q9s/799pUMnOocOHVJNTY1iY2Md9sfGxuqzzz6TJAUFBWnhwoUaNWqUamtrdd999zU74io0NFShoaFejRuAI4vF0q7V1haLHK7XWN8V53Mcj1kUHtL8B3Njr6ml5zTFnee19hqtFebh64WHBPqsCcMTPP378AfhId7/46Y5nfd/Wz3jxo3TuHHjfB0GAADoYHzeR6c5vXr1UmBgoEpLSx32l5aWKi4urk1l5+XlKSUlRUOHDm1TOQAAoOPq0IlOSEiIUlNTlZ+fb99XW1ur/Px8ZWS4PgqiMTk5OSoqKtKWLVvaGiYAAOigfN50dfz4ce3evdu+XVxcrO3btys6OlqJiYnKzc1Vdna20tLSlJ6erkWLFqmiokJTpkzxYdQAAKAz8HmiU1hYqFGjRtm3bSOisrOztXz5ck2cOFEHDx7U7NmzVVJSoiFDhujtt99u0EEZAADAmc8TnZEjR7Y4Y+K0adM0bdq0dooIAACYRYfuo+NNdEYGAMD8/DbRoTMyAADm57eJDgAAMD8SHQAAYFokOgAAwLT8NtGhMzIAAObnt4kOnZEBADA/n8+j42u2OXzKy8t9HAmA1qi0Vqv2VKXDvvLyclU7rV7e0jnNlWv7fKi/3dhznZ9THRLU5LWdy2uu/OZi+a68vE2vrSWNld/SdZorv7HfS0vxtqQ1729ry6//+2jtNbwdb0fgyfe3KbZ7oaW5+CxGS2eY3DfffKO+ffv6OgwAANAKX3/9tfr06dPkcb9PdGpra/Xtt9+qe/fuslgsHimzvLxcffv21ddff62IiAiPlAnv4j3rnHjfOh/es86pI75vhmHou+++U0JCggICmu6JY556slYKCAhoNhNsi4iIiA7zHwKu4T3rnHjfOh/es86po71vkZGRLZ7jt52RAQCA+ZHoAAAA0yLR8YLQ0FDNmTNHoaGhvg4FLuI965x43zof3rPOqTO/b37fGRkAAJgXNToAAMC0SHQAAIBpkegAAADTItEBAACmRaLjBXl5eerXr5+6dOmiYcOGafPmzb4OCU2YO3euLBaLw+Oss87ydVhwsmHDBl1xxRVKSEiQxWLRq6++6nDcMAzNnj1b8fHxCgsLU2Zmpr744gvfBAtJLb9nN910U4N7b/To0b4JFpKkRx99VEOHDlX37t0VExOj8ePHa9euXQ7nnDx5Ujk5OerZs6e6deumq666SqWlpT6K2DUkOh62YsUK5ebmas6cOdq2bZsGDx6srKwsHThwwNehoQkDBw7U/v377Y8PPvjA1yHBSUVFhQYPHqy8vLxGjz/xxBN65plntGTJEm3atEldu3ZVVlaWTp482c6Rwqal90ySRo8e7XDvvfLKK+0YIZy99957ysnJ0caNG/XOO++oqqpKl156qSoqKuzn/OIXv9Drr7+ulStX6r333tO3336rK6+80odRu8CAR6Wnpxs5OTn27ZqaGiMhIcF49NFHfRgVmjJnzhxj8ODBvg4DbpBkrFq1yr5dW1trxMXFGU8++aR9X1lZmREaGmq88sorPogQzpzfM8MwjOzsbONnP/uZT+KBaw4cOGBIMt577z3DMOruq+DgYGPlypX2c3bu3GlIMgoKCnwVZouo0fEgq9WqrVu3KjMz074vICBAmZmZKigo8GFkaM4XX3yhhIQEnX766br++uu1d+9eX4cENxQXF6ukpMThvouMjNSwYcO47zq4d999VzExMUpOTtYdd9yhw4cP+zok1HPs2DFJUnR0tCRp69atqqqqcrjXzjrrLCUmJnboe41Ex4MOHTqkmpoaxcbGOuyPjY1VSUmJj6JCc4YNG6bly5fr7bff1u9//3sVFxfrggsu0Hfffefr0OAi273Ffde5jB49Wn/84x+Vn5+vxx9/XO+9957GjBmjmpoaX4cGSbW1tbr77rs1fPhwnXPOOZLq7rWQkBD16NHD4dyOfq/5/erl8G9jxoyx/3zuuedq2LBhSkpK0t/+9jfdcsstPowMMLdJkybZfx40aJDOPfdcnXHGGXr33Xd18cUX+zAySFJOTo527Nhhij6L1Oh4UK9evRQYGNigB3ppaani4uJ8FBXc0aNHDw0YMEC7d+/2dShwke3e4r7r3E4//XT16tWLe68DmDZtmlavXq3169erT58+9v1xcXGyWq0qKytzOL+j32skOh4UEhKi1NRU5efn2/fV1tYqPz9fGRkZPowMrjp+/Li+/PJLxcfH+zoUuKh///6Ki4tzuO/Ky8u1adMm7rtO5JtvvtHhw4e593zIMAxNmzZNq1at0rp169S/f3+H46mpqQoODna413bt2qW9e/d26HuNpisPy83NVXZ2ttLS0pSenq5FixapoqJCU6ZM8XVoaMQvf/lLXXHFFUpKStK3336rOXPmKDAwUNddd52vQ0M9x48fd/hLv7i4WNu3b1d0dLQSExN19913a968efrRj36k/v3766GHHlJCQoLGjx/vu6D9XHPvWXR0tH7961/rqquuUlxcnL788kvdd999OvPMM5WVleXDqP1bTk6OXn75Zb322mvq3r27vd9NZGSkwsLCFBkZqVtuuUW5ubmKjo5WRESEpk+froyMDJ1//vk+jr4Zvh72ZUaLFy82EhMTjZCQECM9Pd3YuHGjr0NCEyZOnGjEx8cbISEhRu/evY2JEycau3fv9nVYcLJ+/XpDUoNHdna2YRh1Q8wfeughIzY21ggNDTUuvvhiY9euXb4N2s81955VVlYal156qXHaaacZwcHBRlJSkjF16lSjpKTE12H7tcbeL0nGsmXL7OecOHHCuPPOO42oqCgjPDzcmDBhgrF//37fBe0Ci2EYRvunVwAAAN5HHx0AAGBaJDoAAMC0SHQAAIBpkegAAADTItEBAACmRaIDAABMi0QHAACYFokOAAAwLRIdAABgWiQ6ADoti8XS7GPu3LmSpFWrVun8889XZGSkunfvroEDB+ruu+/2aewA2geLegLotPbv32//ecWKFZo9e7Z27dpl39etWzfl5+dr4sSJmj9/vsaNGyeLxaKioiK98847vggZQDsj0QHQacXFxdl/joyMlMVicdgnSa+//rqGDx+ue++9175vwIABrGwO+AmargCYWlxcnP773/9qx44dvg4FgA+Q6AAwtenTp2vo0KEaNGiQ+vXrp0mTJmnp0qU6deqUr0MD0A5IdACYWteuXfXGG29o9+7dmjVrlrp166Z77rlH6enpqqys9HV4ALyMRAeAXzjjjDN066236oUXXtC2bdtUVFSkFStW+DosAF5GogPA7/Tr10/h4eGqqKjwdSgAvIxRVwBMbe7cuaqsrNTYsWOVlJSksrIyPfPMM6qqqtIll1zi6/AAeBk1OgBM7cILL9T//vc/3XjjjTrrrLM0ZswYlZSU6N///reSk5N9HR4AL7MYhmH4OggAAABvoEYHAACYFokOAAAwLRIdAABgWiQ6AADAtEh0AACAaZHoAAAA0yLRAQAApkWiAwAATItEBwAAmBaJDgAAMC0SHQAAYFr/D2Q++YNysY+8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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"
   ]
  },
  {
   "attachments": {},
   "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": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.utils.analysis import extend_trial_data_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function extend_trial_data_file in module skyllh.core.utils.analysis:\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": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 40001/40001 [1:33:15<00:00,  7.15it/s]\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": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TimeLord: Executed tasks:\n",
      "[Generating background events for data set 0.]   0.003 sec/iter (40000)\n",
      "[Generating background events for data set 1.]   0.005 sec/iter (40000)\n",
      "[Generating background events for data set 2.]   0.004 sec/iter (40000)\n",
      "[Generating background events for data set 3.]   0.008 sec/iter (40000)\n",
      "[Generating background events for data set 4.]   0.029 sec/iter (40000)\n",
      "[Generating pseudo data.                     ]   0.045 sec/iter (40000)\n",
      "[Initializing trial.                         ]   0.126 sec/iter (40000)\n",
      "[Get sig probability densities and grads.    ] 2.6e-04 sec/iter (7959160)\n",
      "[Evaluating bkg log-spline.                  ] 3.3e-04 sec/iter (7959160)\n",
      "[Get bkg probability densities and grads.    ] 4.0e-04 sec/iter (7959160)\n",
      "[Calculate PDF ratios.                       ] 1.3e-04 sec/iter (7959160)\n",
      "[Calc pdfratio value Ri                      ]   0.002 sec/iter (3979580)\n",
      "[Calc logLamds and grads                     ] 4.4e-04 sec/iter (3979580)\n",
      "[Evaluate llh-ratio function.                ]   0.008 sec/iter (795916)\n",
      "[Minimize -llhratio function.                ]   0.166 sec/iter (40000)\n",
      "[Maximizing LLH ratio function.              ]   0.166 sec/iter (40000)\n",
      "[Calculating test statistic.                 ] 6.1e-05 sec/iter (40000)\n"
     ]
    }
   ],
   "source": [
    "print(tl)"
   ]
  },
  {
   "attachments": {},
   "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": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-log10(p_local) = 2.89\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": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXiU9bn/8fdtEMEFlE0jIQJCLUEwYDRup6IVhSoC/SGbtlItKEK3X1ulnLaovVqptupPrSJU3FpRjz2IKLX2KBQXCorNKZssBYqBgCxlsYps9++PyfM4GWYmE5LJZDKf13XlSuZ55nnmHhxz57vdX3N3REREAI7KdAAiItJwKCmIiEhISUFEREJKCiIiElJSEBGRUJNMB1Abbdq08Y4dO2Y6DBGRrLJ48eJt7t423rmsTgodO3bkvffey3QYIiJZxcz+meicuo9ERCSkpCAiIiElBRERCWXlmIKZDQAGdOnSJdOhSI7Yv38/5eXl7N27N9OhiKSsWbNmFBQUcPTRR6d8jWVz7aOSkhLXQLPUh3Xr1nHCCSfQunVrzCzT4YhUy93Zvn07e/bsoVOnTlXOmdlidy+Jd526j0RSsHfvXiUEySpmRuvWrWvculVSEEmREoJkmyP5zGblmEJt3TF7Gcs37WZgcXtGlhZmOhwRkQYjZ1sKC9ftYFbZxkyHIXKYlStXsnLlysOOr1+/njPPPLNW9543bx5XXXVVre6RLqNGjeKFF17IdBg5LydbCpMGdGf5pt2ZDkMkq7g77s5RR+Xs35I5Qf91RbLIgQMHuP766+nZsydDhgzhk08+4c477+Scc87hzDPPZMyYMQQzCtesWcNll13GWWedRe/evfnHP/5R5V7vvvsuvXr1Yu3atWzdupW+ffvSu3dvbrrpJk477TS2bdvG+vXr6datG7fccgu9e/fmww8/ZMaMGfTo0YMzzzyT2267Lbzf8ccfH/78wgsvMGrUKCDSAvj2t7/NBRdcQOfOncPWgLszfvx4ioqKuPLKK/noo4/S/K8nqcjKloLWKUgmBWNSdano1BZMGtC92uetXLmSxx57jAsvvJAbbriBhx9+mPHjx/PTn/4UgK997Wu8/PLLDBgwgGuvvZYJEyYwePBg9u7dy6FDh/jwww8BeOedd/jWt77FrFmzKCwsZPz48Vx66aX86Ec/4tVXX2Xq1KlVXvPxxx/n4YcfZtOmTdx2220sXryYk046icsvv5wXX3yRQYMGJY27oqKCt956iw8++ICrr76aIUOGMHPmTFauXMmSJUvYsmULRUVF3HDDDbX4V5S6kJUtBXef7e5jWrZsmelQROpVhw4duPDCCwG47rrreOutt5g7dy6lpaX06NGDN954g2XLlrFnzx42btzI4MGDgcgipmOPPRaAFStWMGbMGGbPnk1hYWSixVtvvcXw4cMB6NevHyeddFL4mqeddhrnnXceEGld9OnTh7Zt29KkSROuvfZa5s+fX23cgwYN4qijjqKoqIgtW7YAMH/+fEaMGEFeXh6nnnoql156aR39K0ltZGVLQSSTUvmLPl1ipxiaGbfccgvvvfceHTp04Pbbb2fv3r0kW5San5/P3r17+dvf/sapp54KkPT5xx13XPhzsudFxxY7N/6YY46Jew9N8214srKlIJKrNmzYwIIFCwCYMWMGF110EQBt2rTh448/DvvrW7RoQUFBAS+++CIAn332GZ988gkAJ554Iq+88goTJ05k3rx5AFx00UU8//zzALz22mv861//ivv6paWl/OUvf2Hbtm0cPHiQGTNmcPHFFwNw8skns2LFCg4dOsTMmTOrfS9f+tKXePbZZzl48CAVFRXMnTv3CP9VpC4pKYhkkW7duvHkk0/Ss2dPduzYwdixYxk9ejQ9evRg0KBBnHPOOeFzn376aR544AF69uzJBRdcwObNm8NzJ598MrNnz2bcuHEsXLiQSZMm8dprr9G7d2/++Mc/kp+fzwknnHDY6+fn53PXXXdxySWXhAPYAwcOBGDy5MlcddVVXHrppeTn51f7XgYPHkzXrl3p0aMHY8eODZOLZFbO1j4a9mjkr63nbjq/LkOSRmrFihV069atXl4rWKNwxhln1MvrQaQlkZeXR5MmTViwYAFjx46lrKys3l5f0ifeZzdZ7SONKYgIGzZsYOjQoRw6dIimTZsybdq0TIckGaKkICJ07dqVv/3tb5kOQxoAjSmIiEhISUFEREINKimY2XFmttjMGmbFLhGRRi6tScHMppvZR2a2NOZ4PzNbaWZrzGxC1KnbgOfTGZOIiCSW7pbCE0C/6ANmlgf8BugPFAEjzKzIzC4DlgNb0hyTiIgkkNak4O7zgR0xh88F1rj7WnffBzwLDAQuAc4DRgKjzaxBdW2JZNL27dspLi6muLiYU045hfbt24eP77jjDrp3707Pnj0pLi5m4cKF4XVDhgxh7dq1lJaWUlxcTGFhIW3btg2vXbJkCaeffjqrV68GYP/+/fTo0SO8x89//vOE9w7s2LGDvn370rVrV/r27Ruuhl6/fj3NmzcPX+vmm28Or1m8eDE9evSgS5cufPvb365S+uL555+nqKiI7t27M3LkyJT/jdatW0dpaSldu3Zl2LBh7Nu3D4jsIdGyZcswjjvvvDO8ZufOnQwZMoQvfvGLdOvWLVwt/oMf/IA33ngj5df+4IMPOP/88znmmGP41a9+FR7fu3cv5557LmeddRbdu3dn0qRJNbo+cPDgQXr16lVlL4zbb7+9yudgzpw5KcebVFAjPV1fQEdgadTjIcBvox5/DXgo6vEo4Kok9xsDvAe8V1hY6Edq6JR3fOiUd474eskty5cvr7fX+uCDD/yDDz5IeH7SpEl+zz33uLv7O++84+edd57v3bvX3d23bt3qGzdudHf3pUuX+qBBg6pc+/jjj/u4ceOqHHvuuee8b9++7u7+i1/8wseMGVPtvaP98Ic/9Lvuusvd3e+66y6/9dZb3d193bp13r1797jv4ZxzzvF33nnHDx065P369fM5c+a4u/uqVau8uLjYd+zY4e7uW7ZsOezaxx9/3CdNmnTY8WuuucZnzJjh7u433XSTP/zww+7uPnfuXL/yyivjxvH1r3/dp02b5u7un332mf/rX/9yd/f169eH/yap2LJliy9atMgnTpwY/rdxdz906JDv2bPH3d337dvn5557ri9YsCDl6wO//vWvfcSIEVXeR/TnIJl4n13gPU/wOzYT6xTiVcAK/0xw9yeSXezuU4GpEFnRXKeRiaSoT58+dXq/oAZRTVVUVNCmTZuw4FybNm3Cc7///e/DEhTJDB06lOnTp3P33XczZcqUcL1CsntHmzVrVhj/9ddfT58+ffjlL3+ZNObdu3dz/vmRagJf//rXefHFF+nfvz/Tpk1j3LhxYZXWdu3aVRs/RP64feONN3jmmWfCOG6//XbGjh2b8Jrdu3czf/58nnjiCQCaNm1K06ZNgUhl2O3bt7N582ZOOeWUal+/Xbt2tGvXjldeeaXKcTML95nYv38/+/fvj1sEMNH1AOXl5bzyyiv853/+J/fee2+1sdRWJrpoyoEOUY8LgE01uYGZDTCzqbt27arTwESyzeWXX86HH37IF77wBW655Rb+8pe/hOfefvttzj777JTuc//993Pbbbfx4x//mFatWlV772hbtmwJax3l5+dX2Sxn3bp19OrVi4svvpg333wTgI0bN1JQUBA+p6CggI0bI1vjrlq1ilWrVnHhhRdy3nnn8eqrr6YU//bt2znxxBNp0qTJYfcEWLBgAWeddRb9+/dn2bJlAKxdu5a2bdvyjW98g169evHNb36Tf//73+E1vXv35u233wbge9/7XthNE/01efLkamM7ePAgxcXFtGvXjr59+1JaWprSewp897vf5e677467491DDz1Ez549ueGGGxIWMaypTLQU3gW6mlknYCMwnMg4QsrcfTYwu6SkZHQa4hOp1pH+ZV/Xjj/+eBYvXsybb77J3LlzGTZsGJMnT2bUqFFUVFTQtm3blO7z6quvkp+fz9Kln08UTHbvVOTn57NhwwZat27N4sWLGTRoEMuWLYtbfjv46/nAgQOsXr2aefPmUV5ezn/8x3+wdOlSDh48yJe//GUgMoaxb9++sALs008/Hfev+eCevXv35p///CfHH388c+bMYdCgQaxevZoDBw7w/vvv8+CDD1JaWsp3vvMdJk+ezM9+9jMg8tf7pk2Rv1fvu+++lN5zPHl5eZSVlbFz504GDx7M0qVLU95r++WXX6Zdu3acffbZh33mxo4dy09+8hPMjJ/85Cd8//vfZ/r06UccZyDdU1JnAAuAM8ys3MxudPcDwHjgT8AK4Hl3X1bD+6qlIFIpLy+PPn36cMcdd/DQQw/xhz/8AYDmzZsftq9BPJs2beKBBx5g0aJFzJkzh7///e/V3jvaySefTEVFBRDpGgq6fI455hhat24NwNlnn83pp5/OqlWrKCgooLy8PLy+vLw83NehoKCAgQMHcvTRR9OpUyfOOOMMVq9eTevWrSkrK6OsrIw777yTm2++OXzco0cP2rRpw86dOzlw4MBh92zRokXYhfOVr3yF/fv3s23bNgoKCigoKAj/ch8yZAjvv/9+GNfevXtp3rw5ULuWQuDEE0+kT58+Kbd+INLae+mll+jYsSPDhw/njTfe4Lrrrgv/3fPy8jjqqKMYPXo0ixYtSvm+yaR79tEId89396PdvcDdH6s8Psfdv+Dup7v7z4/gvtp5TYRIRdVg5hBAWVkZp512GhAps71mzZpq7/G9732PiRMnUlBQwL333su4ceNw96T3jnb11Vfz5JNPAvDkk0+G4xhbt27l4MGDQKSrZvXq1XTu3Dksy/3Xv/4Vd+epp54Krxk0aFC4r8K2bdtYtWoVnTt3rvY9mBmXXHJJuJ9EdBybN28OWyeLFi3i0KFDtG7dmlNOOYUOHTqEVWlff/11ioqKwnuuWrUq/Iv+vvvuC5NQ9NeECdHLrA63detWdu7cCcCnn37K//zP//DFL36x2vcTuOuuuygvL2f9+vU8++yzXHrppfzud78DCBMxwMyZM1NufVRHBfFEstjHH3/Mt771LXbu3EmTJk3o0qVLuL/ylVdeybx587jssssSXv/nP/+ZDRs2cOONNwIwYMAApk2bxlNPPcWZZ56Z8N7RJkyYwNChQ3nssccoLCzkv/7rv4DIdps//elPadKkCXl5eUyZMiUcr3jkkUcYNWoUn376Kf3796d///4AXHHFFbz22msUFRWRl5fHPffcE7Y2qvPLX/6S4cOH8+Mf/5hevXqF7+mFF17gkUceoUmTJjRv3pxnn3027Fp68MEHufbaa9m3bx+dO3fm8ccfByKDwmvWrKGkJG516cNs3ryZkpISdu/ezVFHHcX999/P8uXLqaio4Prrr+fgwYMcOnSIoUOHhtNKp0yZAsDNN9+c8PoWLVokfM1bb72VsrIyzIyOHTvy6KOPphRrdbJyPwUzGwAM6NKly+jov2RqQvspSE1k434Kn376KZdccglvv/02eXl5dRFazpg5cybvv/9+OL6QzWq6n0JWLhBT95FI9Zo3b84dd9xRZRaOpObAgQN8//vfz3QYGaHuI5EUuXvWbTR/xRVXZDqErHTNNddkOoQ6cSQ9QVnZUtDsI6lvzZo1Y/v27Uf0P5lIJrg727dvp1mzZjW6LitbClqnIPUtmEa5devWtL/W5s2bATh06FDaX0sat2bNmlVZKJiKrEwKIvUtmDdfH4LSDA1lgZzkFnUfiYhIKCuTgmYfiYikR1YmBRERSQ8lBRERCSkpiIhIKCuTggaaRUTSIyuTggaaRUTSIyuTgoiIpIeSgoiIhJQUREQklJVJQQPNIiLpkZVJQQPNIiLpkZVJQURE0kNJQUREQkoKIiISUlIQEZGQkoKIiISUFEREJJSVSUHrFERE0iMrk4LWKYiIpEdWJgUREUkPJQUREQkpKYiISEhJQUREQkoKIiISUlIQEZFQTieFhet28MzCDZkOQ0SkwcjZpDCwuD0AE2cuUWIQEamUs0lhZGkhvxjcA4BZZRszHI2ISMPQYJKCmXUzsylm9oKZja2P1xxZWkhpp1Ysr9jNsEcXqMUgIjmvSTpvbmbTgauAj9z9zKjj/YD/B+QBv3X3ye6+ArjZzI4CpqUzrmhBN9LCdTtYuG4HEEkWIiK5KN0thSeAftEHzCwP+A3QHygCRphZUeW5q4G3gNfTHFdoZGkhz910ftiVpDEGEcllaW0puPt8M+sYc/hcYI27rwUws2eBgcByd38JeMnMXgGeiXdPMxsDjAEoLKy7v+iD1sHEmUuYOHMJs8o28sm+gxzbNI+Bxe3VehCRnJDWpJBAe+DDqMflQKmZ9QG+ChwDzEl0sbtPBaYClJSUeF0GFvzin1W2keUVu9mz98Bh50REGrNMJAWLc8zdfR4wL6UbmA0ABnTp0qUOw4oYWVrIyNJChj26IBxjEBHJFZmYfVQOdIh6XABsqskN6mM/hYHF7Snt1ArQIjcRyR3VJgUzu9DM/mxmq8xsrZmtM7O1tXjNd4GuZtbJzJoCw4GXanKD+th5Ld4AtKatikhjl0r30WPA94DFwMGa3NzMZgB9gDZmVg5McvfHzGw88CciU1Knu/uymtzX3WcDs0tKSkbX5LojET3OEExbnbFoA8c2zQPQILSINCqpJIVd7v7HI7m5u49IcHwOSQaTG5pgnOGZhRuYOHMJSzZGWignNGsSnhcRaQxSSQpzzewe4L+Bz4KD7v5+2qKqRjoHmpOJbjUMLG4fth6eWbhBiUFEGoVUkkJp5feSqGMOXFr34aSmPruPYgWthsDCdTuYVbZRSUFEGoVqk4K7X1IfgWSjkaWFYTG9ZxZuCH/WOIOIZKuEScHMrnP335nZ/4133t3vTV9YyWWq+yiR5RW7wzUNGmcQkWyWrKVwXOX3E+ojkJrIZPdRrKCgXvCzWg4iks0SJgV3f7Ty+x31F072iR1jCEpkqOUgItkolcVrnc1stpltNbOPzGyWmXWuj+Cy0cDi9hTlt6C0Uyt+MbgHRfkttF+DiGSNVGYfPUOk1PXgysfDgRl8Piup3jW0MYVosS2HwPKK3eF5EZGGKpXaR+buT7v7gcqv3xGZkpox9VH7qK4E5TKK8luohpKINHgJk4KZtTKzVkQWr00ws45mdpqZ3Qq8Un8hNg7BgLT2gxaRhixZ99FiIi2CoNT1TVHnHPhZuoJqjII1DcH4gmYkiUhDlGz2Uaf6DKQmGvKYQjKx+0EH5TKUHESkocjEfgq1lk1jCtGiy3GXdmrF8ordVbqTnlm4QbOURCSjsjIpZLvowedAUIE1aEGIiGRC0impZmZAgbt/mOx5cuSCMQZt/SkiDUHSpODubmYvAmfXUzw5JbpERmmnVlXKZIiIZEIqi9f+ambnuPu7aY8mx8Rb6KYZSiKSSakkhUuAm81sPfBvIlNU3d17pjOwZLJ19lEq4s1Qij2vRCEi6ZJKUuif9ihqqCFVSa1r0Vt/xiYElcoQkXRLZZOdf5rZRUBXd3/czNoCx6c/tNwWr2tp2KMLMhSNiOSKVKqkTgJuA35Ueeho4HfpDEoSU/0kEUmnVLqPBgO9gPcB3H2TmTW4jXdywcDi9ixct4OJM5doAx8RSYtUksK+yqmpDmBmx1V3gaRHUD8pGISGSMthxqINHNs0D1CSEJHaSSUpPG9mjwInmtlo4AZgWnrDkkSC2UnB94kzl7Bk465whzfQQLSIHLlUBpp/ZWZ9gd3AF4Cfuvuf0x6ZxBVv+8+F63ZUKZkhInKkUq19tAR4E5hf+XNGmdkAM5u6a9euTIeScQOL24eroUED0SJSO6nMPvomsAj4KjCEyArnG9IdWDLZWiU1HYLieiNLC6t0KanaqogciVTGFH4I9HL37QBm1hp4B5iezsCk5oJupejBaO3ZICI1kUpSKAf2RD3eA6hqagMVuyJaq6BFpCZSSQobgYVmNovINpwDgUVm9n8B3P3eNMYnRyhIDloFLSI1kUpS+EflV2BW5XctYBMRaWRSmZJ6R30EIumzcN2Ow1oMGmcQkXi0HWeOiN7ZLXZvaBGRQCrdR5LFoldABy0DjTOISCJKCo1cvBLcgHZ3E5G4Ulm8dreZtTCzo83sdTPbZmbXpSMYMxtkZtPMbJaZXZ6O15BIq6Eov0VYcVUL3UQkkMqYwuXuvhu4isiahS8QWdCWEjObbmYfmdnSmOP9zGylma0xswkA7v6iu48GRgHDUn0NqZlgFfQvBvegtFMrjTGISCiVpHB05fevADPcfUeyJ8fxBNAv+oCZ5QG/IbLVZxEwwsyKop7y48rzkkZBcijKbxF2J6nFIJLbUhlTmG1mHwCfArdUbse5N9UXcPf5ZtYx5vC5wBp3XwtgZs8CA81sBTAZ+KO7vx/vfmY2BhgDUFiovvC6EF1MLyiNEX1OYw4iuaPaloK7TwDOB0rcfT/wbyKrmmujPVVLZZRXHvsWcBkwxMxuThDPVHcvcfeStm3b1jIMgcO7kwLR3UrPLNygloRIDkjYUjCzr8Y5Fv3wv2vxuhbnmLv7A8AD1V5sNgAY0KVLl1qEILFiZypFT10NiuwFzxORxilZ99GAJOec2iWFcqBD1OMCYFOqF7v7bGB2SUnJ6FrEICIiMRImBXf/Rhpf912gq5l1IlJwbzgwMo2vJ0coKJERVFsVkcYtpcVrZnYl0B1oFhxz9ztTvHYG0AdoY2blwCR3f8zMxgN/AvKA6e6+LNWg1X1Uv5ZX7A7XNYhI41ZtUjCzKcCxwCXAb4nsvrYo1Rdw9xEJjs8B5qR6n5hr1X1Uj4ryW/DcTecz7NEFVYrraWaSSOOTSkvhAnfvaWZ/d/c7zOzX1G48odbUUqgf0XWTogUtBg08izQ+qSSFTyu/f2JmpwLbgU7pC6l6ainUj9jZSLFJYuLMJUycuSR8rohkv1SSwstmdiJwD/A+kZlHv01rVNIgxSuuN3HmEmaVbayyBWjQrRT7WEQavlQ22flZ5Y9/MLOXgWbuviu9YSWn7qOGYWRpYbgPdDDeEHtOXUwi2SXZ4rVL3f2NRIvY3D1j4wrqPmo4oscbSju1CktlaOWzSHZK1lK4GHiD+IvYart4TRqJ2C6lZxZuCMcaTmim7TpEsk2yxWuTzOwoIsXpnq/HmCSLBQkiqJkU3XJQF5JIw5e0IJ67HwLG11MsKTOzAWY2ddeujA5tSAJBgb2gyB6g/RpEskQq+yn82cx+YGYdzKxV8JX2yJJw99nuPqZly5aZDENSMLK0MBxrCMYZVHFVpOFKpdP3hsrv46KOOdC57sORxmhgcftwn4boWUnRezdo2qpIw5BKUujm7lU21TGzZomeLBIr0dTVYO+GoNhebFLQOgeR+pdKUngH6J3CsXqjdQrZJ5i6GiSA0k6teO6m84Gq+zZE0zoHkfqXcEzBzE4xs7OB5mbWy8x6V371IVIgL2M0ppB9oveDjicotKdxBpHMStZSuAIYRWQDnF/z+W5pu4GJ6Q1LcpFaBSKZl2ydwpPAk2b2f9z9D/UYk+SYoGtJ+zWIZF4qU1K7BT+Y2TFpjEVyVNC1FAw8i0jmJBtTuNXMzieyqU4g/oigSB0JZihpbEEkM5KNKawErgE6m9mbwAqgtZmd4e4r6yW6BDT7qHGKnqGUqCspmKYa/XxNWxWpO8mSwr+IDCj3qfzqRmTweUJlYrgg7dEloCqp2SvRbm7weXG92DLc0YL1DtE0QC1Sd5IlhX7AJOB04F7gf4F/u/s36iMwaZzibdSTSOwYwzMLN7Bw3Y6EYw9a7CZSewnHFNx9ort/GVgP/I5IAmlrZm+Z2ex6ik8kFNttFO98dOkMEam5VGYf/cnd33X3qUC5u18EqLUg9Sa6PEZpp1ZqBYikUSrbcd4a9XBU5bFt6QpIJFpseYxoWtcgUvdSaSmE3P1/0xWISDzVlceIJyiZoamtIjWn/RIlq8VWWg1EP47dLlSD0SKJZWVS0DqF3JNsKmv0saBLKVHLQpVXRZLLyqSgdQqNW7wEkGwqa/S5ZGscRKR6WZkUpHGryVqGuqAuJZHPKSlIoxYMOif7ha/tQUU+p6QgWSnRGEP08eAXfE3GEBINUIvkCiUFyUqJupiij0evbF64bkdK01NrMvVVpDGq0ToFkWwysLg9pZ1a0aN9ZNvWiTOXxF0EFyjt1CphCQ2RXKGWgjRa0a2G6JLb8WYnlXZqxXM3nQ+g2kmS05QUJCfETlsVkfiUFEQSiG5dgGYkSW5oMGMKZtbZzB4zsxcyHYsIVN3QZ3nFbnUrSU5Ia1Iws+lm9pGZLY053s/MVprZGjObAODua939xnTGI1JTRfktwoJ8wZoHFdmTxizdLYUniOzgFjKzPOA3QH+gCBhhZkVpjkOkTmgTH2ns0jqm4O7zzaxjzOFzgTXuvhbAzJ4FBgLL0xmLSLTgr36IdA3Frk8IZigl2vozED3uoDEHaQwyMabQHvgw6nE50N7MWpvZFKCXmf0o0cVmNsbM3jOz97Zu3ZruWKURC8YLivJbVFmfUJO1CsG4g8YcpLHIxOwji3PM3X07cHN1F1duCzoVoKSkxOs4NskhwXhBrJGlhVVKbKdyH5HGIhNJoRzoEPW4ANhUkxtoPwVJt2T7N4g0ZpnoPnoX6GpmncysKTAceKkmN3D32e4+pmXLlmkJUCTYBrQmYwSp1lcSacjSPSV1BrAAOMPMys3sRnc/AIwH/gSsAJ5392U1vO8AM5u6a9euug9a5AgELQqNK0i2S/fsoxEJjs8B5tTivtp5TRqUYBxCJNupzIXIEQimoi5ctyPptNXoXd0ATV+VBi8rk4IGmqW+Bb/UgxlJqc5Oin2eNvGRhq7B1D6qCQ00S30LBp6rW8xWnaL8FprCKg1aVrYURDIp2UY90XtCB49FsklWJgV1H0lt1GYNQvQ1iX7hKxFINsvKpKDZR1IbifZ3rum1wx5dUG0CqG13k0h9y8oxBRERSY+sTApavCYN0cDi9moZSNbLyqSg2UfSENXVDCWRTMrKpCAiIumhpCAiIqGsnH0kkm2iN/SBqju/fbLvIMc2zSwngxMAAAc+SURBVAufu6ddT4DwfEClMaQ+ZGVS0DoFaQhSXe8Qu6tbUP8oSBR79h7ghGZNKMpvwfKK3exr0y08HySR4LlKCpJuWZkUtE5BGoJU1zvEPi9ICtGthmAXuGGPLqDsI8Lzwc5wsa0GkXTRmIKIiISUFEREJJSV3UciuS7YpyF6kLq6gejovR1qOjZRm2slu2RlUtBAs+S66H0aTmj2+f/GyX5hR19T01/stblWsktWdh9pRbPI57RHg9SlrEwKIiKSHkoKIiISUlIQEZGQkoKIiISUFEREJJSVU1JFGrrYAng18VmL+FM+o4voBfdPdD6QyrqC6DUIgNYj5LisTApapyANWWwBvETn4z1vYHH7pPs+B0XyivJbhAvXYovs1bSIXvQaBEDrEXJcViYFFcSThqy6Qnmx52N/vuPJVxK2FqKL5EWLLrKnInpSGxpTEBGRkJKCiIiElBRERCSkpCAiIiElBRERCSkpiIhISElBRERCSgoiIhJSUhARkVCDWdFsZscBDwP7gHnu/vsMhyQiknPS2lIws+lm9pGZLY053s/MVprZGjObUHn4q8AL7j4auDqdcYmISHzpbik8ATwEPBUcMLM84DdAX6AceNfMXgIKgCWVTzuY5rhEckJQOTW6+F5sNdWgcF505dXlFburPCe2SF90ZdWRpYXh42hBwb5AKlVYY+8b71xwr5oU7AuujS4iGFx/pPetTTy1Ebxu0aktmDSge53fP61Jwd3nm1nHmMPnAmvcfS2AmT0LDCSSIAqAMpK0YMxsDDAGoLBQVRyl8Zl0/ZVVSllD/Mqq0ZKdj1d1NbqaanTFVTj8l3m8Mt3RlVVHlhYyq2zjYRVa9+w9wAnNmlCU36LKPZJVYY29b+y56PvU5JdwbCXY6OuP9L61iac2gvdSdGrNy7KnIhNjCu2BD6MelwOlwAPAQ2Z2JTA70cXuPhWYClBSUuJpjFMkI+JVWa1p5dXqJKq2Gs+wRxckLecd757BNcGxuqrYeiT7U6TzvumKpzqlnVqlpZUAmUkKFueYu/u/gW+kdAPtpyAikhaZmJJaDnSIelwAbKrJDdx9truPadmyZZ0GJiKS6zKRFN4FuppZJzNrCgwHXqrJDcxsgJlN3bVrV1oCFBHJVemekjoDWACcYWblZnajux8AxgN/AlYAz7v7sprcVy0FEZH0SPfsoxEJjs8B5qTztUVEpOayssyFuo9ERNIjK5OCuo9ERNIjK5OCiIikh7ln7/ovM9sK/PMIL28DbKvDcLKN3r/ev95/7jrN3dvGO5HVSaE2zOw9dy/JdByZovev96/3n7vvPxl1H4mISEhJQUREQrmcFKZmOoAM0/vPbXr/ElfOjimIiMjhcrmlICIiMZQUREQklHNJIcH+0DnDzNab2RIzKzOz9zIdT32It1e4mbUysz+b2erK7ydlMsZ0SvD+bzezjZWfgzIz+0omY0wnM+tgZnPNbIWZLTOz71Qez5nPQE3kVFKI2h+6P1AEjDCzosxGlRGXuHtxDs3TfgLoF3NsAvC6u3cFXq983Fg9weHvH+C+ys9BcWWRysbqAPB9d+8GnAeMq/z/Ppc+AynLqaRA1P7Q7r4PCPaHlkbM3ecDsftJDgSerPz5SWBQvQZVjxK8/5zh7hXu/n7lz3uIlOxvTw59Bmoi15JCvP2h4++E3ng58JqZLTazMZkOJoNOdvcKiPzSANplOJ5MGG9mf6/sXsqJrhMz6wj0Ahaiz0BcuZYU4u4PXe9RZNaF7t6bSBfaODP7UqYDkox4BDgdKAYqgF9nNpz0M7PjgT8A33X33ZmOp6HKtaRQ6/2hs527b6r8/hEwk0iXWi7aYmb5AJXfP8pwPPXK3be4+0F3PwRMo5F/DszsaCIJ4ffu/t+Vh3P6M5BIriWFWu8Pnc3M7DgzOyH4GbgcWJr8qkbrJeD6yp+vB2ZlMJZ6F/wyrDSYRvw5MDMDHgNWuPu9Uady+jOQSM6taK6cenc/kAdMd/efZzikemNmnYm0DiCyFeszufD+K/cK70OkXPIWYBLwIvA8UAhsAK5x90Y5GJvg/fch0nXkwHrgpqB/vbExs4uAN4ElwKHKwxOJjCvkxGegJnIuKYiISGK51n0kIiJJKCmIiEhISUFEREJKCiIiElJSEBGRUJNMByCSzcysNZFiagCnAAeBrZWPZwJDK48dIjLtc2G9BylSA5qSKlJHzOx24GN3/5WZnQ/cC/Rx98/MrA3QNFhRLtJQqaUgkh75wDZ3/wzA3bdlOB6RlGhMQSQ9XgM6mNkqM3vYzC7OdEAiqVBSEEkDd/8YOBsYQ2SM4TkzG5XRoERSoO4jkTRx94PAPGCemS0hUnTtiUzGJFIdtRRE0sDMzjCzrlGHioF/ZioekVSppSCSHscDD5rZiUT2CF5DpCtJpEHTlFQREQmp+0hEREJKCiIiElJSEBGRkJKCiIiElBRERCSkpCAiIiElBRERCf1/8bTF09INVrAAAAAASUVORK5CYII=",
      "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}