{ "cells": [ { "cell_type": "markdown", "id": "6af675ed", "metadata": {}, "source": [ "# Hypothesis tests with pyhf\n", "\n", "This notebook will provide you with the tools to do sensitivity estimates which can be used for search region optimization or sensitivity projections.\n", "\n", "## p-value for discovery of a new signal\n", "\n", "In searches for new physics we want to know how significant a potential deviation from our Standard Model (SM) expectation is. We do this by a hypothesis test where we try to exclude the SM (\"background only\") hypothesis. We use a so called **p-value** $p_0$ for this, abstractly defined by:\n", "\n", "$$p_0 = \\int\\limits_{t_\\mathrm{obs}}^{\\infty}p(t|H_0)\\mathrm{d}t$$\n", "\n", "where $t$ is a test statistic (a number we calculate from our data observations) and $p(t|H_0)$ is the probability distribution for $t$ under the assumption of our **null Hypothesis** $H_0$, in this case the background only hypothesis. This p-value is then typically converted into a number of standard deviations $z$, the **significance** (\"number of sigmas\") via the inverse of the cumulative standard normal distribution $\\Phi$:\n", "\n", "$$z = \\Phi^{-1}(1 - p)$$\n", "\n", "The typical convention for particle physics is to speak of *evidence* when $z>3$ and of an *observation* when $z>5$.\n", "\n", "So what do we use for $t$? We want to use something that discriminates well between our null Hypothesis and an **alternative Hypothesis** that we have in mind. When we try to discover new physics, our null Hypothesis is the absence and the alternative Hypothesis the presence of a signal. We can parametrize this by a **signal strength** parameter $\\mu$. The test statistics used in almost all LHC searches use the **profile likelihood ratio**\n", "\n", "$$\\Lambda_\\mu = \\frac{L(\\mu, \\hat{\\hat{\\theta}})}{L(\\hat{\\mu}, \\hat{\\theta})}$$\n", "\n", "where $\\theta$ are the other parameters of our model that are not part of the test, the so called **nuisance parameters**. In contrast, the parameter that we want to test, $\\mu$, is called our **parameter of interest** (POI). The nuisance parameters include all fit parameters, like normalization factors and parameters for describing uncertainties. $L(\\mu, \\hat{\\hat{\\theta}})$ is the Likelihood function, maximized under the condition that our parameter of interest takes the value $\\mu$ and $L(\\hat{\\mu}, \\hat{\\theta})$ is the unconditionally maximized Likelihood. So roughly speaking, we are calculating the fraction of the maximum possible likelihood that we can get under our test condition. If it is high, that speaks for our hypothesis, if it is low, against. The test statistic $t_\\mu$ is then defined as\n", "\n", "$$t_\\mu = -2\\ln\\Lambda_\\mu$$\n", "\n", "giving us a test statistic where **high values speak against the null hypothesis**.\n", "\n", "
\n", " Question 7a: If we want to discover a new signal (using the p-value $p_0$), which value of $\\mu$ are we testing against? Or in other words, what is our null Hypothesis?\n", "
\n", "\n", "All that's left now is to know the distribution of $p(t_\\mu|H_0)$. [Wilk's theorem](https://en.wikipedia.org/wiki/Wilks%27_theorem) tells us that the distribution of $t_\\mu$ is asymptotically (for large sample sizes) a chi-square distribution. For the discovery p-value we use a slightly modified version of test statistic, called $q_0$ where $\\hat{\\mu}$ is required to be $>=0$ ($q_0=0$ for $\\hat{\\mu} < 0$). For $q_0$ the p-value in the asymptotic limit collapses to a very simple formula:\n", "\n", "$$p_0 = \\sqrt{q_0}$$\n", "\n", "The asymptotic limit often matches quite well even for fairly small sample sizes, but it should be kept in mind this is an approximation. Alternatively, one can evaluate $p(t_\\mu|H_0)$ by Monte Carlo sampling (\"toys\").\n", "\n", "## CLs for exclusion of an absent signal\n", "\n", "Now, sadly, not all searches find evidence for new physics. What we still can do in such a case is to try exclude models by rejecting the hypothesis of a signal being present. That usually means we test against $\\mu=1$ or some other value $>0$. The rest of the procedure is very similar with one small detail worth mentioning ... In high energy physics it is very common to use a quantity called $CL_s$ instead of plain p-value. It is defined by\n", "\n", "$$CL_s = \\frac{CL_{s+b}}{CL_{b}}$$\n", "\n", "where $CL_{s+b}$ is the p-value for rejecting the hypothesis of signal + background being present (what would be the \"normal\" p-value) and $CL_{b}$ is the p-value for rejecting the background only hypothesis, but now using the test statistic for $\\mu=1$ (so this is different from $p_0$!). We won't go into further details how to calculate those p-values. `pyhf` has the formulas included and does it automatically for us. The asymptotic distributions for all different variants are described in the paper \"Asymptotic formulae for likelihood-based tests of new physics\" ([arXiv:1007.1727](https://arxiv.org/abs/1007.1727)).\n", "\n", "Just a qualitative explanation of why we use $CL_s$ instead of the p-value: We want to avoid excluding signals in cases where we don't have sensitivity, but observe an *underfluctuation* of the data. In these cases $CL_{s+b}$ and $CL_b$ will be very similar and consequently lead to a large value of $CL_{s}$, telling us the signal is **not** excluded. In case our observations are exactly on spot with the background expectations $CL_b = 0.5$ in the asymptotic limit, so on average we have twice as high \"p-values\" with $CL_s$.\n", "\n", "The typical convention for particle physics is to speak of an **exclusion** of a signal if $CL_s < 0.05$. That's usually what is meant by \"limit at 95% confidence level\"." ] }, { "cell_type": "markdown", "id": "59665ed7", "metadata": {}, "source": [ "## Discovery or exclusion of a signal for a cut & count experiment" ] }, { "cell_type": "markdown", "id": "fbc8e9e6", "metadata": {}, "source": [ "Let's start with a simple case where we only want to count the number of events in a certain search region. We assume a certain number of expected background events `b`, expected signal events `s` and a total uncertainty on the expected background `delta_b` ($\\sigma_b$).\n", "\n", "The likelihood function for this can be formulated as a primary measurement of `n` events and a control (\"auxiliary\") measurment of `m` events that constrains our background parameter within the uncertainty. So, a product of 2 Poisson distributions:\n", "\n", "$$L(s, b) = \\mathrm{Pois}(n|s + b)\\cdot \\mathrm{Pois}(m|\\tau b)$$\n", "\n", "The parameter $\\tau$ can be given in terms of $\\sigma_b$ by asking the question \"How much more events do i have to measure in the control region to get the relative uncertainty $\\sigma_b / b$\". That gives\n", "\n", "$\\tau = \\frac{b}{\\sigma_b^2}$\n", "\n", "Equivalently, we can replace $b$ by $\\gamma b$ and $s$ by $\\mu s$ to fit normalization factors (initialized to 1) and keep $s$ and $b$ fixed to our expectation.\n", "\n", "$$L'(\\mu, \\gamma) = L(\\mu s, \\gamma b)$$\n", "\n", "`pyhf` has a convenience function to create the specification for such a model: `pyhf.simplemodels.hepdata_like`. It also works for arbitrary many bins, but for now let's go with one bin and 5 expected background events, 7 expected signal events and an uncertainty of 2 on the expected background events:" ] }, { "cell_type": "code", "execution_count": null, "id": "29ba6b01", "metadata": {}, "outputs": [], "source": [ "import pyhf\n", "from scipy import stats\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "id": "228261b1", "metadata": {}, "outputs": [], "source": [ "s = 7\n", "b = 5\n", "delta_b = 2" ] }, { "cell_type": "code", "execution_count": null, "id": "5d65cb13", "metadata": {}, "outputs": [], "source": [ "model = pyhf.simplemodels.hepdata_like(\n", " signal_data=[s], bkg_data=[b], bkg_uncerts=[delta_b]\n", ")" ] }, { "cell_type": "markdown", "id": "fd0168e3", "metadata": {}, "source": [ "The model comes with a \"parameter of interest\" (POI) called `mu` that is our signal strength:" ] }, { "cell_type": "code", "execution_count": null, "id": "c9a435b3", "metadata": {}, "outputs": [], "source": [ "model.config.poi_name" ] }, { "cell_type": "markdown", "id": "403d454b", "metadata": {}, "source": [ "In addition, we have one nuisance parameter, the constrained background normalization $\\gamma$, called `uncorr_bkguncrt` here:" ] }, { "cell_type": "code", "execution_count": null, "id": "eea9a3a6", "metadata": {}, "outputs": [], "source": [ "model.config.par_order" ] }, { "cell_type": "markdown", "id": "13498031", "metadata": {}, "source": [ "It's initial value should be 1" ] }, { "cell_type": "code", "execution_count": null, "id": "535f5314", "metadata": {}, "outputs": [], "source": [ "gamma_initial = 1" ] }, { "cell_type": "markdown", "id": "6d6d4528", "metadata": {}, "source": [ "So the expected data in our model scales with `mu`. For `mu=1` we get `5 * 1 * 7 = 12`" ] }, { "cell_type": "code", "execution_count": null, "id": "c27bd580", "metadata": {}, "outputs": [], "source": [ "model.expected_actualdata([1, gamma_initial])" ] }, { "cell_type": "markdown", "id": "de757850", "metadata": {}, "source": [ "for `mu=2` we get `5 + 2 * 7 = 19`" ] }, { "cell_type": "code", "execution_count": null, "id": "e30c4ab7", "metadata": {}, "outputs": [], "source": [ "model.expected_actualdata([2, gamma_initial])" ] }, { "cell_type": "markdown", "id": "c823d879", "metadata": {}, "source": [ "The auxiliary data corresponds to $\\tau b$ in the formula above:" ] }, { "cell_type": "code", "execution_count": null, "id": "041e884e", "metadata": {}, "outputs": [], "source": [ "model.config.auxdata" ] }, { "cell_type": "markdown", "id": "484c7fb9", "metadata": {}, "source": [ "It's given by our background uncertainty `delta_b`:" ] }, { "cell_type": "code", "execution_count": null, "id": "22f7a85c", "metadata": {}, "outputs": [], "source": [ "b ** 2 / (delta_b ** 2)" ] }, { "cell_type": "markdown", "id": "c69d76d2", "metadata": {}, "source": [ "To get the p-value for rejection of the background only hypothesis, we call `pyhf.infer.hypotest` with the test value 0 of our POI $\\mu$ using the `q0` test statistic.\n", "\n", "We want to know which p-value we would get if we would observe an excess of events of precisely the expected signal, so we plug in `s + b` for the data:" ] }, { "cell_type": "code", "execution_count": null, "id": "4d42e783", "metadata": {}, "outputs": [], "source": [ "pvalue = pyhf.infer.hypotest(\n", " poi_test=0,\n", " data=[s + b] + model.config.auxdata,\n", " pdf=model,\n", " test_stat=\"q0\"\n", ")\n", "pvalue" ] }, { "cell_type": "markdown", "id": "69482321", "metadata": {}, "source": [ "We can convert this into a significance (number of standard deviations) using the inverse of the cumulative standard normal distribution $\\Phi$\n", "\n", "$$z = \\Phi^{-1}(1 - p)$$\n", "\n", "The function [`scipy.stats.norm.isf`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm) (\"inverse survival function\") calculates $\\Phi^{-1}(1 - p)$ in a numerically stable way (also for small p-values)." ] }, { "cell_type": "code", "execution_count": null, "id": "78ca9b05", "metadata": {}, "outputs": [], "source": [ "def pvalue_to_significance(pvalue):\n", " return stats.norm.isf(pvalue)" ] }, { "cell_type": "code", "execution_count": null, "id": "6ad12b0a", "metadata": {}, "outputs": [], "source": [ "pvalue_to_significance(pvalue)" ] }, { "cell_type": "markdown", "id": "8ebf754d", "metadata": {}, "source": [ "That would not count as \"Evidence\" yet.\n", "\n", "
\n", " Exercise 7b: How much excess events would we need to observe in our search region (assuming unchanged expected background) that we have potential for finding evidence (3 $\\sigma$) of a new signal?\n", "
" ] }, { "cell_type": "markdown", "id": "35a7c630", "metadata": {}, "source": [ "Equivalently we can test for exclusion and calculate $CL_s$. For that we use 1 as the test value for $\\mu$ and the `qtilde` test statistic.\n", "\n", "We want to know if we could exclude a signal if we would not observe any more data than our background expectation, so we set our data to `b`:" ] }, { "cell_type": "code", "execution_count": null, "id": "38991cf2", "metadata": {}, "outputs": [], "source": [ "CLs = pyhf.infer.hypotest(\n", " poi_test=1,\n", " data=[b] + model.config.auxdata,\n", " pdf=model,\n", " test_stat=\"qtilde\"\n", ")\n", "CLs" ] }, { "cell_type": "markdown", "id": "7a54c716", "metadata": {}, "source": [ "
\n", " Question 7c: Would that signal count as excluded?\n", "
" ] }, { "cell_type": "markdown", "id": "2540c262", "metadata": {}, "source": [ "
\n", " Tip: For this simple number counting scenario, often used for optimizing search sensitivity you don't need to spin up pyhf to calculate the expected significance or CLs. The maximum likelihood estimates can be calculated exactly. Have a look at numbercounting.py for some functions to help you with that.

For a function that calculates the expected significance in one step as a formula have a look at http://www.pp.rhul.ac.uk/~cowan/stat/medsig/medsigNote.pdf\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "35a055cd", "metadata": {}, "outputs": [], "source": [ "from numbercounting import z0_exp, cls_exp" ] }, { "cell_type": "code", "execution_count": null, "id": "0648aabe", "metadata": {}, "outputs": [], "source": [ "z0_exp(s, b, delta_b)" ] }, { "cell_type": "code", "execution_count": null, "id": "f101cab3", "metadata": {}, "outputs": [], "source": [ "cls_exp(s, b, delta_b)" ] }, { "cell_type": "markdown", "id": "d682cdd6", "metadata": {}, "source": [ "These functions support numpy arrays and you can run these functions easily to test millions of different values!" ] }, { "cell_type": "code", "execution_count": null, "id": "5bbd8344", "metadata": {}, "outputs": [], "source": [ "N = 1_000_000\n", "s_random = 20 * np.random.rand(N)\n", "b_random = 20 * np.random.rand(N)\n", "db_rel_random = np.random.rand(N)" ] }, { "cell_type": "code", "execution_count": null, "id": "a79c17fe", "metadata": {}, "outputs": [], "source": [ "%time z0_exp(s_random, b_random, db_rel_random * b_random)" ] }, { "cell_type": "code", "execution_count": null, "id": "941742db", "metadata": {}, "outputs": [], "source": [ "%time cls_exp(s_random, b_random, db_rel_random * b_random)" ] }, { "cell_type": "markdown", "id": "0b1ea2e8", "metadata": {}, "source": [ "## Run an upper limit scan on the signal strength" ] }, { "cell_type": "code", "execution_count": null, "id": "b15741cb", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "4d08b81f", "metadata": {}, "source": [ "Often does not only quote limits on a particular assumed signal strength, but gives a limit on the signal strength itself. This is especially interesting for single-bin (cut & count) search regions, since it is quite model independent (everybody can simulate their own model and calculate the number of excess events in a certain search region to determine if it is excluded by such a limit).\n", "\n", "To do that, we need to invert the hypothesis test, e.g find the value of the signal strength for which we can exclude it at $CL_s<0.05$.\n", "\n", "So, let's do a scan!\n", "\n", "Ok, before we do so, let me introduce another concept: In addition to the expected $CL_s$ values we usually also show the expected 1 and 2 $\\sigma$ bands to get a feeling in which range we expect to actually observe $CL_s$ values when we do the analysis on real (and therefore fluctuating) data. `pyhf.infer.hypotest` can return us as a second return value the `[-2, -1, 0, 1, 2]` $\\sigma$ bounds on $CL_s$ if we pass `return_expected_set=True`.\n", "\n", "For example:" ] }, { "cell_type": "code", "execution_count": null, "id": "d6edf2ed", "metadata": {}, "outputs": [], "source": [ "pyhf.infer.hypotest(1, [b] + model.config.auxdata, model, return_expected_set=True, test_stat=\"qtilde\")" ] }, { "cell_type": "markdown", "id": "8b4e2d54", "metadata": {}, "source": [ "So the first return value is the observed $CL_s$ value (which is in our case the same as the expected one, since we plugged the exact background expectation into our model) and the second return value is a list of 5 $CL_s$ values for the `[-2, -1, 0, 1, 2]` $\\sigma$ bounds. So in our case the third value in that list is the same as our \"observed\" $CL_s$ value.\n", "\n", "Now for the actual scan - let's see what the lowest possible value (and the 1 and 2 $\\sigma$ bands) for $\\mu$ that is excluded with $CL_s<0.05$. Let's scan with 31 points between 0 and 3:" ] }, { "cell_type": "code", "execution_count": null, "id": "1659aa08", "metadata": {}, "outputs": [], "source": [ "mu_scan = np.linspace(0, 3, 31)\n", "results = [\n", " pyhf.infer.hypotest(mu, [b] + model.config.auxdata, model, return_expected_set=True, test_stat=\"qtilde\")\n", " for mu in mu_scan\n", "]\n", "# for this example we only need the expected band (second return value)\n", "# let's also convert this to a numpy array, such that we can slice it column-wise\n", "results = np.array([r[1] for r in results])" ] }, { "cell_type": "markdown", "id": "def14c8f", "metadata": {}, "source": [ "This is often visualized with interpolated lines and a yellow and green band for the 1 and 2 sigma bounds (\"brazil plot\"):" ] }, { "cell_type": "code", "execution_count": null, "id": "9ae235c6", "metadata": {}, "outputs": [], "source": [ "def plot_scan(scan_parameters, results, exclusion_level=0.05, ax=None):\n", " ax = ax or plt.gca()\n", " ax.axhline(exclusion_level, linestyle=\"--\", color=\"red\")\n", " ax.plot(scan_parameters, results[:, 2], \"--\", color=\"black\", label=\"Expected\")\n", " ax.fill_between(\n", " scan_parameters, results[:, 1], results[:, 3], alpha=0.5, color=\"green\", label=\"Expected +/- 1 σ\"\n", " )\n", " ax.fill_between(\n", " scan_parameters, results[:, 0], results[:, 4], alpha=0.5, color=\"yellow\", label=\"Expected +/- 2 σ\"\n", " )\n", " return ax" ] }, { "cell_type": "code", "execution_count": null, "id": "888b1144", "metadata": {}, "outputs": [], "source": [ "ax = plot_scan(mu_scan, results)\n", "ax.set_xlabel(\"Signal stregth $\\mu$\")\n", "ax.set_ylabel(\"$CL_s$\")\n", "ax.legend()" ] }, { "cell_type": "markdown", "id": "4c1fc03f", "metadata": {}, "source": [ "By looking where the red line crosses the expected line and the error bands we can conclude that the minimum signal strength we expect to exclude in case of no excess events is slightly below 1. We would expect that limit to fluctuate between $\\approx$ 0.6 and 1.4 at $1\\sigma$ level and 0.5 to 2 at $2\\sigma$ level.\n", "\n", "`pyhf` also has a convenience function to run the scan and interpolate this for us, so we don't need to read it from the plot with a ruler ;)" ] }, { "cell_type": "code", "execution_count": null, "id": "11059940", "metadata": {}, "outputs": [], "source": [ "pyhf.infer.intervals.upperlimit([b] + model.config.auxdata, model, mu_scan)" ] }, { "cell_type": "markdown", "id": "8b727fa2", "metadata": {}, "source": [ "
\n", " Exercise 7d [medium]: According to this method, what would be the upper limit on the number of excess events in case of exactly 0 expected background events (no uncertainty)? Follow up question: Is that correct? Or in other words: How well does the asymptotic limit do here? Hint: To create a model with negligible background you have to set the background expectation to a small value (e.g. 1e-10) since the $\\lambda$ parameter of a poisson distribution has to be strictly greater than 0.\n", "
" ] }, { "cell_type": "markdown", "id": "51b3c87c", "metadata": {}, "source": [ "## Run an upper limit scan on signal parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "d69d148e", "metadata": {}, "outputs": [], "source": [ "import json\n", "import mplhep as hep" ] }, { "cell_type": "markdown", "id": "b5cffbe2", "metadata": {}, "source": [ "Now, often we are not only interested in one particular signal, but we might have a certain class of signal models in mind and ask ourselves which parameters of that model are excluded.\n", "\n", "What we can do is run an upper limit scan for each parameter and look for which parameters the excluded signal strength is below 1!\n", "\n", "To get something that looks realistic, i prepared a simplified version of a fit to 9 bins using the [published data from the ATLAS SUSY 1L Wh analysis](https://doi.org/10.17182/hepdata.90607.v3). If you are interested how the procedure of simplifiying the model work, have a look at [dump_signal_grid.ipynb](dump_signal_grid.ipynb).\n", "\n", "Lets load it! We have the following background expectations for each bin:" ] }, { "cell_type": "code", "execution_count": null, "id": "8a8b1503", "metadata": {}, "outputs": [], "source": [ "with open(\"example_background.json\") as f:\n", " b_9bins = json.load(f)" ] }, { "cell_type": "code", "execution_count": null, "id": "179f9fad", "metadata": {}, "outputs": [], "source": [ "b_9bins" ] }, { "attachments": { "image.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAEWCAYAAACpERYdAAAABHNCSVQICAgIfAhkiAAAIABJREFUeF7tXQn8VdP2Xw0q/oYQSRQRqUf6G15IGgz9VfJMRWVqoPfMj1eGyFDmvJIGISHDe0qKl1Dml0opQzLUyxAZUvIkif1f333a57fvuefee+58zrlrfT63fnefPX73ueuss/YaqikmEhIEBAFBQBCIFALVIzVbmawgIAgIAoKARkCYt9wIgoAgIAhEEAFh3hHcNJmyICAICALCvOUeEAQEAUEgggjUjOCcZcqCgCBQQQisXr2aLr74YqpevTo1adKEhgwZUkGrT71UkbxTYyNXBAFBoMwIfP/999SmTRvq3r07PfTQQ7R8+XL65ZdfyjyrcAwvzDsc+yCzEAQEAR8EzjvvPDrkkEOoa9eu+mrNmjW1BC7EWAgIgoAgIAiEEYFnnnmGpkyZQkuXLnWnt8UWWxA+QmIqKPeAICAIhBAB+A4OHjyYOnXqRE2bNtUzXLNmDW3cuDGEsy3PlOT9ozy4y6iCgCCQBoGXXnqJFi1aRD179nRrPfDAA9SnT580rSrrkjDvytpvWa0gEAkEHn30UapVqxadcMIJer5r166ll19+WR9eCjkICPOWO0EQEARChQBUI5MnT6a2bdvS1ltvrec2bNgwbS4oVIWAMG+5GwQBQSBUCMyYMUNL2tB3g+bOnUu//fYbHX300aGaZ7knI9Ym5d4BGV8QEAQSEHjyySf19w4dOtC3335L1113HU2dOlVQ8iAgzFtuCUFAEAgVAm+++SZtt912tN9+++kDy+HDh1OdOnVCNccwTKaaxPMOwzbIHAQBQQAIrFu3jurWrUvt27enRo0aUb9+/ejwww8XcHwQEJ23DyhSJAgIAuVB4MsvvyTYeC9cuFC7xAvjTr0PwrxTYyNXBAFBoMQIwBEHNHDgQPfAssRTiMxwwrwjs1UyUUEg3gjAomTcuHF6kQcccEC8F1uA1QnzLgCI0oUgIAjkhwBsu8866yzq2LEjVatWjX7++ef8OqyA1sK8K2CTZYmCQNgRWLx4MV1wwQXUu3dv+sMf/kArVqxImPKPP/5IAwYMkNgmFirCvMN+V8v8BIEKQABhX1u3bq1XCk9KuMIbQlRBSOV9+/bVLvNCDgJiKih3giAgCIQOAWTL+eijj6h27dpUv359uuqqq2jbbbcN3TzLOSFh3uVEX8YWBAQBQSBHBERtkiNw0kwQEAQEgXIiIMy7nOjL2IKAICAI5IiAMO8cgZNmgoAgIAiUEwFh3uVEX8YWBAQBQSBHBIR55wicNBMEBAFBoJwICPMuJ/oytiAgCAgCOSIgzDtH4KSZICAICALlRECYdznRl7EFAUFAEMgRAWHeOQInzQQBQUAQKCcCwrzLib6MLQjEEIGxY8fSaaedRsuXL4/h6sKzJMlhGZ69kJkIApFH4IknntDR/0DIQzl+/PjIrymsCxDJO6w7I/MSBCKGAFKYIeekIWR/FyoeAsK8i4et9CwIVBQCyPSOuNug4447jk4//fSKWn+pFyvMu9SIy3iCQAwRGDNmjBuDe/vtt6eHH344hqsM15IkJGy49kNmIwhEDoFly5bp7DcbNmzQc586dSp169YtcuuI2oRF8o7ajsl8BYEQIfD7779T9+7dXcbdq1cvYdwl2h9h3iUCWoYRBOKIwC233EILFizQS2vQoAHdc889cVxmKNckapNQbotMShAIPwLvvPMOHXTQQbRp0yad8X327NnUrl278E88JjMUyTsmGynLEARKicDGjRu1ugSMG4SkwcK4S7kDRMK8S4u3jCYIxAKBa665hpDVHdSsWTO69dZbY7GuKC1C1CZR2i2ZqyAQAgT+/e9/U5s2bUgpRTVr1tQ67wMOOCAEM6usKYjkXVn7LasVBPJC4KefftLON2DcoCFDhgjjzgvR3BuL5J07dtJSEKg4BOD+ft999+l147By3rx5VL26yIDluBGEeZcDdRlTEIggAjNnzqROnTrpmf/P//wPLV68mPbaa68IriQeU5ZHZjz2UVYhCBQVgTVr1lDv3r3dMYYPHy6Mu6iIZ+5cmHdmjKSGIFDxCJxzzjn07bffahwQdKp///4Vj0m5ARC1Sbl3QMYXBEKOwCOPPOJK3Qg69eGHH9JOO+0U8lnHf3rCvOO/x7JCQSBnBL766ittx71u3TrdhwSdyhnKgjcseCadFStWEAz4v/vuO60Tu+uuu/RpNGIgvPfee1S7dm0d83fEiBG0++67F3xB0qEgIAgUDoEzzjjDZdxnnnmmBJ0qHLR591Rw5n3hhRfqWL54Ujdu3FgHqwFDP/fcczVTB1111VX0xz/+UXtobbvttnkvQjoQBASBwiNw9913uzG6GzVqRPguFB4ECnpgOX/+fP2KVbduXdeI/4477qBBgwbR4Ycf7q76xBNPJLyOQfoWEgQEgfAhgBjdf/vb3/TEEHTqsccei6SgBUHxiCOOoH333VfzpVdeecUXbAiYDRs2JOTgjAoVlHlPmTJFZ40GvfXWW/p/JCPde++9E/CA6gQ0d+7chHL5IggIAuVHwBuj+9JLL00Qvso/w+AzuOCCC+jvf/+7tpSBNqBFixa+jaHeRQ5OuPtHhQo6UzDuVq1a6bW/+uqr+n9j1G8Dgqcc6JtvvrGL5W9BQBAIAQI33XSTG6Mbb9I333xzCGaV2xR23XVX/caAczZk+6lXr15SR7/99puWuCGZ/9///V/S9bAWFFTyNowbi33ttdeoTp061Lp166S1Q70CwmuKkCAgCIQHAcTovvHGG/WEIIWCqdWqVSs8E8xhJm+88YYOXZsqZC0saL7++ms6//zzaauttsphhPI0KSjzNkvA6wluAuiajIrEXt6sWbP01w4dOpRn1TKqICAIJCHgjdENJh6HaIEvv/yyXqsf80aALQTX2nnnnWngwIFJmIS5oKBqE7NQhIzEq0j79u2T1g69EoLZ7LjjjgQzJCFBQBAIBwJgXiZGNwwMzIFlOGaX+yzAvHHo2rZt26ROxo4dq02YH330Ua02iRIVRfI2+m4/5j169GjCgQjsvsHAhQQBQaD8CIDBGesvBJ2CdUkcogX+97//1cYTfvrujz76iK644gqCeTPC3EaNisK8oe8GQedt06JFi+i2226jY445hvr06RM1rGS+gkAsEYCaE2/BJkY3rDNg1x0HSqXvhrFE586dtUHFnXfeGcmlFpx5//LLL4QDSTzpYNRvbogvvviCevXqpWMkPPvss/o1RkgQKBQCeOWHRy9efevXr0877LADIQ7HfvvtR0uWLHGHGTlyJDVt2lTXg54Tb3+wRmjSpAmZsxhT2Vgg4Drqok+MsXr16kJNOxT9XHTRRdrvAoSgU3379g3FvAoxCT9992effUYdO3ako446Sh/IbrHFFoUYqvR9MHMtKLHKBCk21LXXXquee+45ddJJJylm2OrYY49VDz30UEHHks4EAS8C/IPU99/BBx+sOOuL97L7nd/8dD0O0aC+//77lPX4EE8df/zxipm2YqEkZb2oXmBLC40DPvywUyyRRnUpvvNmazfFgqJiO299nQVH1bx5c8XBtnzrR6mw4AeWRt992GGH6ac4PkKCQKkQgKQNLzqYoaYz+4KkDTISeqr5QSqD+u9f//oX7bPPPqmqRbIcjisI9WpowoQJsYoWaPTdcMyB9RvO2xANEW9Yu+yySyT3zJ50wZm30XcfeuihkQdHFhA9BMwhOAKjpSOWsPTlTCoQCCNHHnlk7Bg31g4VJpIsgM4+++zYBZ2CuqxGjRo668/ChQtp/PjxWpUWFyoo84Z5IMwEIaFAohESBEqNQBDmDYnswQcf1FNLx+RhFTVmzBhidV+pl1H08SBlI60ZCIeTo0aNKvqYpR4AAuSGDRtKPWzJxivogSVy2sEN9ZBDDinZAmQgQcBGwDDvdBI1zFThTQfCj3v9+vW+IEJSQxjUyB5o+a6KCAd2MI8DmaBTkE6FooVAQZn39OnT9epxmi8kCJQDAfPGx4eQ2p/AS8uXLydI3rA2MOTH6NeuXUtz5syJVKwL71r9vkNdBJtmPszVly+//PLIBp3yW19FlRXidJX1SYoPitQ222yjT6z5B6T4kEDNmDGjEN1LH4JAYARYbedaTzBTTmrHul1tecDSp1sP96+X2HlDsbehtzjy39mm2V03B51SbNob+TVV6gIKovNGQCrblrainn6y2FAhYHvtQp9tn73AygD3KiLL2aoSr+T9wQcfaO9CxICOE8H1/corr9RLQrCpOASditP+ZLuWgjDvbAeV+oJAsRCwmbfNlHGYfu+999KkSZP00DAjhAkgdN5e5j1s2LDYHeAhql737t0JwadAcQk6Vaz7KAr9CvOOwi7JHAMjAFMwHMLxq3SCJQkCECEkgx1sH4x+5cqVCcwb5zbwvNtuu+0CjxmFitdff722dQYh6BRieuRKOA+AcUIQQux+E78/U32Ez0DfNl188cWEzFtCyQgI807GREoijADUHXB9h/2ykahxeLlgwQL6y1/+krAyqE/AvI25IKRS5F99/PHHI4xA8tQRMQ8JFkDABw82v3DMxpU8uYfylWAfhXn74y/M2x8XKY0wApCo8aM3TBlqkKuvvjppRV6zQgRkQtqsOETTsxcL6dUQLHBgRRMVEsadeqeEeafGRq5EFAEw5U8++UQz7/fff5+23HJLHVDKSzbzXrVqFSFEaFxiWNtrhb47HUFFdOCBB6ar4l7bY489CJ8ghD6Dxsj2S5QQZIxKriPMu5J3P5u1c7RI5ojEGVyzaVWWuoYpg3nfeuutOqaFH9n1brjhBho8eLBftUiXQYds9MhdunQh44sR6UXJ5DUCBXXSEUxjjAAOqP72t0gs0DDladOm0dFHH01bb72177xNPQTrR8jXxo0b+9aLcuH999/vTl+k2yjvZPLcI8W8ERMchy9CZUAAgZx8PBbLMJOMQxrb7j333FMHX0pFJpM41CpxVJcgauDQoUP18qHHP+uss1JBIeURRCAyzBvZnQ866CDq2bOnzsYjJAikQgASNcwFkXghXdIPI3njfkoXPjbVOMUsR6YXjolPf/zjH3X4UqwDAaTatGlDw4cPJ2S/SUdQlfTo0cONGogDW/OwStdOrkUIgai4lv7www+KvcK0ay/+f/PNN6My9XjME3h36hSJtUycOFGdd955GefKpnGKU/JlrFfqChw0ynVhx/3u92EHI8XBtdTHH3+cND2sv0GDBm47ZvhJdaQg+ghUwxKi8qxBGE8TPB6vxO+++66O1StUAgTmziUaMoRoxowSDFaZQyxbtkzbNCObuU277bYb4WMsaLzoQDKHNQ0ckNAH1IuGkBRlypQpsUg+4F13xX8P8/OHY0wojqes2GbXnSa7+LoSBevwwjz9eM0tQpJ3FIFH+jGkZGOG5H7Y5lx9+umnCcthL0TFuWB9pXG7Lf5mr8ooQiFzDohAqCVvWADg0AXBhJAJAwTnC+i+//Of/+jvjz32mNbtCeWAALKobLMNschW1fiyy4g4UH8SwVYYge39LDdq1yZ6+22iBg2SmklBMATatm1LJgvV/vvvT//85z/TBsaCGeS4ceO06d/nn39OX375Je26667aYuaEE07Q8cqD2lgHm6HUCh0CAZl8Waq1bNnSlTDYS8ydA/TdfHquryEELcdoLsv8IjEoJGYOz8txepUaMECpSZMUv8oo9corSnXpojg+arBliOQdDKccas2bN8+9z1kNqMPVCgkCmRAItbUJJAuErgSNGDGCXnzxRf03TuBNYB1I4ggu7xd4P3RPylJPCHpquIWz1QLdfjvR/PnE5jrEifyI/vxnonvuITZBKPWsZDwPAvfdd59bMoTPFVh9IhgJApkRyMTdS3n9jTfeUKeeeqp66qmn3GFvvvlmVyrhgxn11Vdf6WsIIs9M3L2GekIeBNiaQm3cmFi4apVSSD7gLc8EnkjemRDK+TqSIvAvVX/YxC/nfqRhZSGACGOhIQ5V6d7EeJUEcRxmxV5ybjn+NgR1CV4zcdOL+WCRt1GYd9EA5tgi+h4+4IADijaGdBw/BEKlNrETFyNwPFQi8AxDmE7EaQZBdYLYzCCYC5qs1wjnaefm0xWEBIGQI4BkEOzDoGeJ+1lIEAiKQKiYNzzdoM8GwZpkwIAB+m/YsSKTtyEkTTW2sJyTUGcIMW28MZvdRvKHIBBCBOxEBUGj9YVwGTKlMiAQKuaNw0mY/hnHG+TYg2MO6OSTTyb2mtN/I/M1pGyT0oltwV2phb3LYhdMP/B9gch/iN0M08m99ya66CKiH39M3xxpsfgw2CWYCuItx/vp1InozTeTy1GPH6701Vfpx5GrvgiwetAtb9KkiW8dKRQEfBEIoyZowoQJro4bOm1jCvjf//5X/eEPf3CvifmgZ/fuukupp592Cq+8Eq6zSp10UvothqPT9Onp68jVoiHAPgvu/cyxx4s2jnQcPwRCdWBpw2t7UsKqBAeXIHaJdw8p+WmkXnjhBbfZdddd5/4Q7Dbx27YUK9q0qerCunWObTcYOHuq+tLVVys1erTvJSkMhgA7yyh4PcL3gJ1lgjXaXIvPa9z79cgjj8yqrVQWBELLvDnvoIJpIBg0PmDMhuAyb8pt80EweNt80G5TkVv917860vegQcnLHz9eqYEDk8ulJCMCeAPklGmK44m496G5Hzk6IftCDVBLly5N2w+f4SgOReu2f/vtt9PWl4uCgBeB0DJvTBRStfGkxP92JEHWgbs3vtd8EF6X+DF523gXH/vv773nMO+GDZX6/feq5c6YoThARmJZ7MEozALZ2knVr18/iWkb5m3/j/vyaVZj2fFJZs2apTj5b0J7RAEUEgSyRSDUzBuLGcjSoflBsCmVgkQO8krmd0Hfu5n40NO3jVuhkv5o3txh4LDTBkHC69xZqQ0bKgmFgqzVVsuZexLqjv79+yu2gOKIuZ0CMXXTdtttt1W4V4UEgVwQCD3z9npSQhduyJbM4aQDfbghRBw0PxK7TS4gRboN1E3Qe+MAEzEzWOrjJ1+wJc2erRQ/MBW8MiucbFUd7qsbb7zR1xuSA6kpzoeZUTqHJzGbw1Y4qrL8fBAIdVRB/pFogs03Iq3BRBDE1igE+27QoEGDdJJZEFuisDXbm9rUEHXRxkQftNvoypVC77xD1LIlUdOmxKlYiA3m4Q2SevUrVji5KmvUIOJsLjR7NhHKYpjfMTUIiVc44QHts88+biGHcSD2Bs7YHJEBX3/9dVq5ciWxnpyDLjbQcbfPPPNMnRVHSBDIC4F8OH8p26YyH/RK5nYGFejITfYd2+SwlPMOxVi77eZI3//6V3bTYWsILbWvWJFdu5jV5iBo7lucxMiO2eYWejmvveaY57IqTf35z85nczymQg8VerWJvWCv+SAYN8iOccJPMvXkk0+6zezAVrBEMW0KDWSo+4MZGpjw2LHZTVOYt/r1118V57rUzBuWJPxGlx2GUrtyEEAgOAhKtsnotGmO6vGLLwqOQ6g8LDO9QtielHM53CkStILsGCf43q9fP1q1apW+hqzgfOqv/7bb6IJKIGBUt66z0mefrYQVF3SNH374Ia1evVr32a1bt9AlKi7oYqWz/BC48EJinRhxzrqqfrp2Jda5EWfDzq9vn9aRYt4ITgX3eQSrAt3OMapNjG87xgkCWvXu3VvH+PYGtrLb+OARryJkxOGs48QhA3S2nFmziNavj9cai7watmpyRxD39SKDHaXuOcNXAuFciHPqUosWyatA2QMPEMfzSL6WR0mkmDfWicBVgwcP1ksGcwaTBrMG2ZI5mDoCXYHswFbeNrpCHOmFF4j4YI2uv96JR3LEEQ7jFuk7q93++uuv3foSOCor6KJfGQy3fn3iU2aigQOJXnmFOAQkcZhTIg6Ol0CI+wPiVHRJhD44bR1nhk66lE9B5Jg3Fgt1iYk+CPUI1CQgr2TOdrlaVQKyA1vZbfTFuBEsTBA2lwN2ufSnPzl/sgWEUHAEkLHdECxFhCoEAWSZeuklImQ56t+f6NFHidq1c1SQ06YR3XtvIhDLljnfN2f+SrhoymC9VUgquBa9RB3ikNJ4UjIeOsu8IduZAo49cGcG4X98R31vmxJNu/jD4GDk+OOV+uGHxLFw4l2jhmKfbCeHZZCZyIGl2m+//dz7ZdmyZUFQkzpxQOD555NX8cknSqUKHnbNNY5RwJw5ye3gQAiDgWeeSb6WR0kkJW88vHBIOXr0aPc5Zsf4hmRuDilh521ifMP+GzpzkxcTbYwdeCEfiGXrC+FfETYXUsG22yZOA2Fbu3Qh+vlnovvvr7r2229lm27YB+Z0fPTBBx/oacJnQHTeYd+xAs7vmGOSO8ObV/PmyeUoMdK13+9p3TqnDXwnCkl5MP5QNLU9KREu1pgCItelLZnbbsixNB9kkzZ12mlKvfNO6n2BVFCN/bJ23lkpjobHoRqVOuUUpdav928Tccl78eLFatKkSYrPPtTw4cMVx4dXH330kf9aPaXI4F63bl1X6n700UcDtZNKFYoA4tNAuubYN0l02WXOtSVLki7lUxApO2+/hXpVIXaMb9h784NOf8DITfJib15MxE+JPCFmySuvZF4GXN45V6hq3VqpHj0Q/St1m6FDnZuOmWCU6B//+EdCdElzD5j/27dvr6ZMmZJySWD4NuPuAZyEBIF0CEAoAPN+8MHkWn36KGZABQ8EF3nmDaRsT0r8QO0Y3/C4ND9aRHkzccFtyRzRB+02yehXUMnPPzsS/J/+5DgccHJczozrSOhDhoQaCLbHVieccIK7316m7f3OLuo6DgmYPaL/4W/7TAT1u3TpEuo1y+QKhMCnnyrVu7dSxx6rVIsWSo0YkbljxAqyA4vtu69SCMPspXbtlOrZ01ua9/dYMG+gYKtCEOPbRB/0Zt9BPUO2ZG63yRtV6aDkCEDNsRt7t9kMGg/rofz2MI293B5//HE1hB8+HJMkMHP/q98PseQrkwFLgsDZZzthIKBKbNPGkaJHjkw99Nq1SnXsqDi6WFWdmTMVJyFQHAO4quyll5TaaSelinDYHYnAVPyDzEiw3z7uuONcpx2YBjJz1u2QrJjTTemclzisfPXVV11Tw/PPP5/GjRun69ltMg4oFUKDAIKQIVDUOzCRZGrHJl3IfdoY9v+wy0VgLuNlyte//PJLtqQcqz8cBTBpHW3atNGeuV3hHZeO4Jjx6adVNRC8a4890rXI/tqiRc4aXn7Z6Rv9e9aTfaeeFhhj8WInABkuAasDDyQ66qi8u3Y7wBhPP+18xd/oH+Ow12rBMctl1jhoNAeK8I/ge0DbbH/2WVW56RfONqecQsSmyMxYEkeDeeHddzuB4FAPfgLwcm7WLJdZpW+T+tESvSveQ0rbfBDxvhkJ/fGaD9p5Me020UOgMmcM6drs7RFNmyq1xx6O5AQdpPnw4aOCdOUJw4oUZvfdd59iz1sdWztQHknoNQ88MHkMjIVyP71ntluTbgyMg7Xkm33nqadSrwNjAMd814I5Qm1g74X37xNPDG6+mi2Oudbn3KJ6zs89l9wDsA9B3leeXbzITsSASIJ2jG+8RpsfOaxUDKGOHX3QbhMvdOK5mt13393d1+VexuD3nVOY5URr1qRndvZYYFiony2hDZiD37y9ZXggTZiQ7QhOfZxfePtL9T1X5oq5YY6p+rXLUS/fh1FuSPi3MpZWZ5yReD1EeV9jx7yBdDrzQTsvpm0+aEvmtsmh/85KaVgQQK5I80DuYTODo45SnPjU+eBBjYNX+7qVEzXQWsBUvRJ9t25VY6A/fLfHgBSeLQMHo/SuA04e0J1CUr74YqUaN06sg2vZkHEaMeMAKzBa9GPGAWb2PDCvbAj9eNeBMQyDxnXMw14LGHi2CSqyrR90DXBqY0MGxdmO2P7YaRWyvK+x0XnbyiFvIgY2HyROGKurIObJMZsN8OFOv2DBAu3wA0K5CXRlt7H7lr/zQACxIdaudXSr0Hlut52j+8y1S+5r2hVXUDe4MDNxJAo6BzrhE09M0HG73SP+xHPPVY2GhB5GRw2dcjpaupQ4VKVTY+utHR0m/vcSJ10g1MX/IDhHBdV3IvaFccdHIDFOLmLr6hOGQj0TKwN1W7d2go955+P9nmmd3vpR+N6nj+PGXui5wh0e9+zMmQikRPTII05cEwR7S0e//uo4wiE2yrx56Wrmdy3ogyhq9dKZD8IWnFHTH8T4ts0HbclczAcLsOuQsKAGSPX6DGk2iBoDEiziJeN1H1LgZil4+OZ9xF7yzyzYK7rUixdOe+1VgBvVpwuYC+JeYXWrDjmRKe/rs886JrUwOWzWzLlHfbotVFEsJW/zOLvlllvoyiuv1F8RWfDtt9/W/8PqpG3btm7QKgSwYjMyXc+WzO02pk/5PyACkLAR0XDzG0/GVpDEEcIW/4Mg8cAqwf74dDKHy5CQDLHcVvpcr8giP2uUTZscixLzRrDjjo5kHwQgSPgm8BKkfIQ4tax33C4wBuqZN5Q6dRxrDLTJRBs2EL8GE6EPEN5Y9t3XvxXGMG8dtWsTTZlCdPzx/nXzKcU4e+9NtOWWzng77BC8t1NPJXrrLeRwDN4m25qFegqEsR+vJyUOLA3Z2XfgpANJ3ZAtmdttwrjGUM7J72APuk3oa6H3NLpb6FptKbh2baUOPjiYVIj+mjTRdZm1OJKR0dmm+t+Ahfm1bFk1Dt4KIP3b+lNI+d6Dw2z15Cb5s1kj3himTnX0vhjLjGG/lQCTbAg6ZFufjzeSSy5x+kb/GM/75gPsstXFe/cKB7J4C0L/wA5jet+usj2AxH1h3w/oD3PHOPj4jZHrgW0QjGHLzUYPek58tpIVIewE9qKIxLOKN3nNB3EwacjOiwnzQePYg/gotvmg3SbeaBVodfahGxgLfpSGwLTwg8ePMZMJmfkhow9zAIm+0AcYhn0924MrMBabgdtMw+9vyzopK5TwwPLrz68Ma8yWqWIyXgbn0DBMAAAgAElEQVTu17cpw5qzZaoYA/PyHmKmGgf7lcsYGAf76z1c9hsHdYrJuDdudNRz5h69886stl3HDBLmnR1mfrVtT0qYBNqmgHZeTPxtCHVgashvMtqMUMwH/ZD1KcMPyvzYcFKPMJhGT+2VzPx+lKasVSvFxtfJ1gdg/F6mn+uPOAhDKgSTAEPyWojYa8cY2Ur1XujNWlIxPpTjQZLLw8EeC1ineugVagw8iDFXv7WgDA+RbB/WXrwyfe/bVyl4TCIHJfYq2zeiEjDvWOu8bRWS7UnJUjWxmoQQIhZZeOB9aULDsjROSKkGggfegAED9N92G7vfyP5tPPeMt1u+lh8AAnrqvn2rLCaCgAP9LPTc+MyeTTR9emIr4+2HUswVunSboCffvF9BhvOtgz7ZI5PgMYkxzJiwRMm3b3tA9D11qlOCMc04sI4pFKFfc05gjwHLiUKS2QtgZix2Cj0G5ov+8QFhHDOWU1Kcf3H+BW/Zc85xQijjfOCXX/hQhU9VoIsPQqLzzvR4DH7dG+MEAasMQd8NvTfviZa2oQ83xC7zuhwfeF++xPrUl1mXaD7M9INPotw1MVfoDb32ykYKRDmk5CDSGV6L4X0HnWQqb0OvZA3pE7bQxm7ZDw9IoH4Sl7evXF///caUMkHAIIA3C/wGbDJqwHvuCY5TCSTv2Ou8bbRtT0owY6hTDNnZd2A+aOKCz+YQqttss43LwA0jt/9H+NATeYM5ngbzvTXBN7iUNaEjDqq2ABO3dZZg+mDUYPxelYWXqZrvOOgxTiw4QMwGF4wHBu89JAPDxut0rvrUUuItY0UPAYRH7tcved7//KejOkEY5aAkzDsoUsHr2Z6U3hjfYNqGKbOahflUu7RMOxUjBxMPFXmtJozeEAwS+lhIu15GWaeOY6salOEjjKZh3Jw6TEgQiBQCfMalDxlxUOkleFgigQnu77feqrq6aZO3ZtX3EjDvitF5M6N1yfakRLq0mexBxWoTrfeG/ttko7fb4G8kPfamwmJ1jLYX/8aTXPRA1uEijdYepdDReSdqf4ed9aWXVpUgEtpmm3atP0Y0OXjdQYcJu9QgGa4Rbc7oqfE/+kGCY6OP5nVrL0chQSAKCHCUSZ0+cNKk5PSBZv633eZkkGd+Qc8/T/yDd77jrMSPSqDzrkjmjezxrVq1Yl8Cx92ZY3zToEGD9BZ06NCBk0a/5G7HzjvvTKeffjp17NjRb4vcsq859CPyY7KaxS1jdYruC4y8LISDns2u/3p84zBjDrTwfxBq1IjotNOIOnd2MmibNmg/YoRz2GfKzjor8XuQ/qWOIFBOBBA6mkMK6xCwqQhyN0JHjxrl/KaQI3bo0NQHqMceS7RwIXHMYaJM7vSpxsxQXpHMG5jYnpQmxjeYOGdUcSHry5YTnJklA4SJl8HEhw0blpDYGJ6daRk4mCziQhvrAHSZbyxlMFZ4OBrrhiCrwAm7kai/+47onnuSWxmLAtsKwNQSxp2Ml5RUDgJg1pDQkQh8zhwixAhnIZHq13fif0MaLyBVLPMGhpdccgkLjiw5Mm3LT9J1m7M8w4SQY0QnqUiywR2BsIwUDgkcDDxBhQJGPXGiI6X6ScAwI4PqAWqOTKoXMFKoP9APVCBBgg8Z00AwYzBs/O91eQbjh6kcEhqkI/SF10dRlaRDSa4JAgVFoKKZN2KcQMeNTDuGCsG4TV82A4fkDQauCcwVNqTGftUd3ecPMNS77qqyN7b11EZXbXTNPs3dIltPDUad6YFgGqJvMHEwZ9hxG4KUjn7AsIVpV+EifwkCJUKgopk3MF6yZIl2wOFjYg05QsFm0m8H3RscZl599dWJDkBoDMZtE1JBGXWFkZ7xvy3xIqwogvcEYfj77ecECEIaJoQt/fzzZKk66CKkniAgCIQSgeqhnFUJJzWP4+0axo3DykIxbixha2aceBgYup4ZeQLjBtNGnkVItrAAgTTMMcb1ISMOCW1CjGg/xg0JGLpmSOc4aEUUM0RaA+MGIdaxVx2S2LN8EwQEgQgiUPGSN9QZi6EvZho/fjyfLfDhQoHJVp88xX2zosExTTrssCpX5iB6arSDfrlTJyJOnKyldcOYwdihQ4dFiVGjwPUc/QrzLvCOSneCQPkRqGjmvYIZnsmiAxtuqDiKQezZ6fbNMjI9uMUWRMi2kY7sA8Xvv3d0zjjFtslYfoBZew89hXGnQ1euCQKRRyBAlPTIrzHlAjg+iXtt//33T1kv3wvoGwehSM+mj/z8GLfX8cVrG37uuc6h5ea3BD2nVNK67YiT7+SlvSAgCIQSgYpm3pC8DRkJvFi7hP5h1aJHtCPpGTO9TAOjHqRrWH1ARw7GbR9oQn9uLD9ETZIJTbkuCEQeAWHem7ew2MwbnpqG1jLjrWt01HAKgjMNCAeWUIWAEadiwLC7LmSYUjMp+V8QEAQihYAw783bBcuQYpJ9ELrokEOonckS7h3UxEoAg4YFSSom7m0n3wUBQaCiEKhoU0F4PpaKYPNtaI9UjNueDJg4TAazcW8v1WJkHEFAECg7AhXNvO14I7AIKSaZTD0YYw/8Azvu228ndrusiogNG21I27DdBsGKBNH6vJYkzlX5VxAQBCoYgYpm3rbk7Q3pWuh7wvTPVtoOffYZ0RVXELVv73yg94YlCQ4dwazheGMIdfwcdAo9SelPEBAEIoNARTNvTrbgblQxJW9Oq+bG+64acfPQkK5hOQIPSzBuqEoQiQyWJMZcEHVSxQ2OzK0mExUEBIFCIlDRB5ZQmzRmFcWnHI4VCRWgly7GwSX6NnQiEubCqgTStfnYttuoCCnbK2nfdBPRxx8THXqow9RzCRkLyxYTzMpMyJgqwsJFSBAQBCKDQEV7WGKXhrDEe/1mUz0kXcCnkIQHQr9+/bSDznbsNQnb8qSDUuMhaRgrmDrie2ciw3jxv7Ed92uDQ09k0/E+EOy6OLzlELk6BK2QICAIhB6BimfeYKZIbYbgVIUMB2t2HnHBjeTNSY71wyIQGXXKm28S3XproCa6ElRBhqmDoePBlI3KBW3xdmBUNsFHlpqCgCBQQgQqnnmfcsopNHnyZBdyMPKbWEVRCPXJtGnT6L777tN9p5S6g2w2bL4RdArUpg1xEHJH5WLH1w7SD9rizeL44x3VDdoYr00weOOxCbUOrGBKaEoZZPpSRxAQBKoQqGjmPXbsWBowYIBGo2bNmrRp0yb9dyEY+KxZs9wsPegTuSztA1I9UFCCugOSsGGuONiEdAzmavTmRuXi1Z/7jQHmjP6MhA79OcZAv6Y9rpnkEX59SJkgIAgkI/D66465b716YCrO9cGDiXbZJbluniUVy7xhd42AUdBFgyAh33333W54WLizX3XVVTmlQkNo2enTp7tbM4EZ7dn5urSDScMKxRAYt4llAuYLNQskcdSDtG+cgqqzQdHvv2e+TcCsmzd3+li50qkvWeAz4yY1BAGDAH5/vXo5+St3280pBR9ATP/XXiNq2LCgWFUk80b6s7Zt27q66O7du9Pjjz/OfG8R9ejRgz788EMXZCRn6MMJDYKoUd5k/TQeArbNeFZ67kxbC+kaD4Egh5noC2Fl0cZYtxjpHP9nykuJ9nhAYDww9nQHoqjrR3igwMLFvB2YPtGfyR7k107KBIEoInDAAURduzpZ5e35I/4+fj/ZnF0FWH9FMu9BgwYxjs4hIAJSLViwgBPYbK+l8EM47sgHH3yQBF3r1q11ujSoVEA44KxWrZpuA6YNO3GbaUPHjSQMeUvc3pmAIeLQ09ZRe+vgO572qJdKbw01iVflkomhoy+jajEMHQ8GP0JiZ4yP+aYitIVHqeTATIWQlEcFAfye4KMxaRLRGWckzvqvfyV66CEivNHWqlW4FTETqih64YUXVPXq1RUjqP9nxuuu/7zzztPl+DRv3lwxA3a/m/Ig/3fr1k2xWqb4uD71lFLXXafUWWcpxWPqvydMUGrNmtzGfvttp/3JJyOjZ7DPHnsodeKJSg0ZotTLLyu1YoVSBx4YrK0ZA+1znXNuK5VWgkB+CHzzTWL7xx5z7vmXXkru99ZbnWvLliVfy6MEEmTF0Pfff6922WUXlyGzSsNd+5NPPumWo85XX33F/GSNuuuuuwIxcTD6s5iJ8sFkPPBs3LiKAZ9+uvOAaNkyO6YM5nzKKUrhIWMeZmDSeMDgYWM/IMDAhQSBMCFw//1K7byzUk2aKPW3vznCydq1Sj30kFJnnpk405tucu7nN95IXgHzEH2vz5mTfC2Pks3HoYWT5MPcE5xlVq1apaeItGfXXnut/htluGbo4Ycf5sNh53R49erVrB5mV3Wmvn378plD4qEDHG7gqZmzJYk7asj+QC5MBMUCPfZYlaoG36EzD+pQ9OSTRPhAzWKrXKBLRx9QmQBfOBLBSQjjCgkC5UbgnnuI9aHO4f+SJUSjRhHddpszKzYvpkceSZzhhg3OdxgIpCLmJQWlPBh/pJqOGTPGlaxZv6043og7/6OPPtq9xtne3fJXXnnFLecDS7Vy5cpIrTnvyRqJwUjIUJFccokjgaBzqFkefFCpM85IlKJr1Agmobdrp1TPnkpttVVV/UKoTyDlQ4UDVQ4+mCPmWmjCGPigf/xfjDEMzmasQuCTCgfgVqw1pBozrOXPP588s08+Uer995PLUXLDDc49/Prrydevv965NmNG8rU8SiiPtpFpyoeJir0nXUb8GPRTmwlqEX4a6g8fSKpffvlFX1m3bl2CigVqlYokMHBWCSWoONLpw6FaAYMBI4C6BA/Do44K3h4PCDDcqVOz14OD8UD9kmp+6Pvvf89vG7E2zA99+Y1TiDEwQ6jfzj479RiYQ76MHO2Bh99azD4U4uwG+4K1eMfB2QjWUYgx8tvV/FtPnOjs1YsvJvd12WXOtSVLkq/lUUJ5tI1EUzBjMGXDoKGXNgSmXqtWLX0NzB3fDZ3Mh3amzbnnnhuJtRZtkvhxeXXUXsYFBm+dIfjOxRyIgqEH1Z/jB28fiPp2zIVgQt45pfqe6wEpcAh6GIt6uTJXMLRUc7fLgU2ukjLaeZmp35h16zoP4VS4ZypP9QDyjoVzkCjTRx85e4a3MC/16aMUv+2r33/3Xsnre+xNBS9hPeoImK0xwSwQJn0mkzvM/5AUGMRqFTr//PP13/fff7/Wb4PsNrqgkglmf9BNwyzK2I8bb81czf2M7hymj0G8Q4G/0Z+b/5HE4pxzqnYGNrXGPh2lJqQubM4NYb5wQgpK6AOmYLbpI2KuY/2IJ4N1ABcTxgD94lq2YQag9998v+qpITGH7eDlHQPmm+y9qzEJSjARRYx4ey1w9EIf6A9r8YZegEdvto5mcCrDWIbgd2DPE9ds81SYt0b5zKNZM6IuXYjuuCNxJ4A1zsq8evKg+5WiXqyZ94svvkjHHHOMXjpL2PTqq6/qg0oQGPW4ceP03yxl85kaH6ox2Z6XcJmfN28eOzbyTShUXASMnSxGQT7Rnj1xkuz8+NM4Ja3g6vzToD34wxbj5LIwME58bILXKRjG5jAItPfeRMYTLrFm8je0M8wO82Obf6pTJ7keDq4gEBgPV7hJo24Q+u47py3qwrUa8/Nzq8b8MYaZD+qAsfrNxzsu5sdhIdy2wAgPMq8/AOo991wV80XfGMNvPt4x8B1tceAHQls4qngfMBgDdfCwMMROcgQm6CUr9r57CffHsmXErtDe2tl/R3KUa64h+vprxx67f3+iiy5K38/nnxO98QaxZ59T7/nnnSQqCAGNTFkgrO2005x1bvYRSd9pFlfzkttD3BimfrZZ4M033+zO1msWCBNC0K+//qqYUfMbq6MDZ0eeEK8whlODOsW8Ttuv60Z/DrUM9OeWDv66zXtl9uxs/r7G+0ou34OpYaKKE1RU+RLUO/BR+O03pdq0cfAaOTJ1rzAZ7NgxWV8/e7ZSf/qTY1qIw32Y2X7wQep+8rgSW523bUGCv3/DpjCBqcPaBD92OOnAacfQFVdc4TJudp9nFVVhdVR57FPlNPXq1mGRYluLQFcL/Xb79voHBkbNShJ337CvdWvWVENgpw5G7/0E1bVHlZFV4rwLwbw3bar6jcFiBDjuuqtSdrmpAaOGrl2Veuutsv4uY6k2gVv6pUg+wAS39yVspwm7bWbGdNxxxxHUKaCBAwfSLbfcov+GSgW22rwbCW30RaHSIQBVANzqbb1vptFZ9/wg7x3ON4xNPprswSqBIdwXH1In9wBVga0Dh0rAVh9AjQM9vP1KD1d+6KSDEtaxOdGHbmLGgE4eagro+HGGgHGMCgTXbD1xprEwP+hUDWENGAf9YA0gvNpDLWj3y8HTtFomCEE9gTOgzUHctPoDts5HHFHVGmoihFbGOEZlxOpIuuCCICM4dTAG1CAgzA3rsMdAOdYAG+xPPnHq7bQT0fz5VUm7ndL8/z34YOK4GY76h3lGAuF8BWuDfrucVNZHRxEGty1IGFdlm/hBdYIyfFj37ZoFej0vn3nmmSLMTLrMCgGYymWSkiFdwxxxM8EjFl6zZo/N//xQTvZ8hSWIV8pPJ7VmsqRJtThbFZSuf1wzZpap+kpVDkuNTH2b61A55WLZgf3wMxmFxQpUXN7x8caTLUE95h0DfePtCx/vOKibq7VNprmNGuWsCT4MNl19tVKjR2dqXZLrsVKbcMqxBLNAxCoxhBgmtlmg7aTTuXNn9wdvtynJDsgg6RHADxq25sZeHKaeYKRpfrSIKwOTUC8TZ8k8eSz0ZYcC8GNC+YY8wAMm3RhgQpZzWPIkA5RgjunGwLpwPQ1uGUdB20wPVIyT64MOE8BDNYhfAOaBe6NYxOpV1qsqte22iqU8Z5Tx45UaOLBYI2bdb6zUJrYFCSIAItqfMQtE7G5YkoDs+NowEfzzn/+sy/fZZx8dz7tOkFN73UIozAi8zCoFqE3YU9adZsqkGFA/2BYlsMKAhYPXYiWfBWMMfKAiwQd9p7L2yHUcM4ZZix0J0s9iI5dxMAZUPVDL4GMwMiobr+VKvmOY/YM5I/o26qdc+s2mDfDC2DNnOjHxYerHoTM4nGj6Xn79FfbGRA88QGyulr5uPlezZvchbWBbkEDCth1ubCmMY3e7K+C43YoZtZbQvG1CukyZVg4I8MOahc7GLDS21MHGhASBQAiMGOGoTtjgQR1/vFIbNqRv9uyzTiC23r2VatbMcYJK3yKvq7GQvBFYikO4Ev8w9XOMXd714RUISRZMRng7djcSMhzEuSCNk87w4cPdQ858Hoaxawsbaxy8IRsPnCzWr3cOr9jBSUgQiDUCODzFwemWWxJ98QXRDjsEX+6ppxK99RYcR4K3ybZmXqw/BI1hAug1CzTTgl7bNgu0Y3cjABVjpT9oL+SDwKpVSu2+u1J2kJ7PPlOqYUPFNpY+DaRIEIgRArDl5rAZzCSUWro0u4UhFDIOc4tI1bNl9mGrfxuHaTSmfzAHRDhXEMwCIXEbaXwwJwE13pWob1zm67EHHAeqCtuywjGfoUOJdt+d2E21aj74DvOvq68OxxwLOAvox69n0761ttt4AfuXriKEAPTW0K1zZi1Nzz4bvskX8cFQ9K5tCxJGNsHhxjYZg1mgcdL59ttvFTNsV+q2nXSKPuGoDQDrhH79kmc9bZojjSxcmHwtoiVIooF7CB+O0c5+QOwIJFS5CPTtq9TMmUqZez1b08cSSN6RNRWEWSDrsN0fnB2HG0zdpDpLF7v7wgsvrNyb07tyhAjg8AAuwVQqldnX3LnONWQaiQnBvNCb9o6dfJLtw2OyXllGGgRg6vjAA06F9euV2nJLx2wQv4mgVALmHVm1yV85qacx/YNZINQnIKhJoC6B2gQ0evRoHRkQNHLkSFfFgjZ3eKN/6VoxJgTMQZCkHXckto8kevRRx2SNvUvpzDMTo8wZTzcO6JVEpuybb5IuRbUA3pgr2OyNhQB3Cfjenr0X8VmUjddjVEGQeTsmkDAFNFEqcVgJD0vwkylTwoVQ0AdJmOrZZoHeONwwBWSE9SdV7G6YB8JMsKIIiZYRSOeJJ5SCxHzwwY70DOm6RQulPv00EQ4Elcc1K6CXWwHOGrh2+eWxhBBS+FH8mmzuI/P/2Ry8SEwNY7nlzqJwCO+nJvznP537vXXr4IsvgeQdObWJHVgKPyqkNzMEe17zQ4NKBaoV0M8//6zYAce9Ngqur5VGSKO1cWPiqmFNAr21txy1Xn3VuWGRWNVLnB5OM+8QeZt5p1iI79CDwzbcZuLQh/PBZiG6lz7ChAASsYDh+v0W4GGJRMS45+1gVH5Bq8yahHkn7m4ms0CT6gz6btsskD0o3R8gXOGFAiAASRw36zXXJFc2hzghifGQPMHClkAo8OrDkT5PKCYIIDdtly5K/fBD6gUhPDR+DzArRrRRCD5WVq6khsK8EyGxA0shVjekcBBSncGixEhIsDQxhCBTphxWJiZ2dxLYUpCMAEJiIs6xlx5+2LmRFy3yXontd2/QK0jlQjFBAGqRTMnFwbDxlg8VIxg9Alali61yzDFK7bhjwVOf2YhHRm2SziyQQ7u6DNo2C/SqWMQsMMsfW//+Sh10UHIjqA0aNfKPdZxcO1Yl0Ie/nU9wp1ihIYtJQGDBAqUQfgOu9MhZiaBWMDE87TSl/vGPgoMVCeadziwQDNnPLBCJFJBQwUjdbJ1ScPBC3SHiMFx0kXMz7bWXUjCLXLcu/ZSh27Ptm7/4QqkGDZRCzAZDX36p1G675ZeUNv0s5KogIAgEQCASzNsOLIVM8FCTgLxxuNlT0l3y7bff7jJuu00ATOJRBTrZp5921nLllY6a46ST0q8NOrzp0xPrLFumVK9eTshSHFAixdOMGen7qfCrkM6nTp1a4SjI8ouNQOiZNxiykZ6zMQvk5MG6XUWaBeKusU/CIXHvsIPDwFPl0wtRkPli3/TF7t8cbh7I6blEN15stCu3/1AzbzuwVK5mgffee2/l7q69cqiNcFo+aFAyHiELMp88weiU4GDTa5ly4okn8tnWf6KzCJlpJBAILfOGWaBtQXLyySe7gIKpG7NAxOG2zQL7ckwCI6mLWaB1D773nsO8EREQJ+eGoAKBWsQui8StG95J4kDTz8kH9uHi5BPefYvazELLvO3AUjALNCZ+XrNAmA8ass0C7TZR25Sizbd5c4eBw9sSBKsJ2L2nCzL/0UdKnXOOUqefrhSydJ9wguPYU7RJxqfjpzj9GZJAeJ18JOhVfPa4nCsJJfO2A0vBksQ28bPNAhGH20QLtM0Cq1Wrpjj1VTlxDefYsH+H6gQHmIjL3aEDTn1Tz/X99x3mjrjGIH4b0u7DyO0HRx2hQAjAocerSpGgV4Ggk0ppEAgd84aEbUcLBLM25DULNE46XrPAQX563TQgVMylxYsd5t20qRPnhNVPaQnOCGDyNsHSB6aC7PCUVmJP23HlXYS6xE4AYqRxeG8KCQK5IBA65m0HloLOO5VZIIJTGRo6dKj7atqqVSuObGqFNs0FlTi3AeOF9P2vf6VfJaxVEArz2GOT68F5B33Mn598TUrSIoCDS1sfbocyTttQLgoCHgRqsgQQGnqQM1I/8cQTej7I+o4MN3wgqb/369ePkKsSdN555xEfYOq/+XCIWD+u/0abyZMnE5sJ6u9CPgggPC7y8X32mc9Fq6g6Rwtu1szJW+mtabKDb94P72X5nhoBhJ5FVnvzMblWU7eQK4JACgTC9DizPSKZkbtTQ+RAnr7+wOHGRAv0el7KK2iG3Rw8WKmuXR2pGf/nSnD5ZSsfDi6Taw/SThAQBPJEIFTJGCBdg3r27EnsVan/Rnb3yy+/XP8NKRzSOCRsEEcLdBMyQBLneMu6XMgHgQkTiBBkfuJE4lcTolmznEzwPlXTFqHdK68Qv+4Qcc5QIUFAECgTAnky/6I2h74bkjZDoz92GE47IcPunOF8Xaa4HUWdacg7R/b3Pn2qJgnJGTrrbIPl4PAScY39EjSEHIIoTs8ccCIJhDj5RHEHizvn0B1Y2su1T+dhFmjoM2Yi22yzjWbo3tjdxYUrgr3DwgQxTewg8wg+BeZ96qnBF4T43vwgVRMnBm8jNfNCwDYvNEkgxMknL0hj1TjUzHsHjscBBm3H7vZ6Xl577bWx2pCCLgZRARGe0htkHrrqGjUcaxI2YctICE7FERp11h1DMDPMZGqYsWOpkA4BxEXxOvnAPtw+D0rXXq7FG4FQM++JLOXBrMp20oGLsVGj2LG7471NOawOaiR4T4KB+1G3bo70fccdVVf90jotXaoUQhN8/XViL0iPZqeE8htDyvJGwCSB8Dr5tGvXToJe5Y1utDuohumXSd2e9bBzOfv54YcfrjPDs9qE3n//fWJ9d9b9xL7Bpk049SW65hqi/ff3X+6bbxKDSbTTTkRLlhBtvz1R9+5EDz1EhIzZID4s1pmzOXu6W4bbZeNGouefJ1q+nGirrfz7l9KCIoBM9iy48HkzHzhbxEGviM+CCCaIQhWGQJSePXhl5O3RHzt2d5TWUJK5ImZJkPAAs2crdfjhSiErdo8eil9xqqYHT0p4UUI37vdp0qQkSynVIEjWcdhhh6mm7H26LWdAsYOdlWoOQcZJFfRKzGSDoBevOpGRvPmQklj/px+tvXv3ZgGRJUQhQaBACKxcuZJ9l74gPhinn3/+mVavXk2sqihQ74XvhpM9EBx8Pv30U905qxe1449Q5SAQGeaNLZk5c6b2oOSobPy2Lq/rlXOblmala9eupR133JH+93//l+bPn1+aQfMcBaoUMG34OIifQ55gRqx5pJh3xLCV6UYMgaeffpqgQ4ZTGKfRi9jsZbqVhkCoPCwrDXxZb7gQYNM8PSGoIIQEgbAjIJJ32HdI5lcyBFq2bKnDMXBY4lDru7MFhGMD6eBt9evXp3r16iU1Z7PDpDK/AnYUIs7L6XcpqQxYor5Q8RAQ5l08bKXnCCGAA2B01JIAAAyCSURBVMqd2GwS+u633norQjPPPNXt2QwU+vywUqqHB5g/YhxBlSWUjIDETk3GREoqEAHOvASHNQIjWbp0Kd1yyy3EDjJs3r4lcdIPuvXWW6l169aRRAYHmTjkDyuls5LhmC7CvFNsnDDvFMBIcWUhYPTdHGaYOC8q3X333cT23hoEmKZ27dqVPvnkk0iqU+DEAyaIA1kQonHec889CRu8aNGiwNJ5OmZrdwppH/0GIdT74YcfkqqK1J0EiVsgapPU2MiVCkKAo1dqj128piMpiE133nmntkCZPXs2O5uyt2kECW8R++67L3377bd69nhYpVJXRHB5FTllsTapyG2XRdsIfPPNN5pxN2jQIEkiRT048IC22GKLyAIHvffDDz/szv+MM84gDqMc2fXIxIkKzrwRg6FXr17UqVMn+stf/sJhMDbSJo61cdNNN1GPHj20ZHPSSSfR559/LvgLAqFAwKgBcN+aRB/2xObMmUM1atQgSOdRpuM4Ts0555yjlwA9/vnnnx/l5cjcC+3t34UzjiMSGrvt6hgkN954I+cB6KPeeOMNd6grr7xSsZTDkUp/KPTw0p8gkDUCzMT0vfrcc88ltf3uu+90zHjEPYkDIXVgo0aN9HrxYTf7OCyrItdQUMkbLsXNOGktTHwYTf1kvOOOO2jQoEE6GqAhHELgyT9ixAi3TP4QBMqFACRvqETatGmTNAXO2KSjWJqE10kVIlZgEntXQ0o8JkjiRg8esaVU/HQLyrynTJlCp512mgbV2MoOGDCA9t577wSga9eurb8jxKuQIFBOBFatWqVNAw899FBflckjjzxCLHlrlR+IJXFCmygTBKlLL71ULwEHmbCmEYoeAgVl3mDchxxyiEbh1Vdf1f9D9+0l6MVBOCgSEgTKiYAxEfSzIlm2bBm9/vrrdMwxx1DDhg31NPEmGQdJFeaQeEsGIeDbuHHjyrkNRRsb5xUdO3YkeHzuuuuu+jwuLlRQ5t2qVSsXl9dee43q1Knj69hgIraZH0RcwJR1RA+BdMwb9zDICCCcEJs+/vhjzm+RIsFFhJZfq1YteuKJJ6hmTcfVA+Fl8bCKG8Gdf+TIkQTeBFWteWDFYZ0FZd4GEJggvfPOO3TEEUeQUZHYYM2aNUt/7dChQxwwlDVEGAEwb9yj9pmMWQ7Cw4LgMg+66qqrODkRZyeKCR1wwAE65glow4YNnEipu9bvx4ngIduiRQut+gLFKehYUZj3v//9b+JEwb4ODV9++SXNmzdPx02GramQIFBOBHBQCakTb4leYsspuuCCC7SZK0xckajBfrv01o/idzyQDjroID31BQsW0G233RbFZWScMw6lscc424gLFcU93ui7/fSIo0eP1k93xI4wkk1cwJR1RA+BJcjfmYJgkQE3+TgTJFKoT2DDDul78ODBdPzxxxOk8rgQsnAhPAD4kZ8mIKrrLIrkbXSFXmkG8QvwZMcBENt+RxUzmbcgECsE9tprLzf5BBzqoD6Bc11cyDhhxUllgr0pOPPGoQ4OJPEkh9Ri7L2RHxAnvTBLevbZZ8nYmcblBpF1CAJRRgDqIRPrBKaTkMDjQnFl3gUPTAWpu23btnTttdfqQ6B7771X289+/fXXLvOOy00h6xAE4oSAscaAwQGEK5hJ+h3kRm3NTZo00fFp3n33XR0ad/369TquC8qHDRtGsLyJJBXar5QPd7Tb7YwZMwrdtfQnCAgCRUaAnZJc13m40cOdPspkwnRwog3FXrKKmbheDquH1MEHH6zLokoFV5sYfXecTnUj+VSWSQsCOSDQs2dP6tatm26Jgz7jiZlDV6FoYlQmiBgJb1k46oAQaAzWRJMnT6aFCxeGYq7ZTqKgzBvmgTAT3GeffWiHHXbIdi5SXxAQBEKAwIQJE3RKOND48eO1B2ZUyTBvGEp4DSjMd5guR5EKyrwXL15MP/74o+siH0VAZM6CQKUj4I39DSMDxECJIoF5b7XVVr4OgTAfBEU1TEdBmff06dM1GE2bNo3iPsucBQFBYDMCiP3dv39//Q2xXEwc8CgBhJwBYNAcztc3kYYJjBfVMB0FYd5vv/02NW/enJAuCk9txBKAqSDHR47SXstcBQFBwEJg+PDhxIeWugT5LydNmhQpfAxz9gv1i+iQsD4BRTVMR0E8LOEynM5TLVI7LpMVBAQBjYCJ/Q3mxxYZOnExGB0O/6JAy5cv19M07v/2nKdNm6ZDeODQcs8994zCcpLmWBDJO6lXKRAEBIFYIAA77yuuuEKvBbbRUYpHVK9ePT1v8/Zgb8jYsWP1AWaUE8II847FT0wWIQgUDwFOZeiGUsUB4KhRo4o3WAF7RjhYkPewdeLEidoLfOjQodpRJ6pUcA/LqAIh8xYEBIHUCCDEM9QPiH0CifW9994jxEQJO3Xu3FlL3mPGjNFTfeGFF3QCdDyAECkyyiTMO8q7J3MXBEqIACKBcvJwPSIYOeyjTZzsEk4jq6F++uknuuyyywgxlxDNFN8HDhwYi9CwwryzuhWksiBQuQiA+R155JHaEQ90ww03xCqAVdR2Vph31HZM5isIlBEBuMzDLBgSLFKoIYFDnGJ/lxHarIeWA8usIZMGgkDlIgD9Mfw5QHGM/R2lnRXmHaXdkrkKAiFA4LzzziN4YIIQ+3vQoEEhmFXlTUHUJpW357JiQSBvBOAyv++++2ozvDjF/s4bmBJ2IJJ3CcGWoQSBuCCAqIOIPgiC9+Xpp5+u9eBCpUNAmHfpsJaRBIFYIYC432DaIBxkIpWaUOkQELVJ6bCWkQSB2CEAl/lmzZoRUqiBEIzO6MNjt9iQLUiYd8g2RKYjCEQNAbjMt2/fXk8b6pT333/fTeYQtbVEab6iNonSbslcBYEQIoCs84g4CMJBJpI3CBUfAZG8i4+xjCAIxB6BDRs26Bj+y5Yt02vFYebZZ58d+3WXc4HCvMuJvowtCMQIAXhbIvE43OgRCxwx/v3CscZoyWVdiqhNygq/DC4IxAcBBKsygatgNnjXXXfFZ3EhXIkw7xBuikxJEIgqAjfddJPOGYmwsWJ1UtxdFLVJcfGV3gUBQUAQKAoCInkXBVbpVBAQBASB4iIgzLu4+ErvgoAgIAgUBQFh3kWBVToVBAQBQaC4CAjzLi6+0rsgIAgIAkVBQJh3UWCVTgUBQUAQKC4CwryLi6/0LggIAoJAURAQ5l0UWKVTQUAQEASKi4Aw7+LiK70LAoKAIFAUBIR5FwVW6VQQEAQEgeIiIMy7uPhK74KAIJAFAnPmzKGOHTtSy5Ytadddd6VevXpl0bqyqgrzrqz9ltUKAqFG4MADD6SRI0dSq1atdHYeZOkR8kdAmLc/LlIqCAgCZUBgyy23pBYtWlD16g5rOuqoo8owi2gMKYGporFPMktBoKIQaNKkiZa8165dS7Vr166otQddrEjeQZGSeoKAIFASBJCJ/j//+Y8OLSuMOzXkwrxTYyNXBAFBoAwIIKExSFQm6cEX5p0eH7kqCAgCJUZAmHcwwGsGqya1BAFBQBAoDQJg3rVq1dKmgshKv379elq3bh1BDz5s2DB9TYhIDizlLhAEBIHQIAB9d+PGjWmnnXaitm3barNBMPHffvuNWrdura89+eSToZlvOSciapNyoi9jCwKCQAICRmXSoEEDeuSRRzTjBtWoUYO6dOlCkydPpoULFwpqjIAwb7kNBAFBIDQIGOZ922236STGNpnv8+bNC818yzkRYd7lRF/GFgQEgQQEwLy32mor6tChQxIyMB8EffPNN0nXKrFAmHcl7rqsWRAIIQKff/65a9+9xRZbJM1w7ty5uqxhw4ZJ1yqxQJh3Je66rFkQCCEChjm3adMmaXbfffcdvfvuu7rcTypPalABBcK8K2CTZYmCQBQQWL58uZ7mQQcdlDTdadOmaYsTHFruueeeSdcrsUCYdyXuuqxZEAghAvXq1dOzatSoUdLsxo4dqw8wR4wYkXStUgvESadSd17WLQiEDAGEgwWtWbMmYWYTJ06k+fPn05133qkddYQcBMRJR+4EQUAQCA0CnTt31pL3mDFj9JxeeOEFOumkk2jUqFF01llnhWaeYZiIMO8w7ILMQRAQBDQCP/30E1122WX0yy+/0O+//66/Dxw4kA499FBByIOAMG+5JQQBQUAQiCACcmAZwU2TKQsCgoAgIMxb7gFBQBAQBCKIgDDvCG6aTFkQEAQEAWHecg8IAoKAIBBBBIR5R3DTZMqCgCAgCAjzlntAEBAEBIEIIiDMO4KbJlMWBAQBQeD/ATz37ctwO23yAAAAAElFTkSuQmCC" } }, "cell_type": "markdown", "id": "803b0cae", "metadata": {}, "source": [ "And we have a long list of signal models for a SUSY process with a chargino/neutralino pair ($\\tilde{\\chi}_1^{\\pm}/\\tilde{\\chi}_2^{0}$) that each decay into the lightest neutralino $\\tilde{\\chi}_1^{0}$, while emitting a W Boson and a Higgs Boson.\n", "\n", "![image.png](attachment:image.png)\n", "\n", "The model parameters are the $\\tilde{\\chi}_1^{\\pm}/\\tilde{\\chi}_2^{0}$ mass (assumed to be the same) and the $\\tilde{\\chi}_1^{0}$ mass. We call them `x` and `y` in the following." ] }, { "cell_type": "code", "execution_count": null, "id": "bc16a5ee", "metadata": {}, "outputs": [], "source": [ "with open(\"example_signals.json\") as f:\n", " signals_9bins = json.load(f)" ] }, { "cell_type": "code", "execution_count": null, "id": "4cf37b69", "metadata": {}, "outputs": [], "source": [ "signals_9bins" ] }, { "cell_type": "markdown", "id": "f418bad5", "metadata": {}, "source": [ "Let's plot a few random signal models against the background expectation for all 9 bins:" ] }, { "cell_type": "code", "execution_count": null, "id": "3d96799c", "metadata": {}, "outputs": [], "source": [ "hep.histplot(b_9bins, label=\"background\", color=\"black\")\n", "for i in np.random.permutation(len(signals_9bins))[:3]:\n", " signal = signals_9bins[i]\n", " hep.histplot(signal[\"data\"], label=f\"i={i} x={signal['x']} y={signal['y']}\", linestyle=\"--\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "c19761fa", "metadata": {}, "source": [ "So by running the previous cell a few times you can see the distribution of the different signals is quite different, which makes a multi-bin search very powerful!\n", "\n", "A little background information: Originally those 9 bins originated from 3 different, statistically independent signal regions. The first 3 bins are optimized for signals with low masses, the next 3 bins for signals with medium masses and the last 3 bins for signals with high masses. Here we don't care about the x-axis though, since we just want to do a statistical analysis.\n", "\n", "So let's start with a 1D scan of all points with a neutralino mass (`y`) of 150 GeV:" ] }, { "cell_type": "code", "execution_count": null, "id": "a88481c7", "metadata": {}, "outputs": [], "source": [ "# filter out all signals with y=150 and sort by x\n", "signals_150 = sorted([i for i in signals_9bins if i[\"y\"] == 150], key=lambda k: k[\"x\"])" ] }, { "cell_type": "markdown", "id": "12a5dbae", "metadata": {}, "source": [ "We will run an upper limit scan on the signal strength for each of these signal points (so in some sense we are actually doing already a 2D scan)\n", "\n", "We will again use the `pyhf.simplemodels.hepdata_like` function and completely neglect any uncertainties:" ] }, { "cell_type": "code", "execution_count": null, "id": "bf27006d", "metadata": {}, "outputs": [], "source": [ "limits = []\n", "for signal in signals_150:\n", " model_s = pyhf.simplemodels.hepdata_like(\n", " signal_data=list(signal[\"data\"]), bkg_data=list(b_9bins), bkg_uncerts=[0] * 9\n", " )\n", " limit_obs, limit_exp = pyhf.infer.intervals.upperlimit(\n", " list(b_9bins) + model_s.config.auxdata,\n", " model_s,\n", " np.linspace(0, 5, 51)\n", " )\n", " # for now, we are just looking at the expected limit, since we didn't input any real data \n", " limits.append(limit_exp)\n", "limits = np.array(limits)" ] }, { "cell_type": "code", "execution_count": null, "id": "1e337758", "metadata": {}, "outputs": [], "source": [ "ax = plot_scan([i[\"x\"] for i in signals_150], limits, exclusion_level=1)\n", "ax.set_xlabel(r\"$\\tilde{\\chi}_1^{\\pm}/\\tilde{\\chi}_2^{0}$ mass [GeV]\")\n", "ax.set_ylabel(\"Upper limit on Signal strength\")\n", "ax.legend()" ] }, { "cell_type": "markdown", "id": "c4c2fb43", "metadata": {}, "source": [ "So the point at 300 GeV is at the boundary of exclusion and then the expected exclusion reaches up to around 850 GeV for these signals." ] }, { "cell_type": "markdown", "id": "be4c7091", "metadata": {}, "source": [ "## Let's go 2D\n", "\n", "As you have seen this example has 2 parameters, so we can run a 2D scan against both parameters. For this we don't do an upper limit scan for each point, but just calculate $CL_s$ for $\\mu=1$ and look at the contour of $CL_s=0.05$ in the grid of the 2 parameters." ] }, { "cell_type": "code", "execution_count": null, "id": "22016ff0", "metadata": {}, "outputs": [], "source": [ "x = []\n", "y = []\n", "cls = []\n", "for signal in signals_9bins:\n", " model_s = pyhf.simplemodels.hepdata_like(\n", " signal_data=list(signal[\"data\"]), bkg_data=list(b_9bins), bkg_uncerts=[0] * 9\n", " )\n", " cls_obs, cls_exp = pyhf.infer.hypotest(\n", " 1,\n", " list(b_9bins) + model_s.config.auxdata,\n", " model_s,\n", " test_stat=\"qtilde\",\n", " return_expected_set=True\n", " )\n", " x.append(signal[\"x\"])\n", " y.append(signal[\"y\"])\n", " cls.append(cls_exp)\n", "x = np.array(x)\n", "y = np.array(y)\n", "cls = np.array(cls)" ] }, { "cell_type": "markdown", "id": "9fc256c4", "metadata": {}, "source": [ "For better interpolation, we convert the $CL_s$ values to significances since they change much more linearily and are therefore easier to interpolate between:" ] }, { "cell_type": "code", "execution_count": null, "id": "1370d8d4", "metadata": {}, "outputs": [], "source": [ "z = pvalue_to_significance(cls)" ] }, { "cell_type": "code", "execution_count": null, "id": "3d342cdf", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(ncols=2, figsize=(12, 3))\n", "ax[0].plot(sorted(cls[:, 2]))\n", "ax[0].set_title(\"Sorted CLs values\")\n", "ax[1].plot(sorted(z[:, 2]))\n", "ax[1].set_title(\"Sorted Significance values\")" ] }, { "cell_type": "markdown", "id": "21b3ea17", "metadata": {}, "source": [ "The significance level corresponding to a p-value of `0.05` is" ] }, { "cell_type": "code", "execution_count": null, "id": "26c38a64", "metadata": {}, "outputs": [], "source": [ "level = pvalue_to_significance(0.05)\n", "level" ] }, { "cell_type": "markdown", "id": "6ef6149c", "metadata": {}, "source": [ "You'll have the number `1.64` in your head after some time working for new physics searches :)\n", "\n", "So, let's draw the contour in the 2D grid of our signal mass parameters. For such grids we don't overdo it and leave out the 2 sigma band. We use the `tricontour` functions of `matplotlib` to do a triangulation of our points in 3D (`x`, `y`, `z`) space and draw contours along that interpolated hill.\n", "\n", "Reminder: the columns of expected $CL_s$ values are for `[-2, -1, 0, 1, 2]` sigma." ] }, { "cell_type": "code", "execution_count": null, "id": "c2c72c94", "metadata": {}, "outputs": [], "source": [ "opt = dict(levels=[level], colors=[\"black\"])\n", "plt.tricontour(x, y, z[:, 2], linestyles=\"dashed\", **opt)\n", "plt.tricontour(x, y, z[:, 1], linestyles=\"dotted\", **opt)\n", "plt.tricontour(x, y, z[:, 3], linestyles=\"dotted\", **opt)\n", "plt.tricontourf(x, y, z[:, 2], levels=100, cmap=\"Spectral_r\")\n", "plt.scatter(x, y, c=z[:, 2], marker='o', edgecolor=\"black\", cmap=\"Spectral_r\")\n", "plt.colorbar(label=\"Significance of CLs\")\n", "plt.xlabel(r\"$\\tilde{\\chi}_1^{\\pm}/\\tilde{\\chi}_2^{0}$ mass [GeV]\")\n", "plt.ylabel(r\"$\\tilde{\\chi}_1^{0}$ mass [GeV]\")" ] }, { "cell_type": "markdown", "id": "fd7358f2", "metadata": {}, "source": [ "
\n", " Question 7e: Which values of the signal parameters are expected to be excluded?\n", "
\n", "\n", "
\n", " Exercise 7f [medium]: How far would the expected limit go if we had twice the luminosity? Bonus Question [hard]: How far would the expected 3 $\\sigma$ contour, based on the discovery p-value reach?\n", "
" ] } ], "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.9.2" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }