{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# SMIPP 21/22 - Exercise Sheet 7\n", "\n", "## Prof. Dr. K. Reygers, Dr. R. Stamen, Dr. M. Völkl\n", "\n", "## Hand in by: Thursday, December 9th: 12:00\n", "### Submit the file(s) through the Übungsgruppenverwaltung\n", "\n", "\n", "### Names (up to two):\n", "### Points: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7.1 Unbinned maximum-likelihood fit with Gaussian parameter constraint (10 points)\n", "\n", "This problem is based on the [example of an unbinned maximum likelihood fit](https://nbviewer.org/urls/www.physi.uni-heidelberg.de/~reygers/lectures/2021/smipp/ml_fit_example_iminuit.ipynb) from the lecture. Here we consider the case that we have an external Gaussian constraint on the parameter $a$. The likelihood function can then be written as \n", "\n", "$$ L(a,b) = \\frac{1}{\\sqrt{2 \\pi} \\sigma_{a}} \\exp \\left( -\\frac{(a - a_\\mathrm{meas})^2}{2 \\sigma_a^2} \\right) \\prod_{i} f(x_i; a, b)$$\n", "\n", "The external information on parameter $a$ is $a_\\mathrm{meas} = 0.5 \\pm 0.01$. Perform a unbinned maximum likelihood fit for this likelihood function. Does the constraint for $a$ significantly improve the uncertainty of parameter $b$ in this example?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7.2 Kolmogorov-Smirnov test (10 points)\n", "\n", "The code snippet below reads 100 measurements drawn from the same underlying distribution.\n", "\n", "a) Plot a normalized histogram of the data along with a standard normal distribution ($\\mu = 0$, $\\sigma = 1$), see [*scipy.stats.norm*](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html). Plot also the empirical distribution function (EDF) along with the cumulative distribution function of the standard normal distribution. Hint: convince yourself that `plt.plot(np.sort(x), np.linspace(0, 1, len(x), endpoint=False))` plots the EDF.\n", "\n", "b) Perform a Kolmogorov-Smirnov test using [*scipy.stats.kstest*](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html). Can the hypothesis that the data come from a standard normal distribution be rejected at 95% confidence level?" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "x = np.loadtxt(\"https://www.physi.uni-heidelberg.de/~reygers/lectures/2021/smipp/data_ks_test.txt\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7.3 Significance of a peak: the 750 GeV resonance (10 points) \n", "\n", "In this exercise we would like to test a theory predicting a new particle, with mass 750 GeV, which decays into a pair of photons. We have collected a data sample of 1000 photon pairs in the invariant mass range $500$-$1000\\,\\mathrm{GeV}$ \n", "measured in proton-proton collisions at the LHC. Theory predicts that 30 of these events are signal events from the decay of the new particle. Due to energy resolution the signal events follow a Gaussian distribution around $m = 750\\,\\mathrm{GeV}$ with a width of $30\\,\\mathrm{GeV}$. In the considered mass range the background can be described by the probability density function $f_\\mathrm{bck}(m) = \\frac{1}{\\xi} \\exp(-\\frac{m-m_\\min}{\\xi})$ with $\\xi = 100\\,\\mathrm{GeV}$ and $m_\\min = 500\\,\\mathrm{GeV}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first step consists in choosing a test statistic. For this exercise, we will use the logarithm of the likelihood ratio (according to the Neyman–Pearson lemma):\n", "$$ \\ln Q(\\{x\\}) = \\ln \\frac{L(\\{x\\}|H_1)}{L(\\{x\\}|H_0)} $$\n", "where $L$ are likelihood functions, $H_0$ is the background-only hypothesis (null hypothesis) and $H_1$ is the background-plus-signal hypothesis." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step1: We create a pdf for the background-only model and one for the brackground-plus-signal model (both normalised)." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEOCAYAAABfM7oIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xUVfr48c8zk96AQBJK6B2kGqoIKhZQFLtYVnRZWSy7bhe+u+vq76uru+7u116wu7ZFF3sXxQJSgtRQpEgJNQklpGcyz++PO0B6JpBhkszzfr3ua2buueWZ+4J5cs499xxRVYwxxphAcgU7AGOMMc2fJRtjjDEBZ8nGGGNMwFmyMcYYE3CWbIwxxgScJRtjjDEBF9BkIyITRGSDiGwSkZnVlIuIPOwrXyUiQ8uVPSci+0RkTaV9EkXkMxHZ6HttFcjvYIwx5sQFLNmIiBt4DJgI9AOuFpF+lTabCPT0LdOBJ8qVvQBMqObQM4F5qtoTmOf7bIwxphELZM1mOLBJVbeoagnwOjC50jaTgZfUsQhoKSLtAFT1a2B/NcedDLzoe/8icHFAojfGGNNgwgJ47A7AjnKfM4ERfmzTAdhdy3FTVHU3gKruFpHk6jYSkek4tSXCouNPHdSvV/2iN8aYELds2bJsVU1qiGMFMtlINesqj43jzzbHRVVnA7MBItv11PfnfUvbFlENcWhjjAkJIrKtoY4VyGa0TKBjuc+pwK7j2KayvUea2nyv+/wJ5rO1e/zZzBhjTAAEMtksBXqKSFcRiQCmAO9W2uZd4Hpfr7SRwKEjTWS1eBeY6ns/FXinrkAiw1x8krG3ftEbY4xpMAFLNqrqAW4DPgHWAXNUNUNEZojIDN9mHwJbgE3A08AtR/YXkdeA74DeIpIpItN8RfcD54jIRuAc3+daJUSHs2hLDgcLShro2xljjKkPCYUpBvoNHKIF59/DP68YxGWnpgY7HGNMM1JaWkpmZiZFRUXBDuW4RUVFkZqaSnh4eIX1IrJMVdMa4hyB7CDQaMREuGnRIopPMvZYsjHGNKjMzEzi4+Pp0qULItX1eWrcVJWcnBwyMzPp2rVrwM4TMsPVnNsvha83ZlFYUhbsUIwxzUhRURGtW7dukokGQERo3bp1wGtmIZNszuvflqJSL1/9kBXsUIwxzUxTTTRHnIz4QybZDO+aSMuYcD7JsC7QxhhzsoVMsglzuxjfJ4V56/ZSWuYNdjjGGBNSQibZAJzXP4XcIg+LtuQEOxRjjAkpIZVsxvZKIjrczcdrrCnNGNO8PPXUU9xyyy0V1vXv35/169cHKaKKQirZRIW7OatPMp9k7KHM2/yfLzLGhI5Vq1YxZMiQo5+LiorYvn07PXv2DGJUx4RUsgG4YGA7svNKWPyjNaUZY5qP1atXM3To0Aqfe/XqhdvtDmJUx4TEQ53lndk7mehwN++v2s3o7m2CHY4xphm5+70M1u7KbdBj9mufwF8u7F/ndhkZGVx66aVHuzHn5eUxadKkBo3lRIRczSY6ws34vsl8vGYPHuuVZoxpBnbs2EFSUhLbtm1j69atbN26lSlTpjBw4EDy8/OZOnUqN910E6+88krQYgy5mg3ApIHteH/VbhZt2c+Ynla7McY0DH9qIIGwatUq+veveO61a9cyefJk5s6dy+WXX86FF17IVVddxbXXXhuUGEOuZgNwRu9kYiPcfLC6rqlzjDGm8Vu9ejX9+vWrsC4jI4OBAweSmZlJx47OtGHBvH8TkskmKtzN2f1S+GjNHnvA0xjT5FVONvv370dVSUlJITU1lczMTAC83uD93oVkMxrABQPa8c6KXSzcnMO4Xg0yxbYxxgRF5XsxiYmJ7NvnTGJ86aWXctttt/HBBx9w4YUXBiM8IFSSTV7VmaPH9koiPjKMD1btsmRjjGm2YmNjef7554MdRog0o+VnQaXqY1S4m3P6pfDxmj2UeKwpzRhjAik0kk1ZCexYXGX1BQPbkVvkYcGm7CAEZYwxoSM0ko24YPUbVVaP6dmG+Kgw3ltlvdKMMSaQQiPZRLWAjLegrLTC6sgwNxP6t+XTjL0UldoMnsYYEyihkWyiE6FwP2z+okrRJUM6kFfs4fN1e4MQmDHGhIbQSDaR8RDdqtqmtBHdWpOSEMnby3cGITBjjAkNoZFsRKDfxbD+AyjJr1DkdgmTB3dg/oYs9ueXBClAY4xp3kIj2QAMuAJKC2D9h1WKLh7cAY9X+WD17iAEZowxzV/oJJtOoyAhtdqmtL7t4umVEsc71pRmjDEBETrJxuWCAZfB5nmQX3HiNBHh4iEdSN92gO05BUEK0Bhjjp9NC92YDLgCvB5Y+1aVoosGtQfgnRVWuzHGND02LXRjknIKJPWBVXOqFKW2imF410TeXrETVQ1CcMYYc/xsWujGRAQGTYHP74KczdC6e4XiS4Z0YNbc1azZmcuA1BbBidEY03R9NBP2rG7YY7YdABPvr3Mzmxa6sRl4lTN8zYpXqxSdf0o7Itwu3rKOAsaYJqS2aaG3bNnCtGnTuPzyy4MaY2jVbAAS2kP3s2Dl63DmH52OAz4tYsI5q08y767cyazz+xDuDr1cbIw5AX7UQAKhtmmhu3XrxrPPPhv0ZBOav6aDr4HcTNj6dZWiy09NJTuvhPkbsoIQmDHG1F9t00I3FqGZbHpfAJEtqm1KO6N3Em3iIpmTviMIgRljTP3VNi10YxGaySY8ynnmZu27UJRboSjM7eKyoR34cv0+svOKgxSgMcb475VXXuHqq68++rn8tNA5OTnMmDGD5cuXc9999wUrxBBNNgCDrwVPIax9u0rRFWmpeLxqg3MaY5q81q1b8+STT7J582ZmzZoVtDgCmmxEZIKIbBCRTSIys5pyEZGHfeWrRGRoXfuKyGARWSQiK0QkXUSGH1dwHU6FNr2qbUrrkRzP4I4tmZO+w565McaYBhCwZCMibuAxYCLQD7haRPpV2mwi0NO3TAee8GPfvwN3q+pg4E7f5+MJEAZdDdu/c565qeTKtI78sDePVZmHjuvwxhhjjglkzWY4sElVt6hqCfA6MLnSNpOBl9SxCGgpIu3q2FeBBN/7FsDxz+k8aIrzzM3K16oUTRrUjqhwF28ss44CxhhzogKZbDoA5X+pM33r/Nmmtn1/BTwgIjuAfwDVNkKKyHRfM1t6VlYN3ZiPPHOz4lXwVpwWOiEqnAn92/Luil02ZbQxplZNvbn9ZMQfyGQj1ayr/I1q2qa2fW8Gfq2qHYFfA89Wd3JVna2qaaqalpSUVHOUQ6dC7k7Y9HmVoivSOpJb5OHTtTZltDGmelFRUeTk5DTZhKOq5OTkEBUVFdDzBHIEgUygY7nPqVRt8qppm4ha9p0K3O57/wbwzAlF2XsixCbDsheg13kVikZ1a02HltHMWbrj6KjQxhhTXmpqKpmZmdTYgtIEREVFkZqaGtBzBDLZLAV6ikhXYCcwBbim0jbvAreJyOvACOCQqu4Wkaxa9t0FjAPmA2cBG08oSnc4DLkWFjwMubucpjUfl0u4alhH/vXZD2zLyadz69gTOpUxpvkJDw+na9euwQ6j0QtYM5qqeoDbgE+AdcAcVc0QkRkiMsO32YfAFmAT8DRwS237+va5CfiniKwE/orTi+3EDL0etAyWv1Kl6KphHXG7hNeWWEcBY4w5XtJU2xnrIy0tTdPT02vf6KXJkLMFbl8BrorzP0x/KZ1l2w7w3azxRISF7nOwxpjQIiLLVDWtIY5lv5xHDJ0Kh7bD5i+rFF0zohM5+SV8krEnCIEZY0zTZ8nmiD6TIKYNfP9ClaKxPZNIbRXNq4u3n/y4jDGmGbBkc0RYhDP1wIaP4HDFGozLJVw9vBPfbclhS1ZekAI0xpimy5JNeUOngtcDK6p2FLgiLZUwl/DaEqvdGGNMfVmyKa9ND+hyuvPMTaURBZLjozi3fwpvLMu0EQWMMaaeLNlUNvwmOLgdNn5apejaEZ05WFDKx2uso4AxxtSHJZvKel8ACR1gyewqRaO6taZL6xheXrQtCIEZY0zTZcmmMncYpN0Im7+A7IqDE7hcwnUjO5O+7QBrdtrUA8YY4y9LNtUZegO4I2DJ01WKrkjrSHS4mxcXbj3pYRljTFNlyaY6cUnQ/xJn6oHiwxWKWkSHc9mpHXhn5S7255cEKUBjjGlaLNnUZPh0KDkMK1+vUjR1VBdKPF7rBm2MMX6yZFOTDqdC+yFOU1ql8eN6psQzpkcbXl60DU+ZN0gBGmNM02HJpiYiTu0mewP8+HWV4htGd2H3oSKbWM0YY/xgyaY2/S+F6MRqu0Gf2SeZjonRvLBg68mPyxhjmhhLNrUJj3K6Qa//APZvqVDkdglTR3Vhydb91g3aGGPqYMmmLsOnO7N5LnqiSpF1gzbGGP9YsqlLfFsYcAUsfxkK9lcoOtoNesUu9h0uClKAxhjT+Fmy8ceoW6G0AJY9X6Vo2phulHq9vLTQhrAxxpiaWLLxR0p/6HYmLJ4NnooPcnZtE8t5/dry70XbyC/2BClAY4xp3CzZ+Gv0bZC3B9b8t0rRTWO7caiwlDfSdwQhMGOMafws2fir+3hI7gffPVrlIc9TO7fi1M6teObbH+0hT2OMqYYlG3+JOPdu9q6BLfOrFE8f243MA4V8nGFz3RhjTGWWbOpjwBUQmwwLH6lSdHbfFLq2iWX211vQSjUfY4wJdZZs6iMsEkbOgM3zYPfKCkVul/Cz07uyKvMQi3/cX8MBjDEmNFmyqa9hP4PIBPjmX1WKLhuaSuvYCJ76anMQAjPGmMarzmQjImki8msReUBE/p+IXCkiiScjuEYpqgUMvwnWvlNlJs+ocDdTR3fhyw1ZZOyyIWyMMeaIGpONiNwgIt8Ds4BoYAOwDxgDfCYiL4pIp5MTZiMz4manSe3bB6sUTR3dhfjIMB7/0mo3xhhzRFgtZbHAaapaWF2hiAwGegKhN4NYXBIMnQrpz8IZM6Flx6NFLaLD+cmozjzx1WY27cujR3JcEAM1xpjGobZmtMU1JRoAVV2hqvMCEFPTMPoXzms1PdOmjelKZJiLx+dvOslBGWNM41RbsnlaRDb67tP0O2kRNRUtO8LAKfD9i5CXVaGodVwk147ozDsrdrFjf0GQAjTGmMajxmSjqkOASUAZ8KaIrBCRO0Sk80mLrrEb8yvwFMOix6sUTR/bDbcIT1jPNGOMqb03mqpuUNW7VbUfMBVoCXwhIgtOSnSNXZue0G8yLHm6yvQDKQlRXJGWypvpmew5ZNMPGGNCm1/P2YiIC0gGUnA6DmTVvkcIGfcHKDnsjJlWyYxx3SlT5amvrXZjjAlttSYbETldRB4HMoHfA98CvVX14pMRXJOQ0h/6XwKLn4L8nApFHRNjuGRIB15dvJ19uVa7McaErtqes9kB3A+sA4ao6rmq+pyq+v20oohMEJENIrJJRGZWUy4i8rCvfJWIDPVnXxH5ha8sQ0T+7m88ATNuJpTkw3dVe6b94qweeLzK4/OtdmOMCV211WzGqOppqvqIqu4Vkdj6HFhE3MBjwESgH3B1Nb3aJuI8q9MTmA48Ude+InImMBkYqKr9gX/UJ66ASO4Dp1zmTK6Wn12hqHPrWK5MS+XVxdvZebDGnuTGGNOs1dYbbRuAiIwSkbU4NRxEZJCvaa0uw4FNqrpFVUuA13GSRHmTgZfUsQhoKSLt6tj3ZuB+VS32xbnP3y8bUOPuAE8hLHioStFtZ/UE4NEv7LkbY0xo8qeDwIPAeUAOgKquBMb6sV8HoPzUlZm+df5sU9u+vYDTRWSxiHwlIsOqO7mITBeRdBFJz8o6Cf0Zkno5UxAsfQbyKua/Di2jmTK8I2+k72B7jj13Y4wJPX71RlPVyvMdl/mxm1R3KD+3qW3fMKAVMBKn08IcEamyvarOVtU0VU1LSkryI9wGMPYP4CmqtnZz65k9cLuEh+ZtrGZHY4xp3vxJNjtEZDSgIhIhIr/D16RWh0ygY7nPqcAuP7epbd9MYK6v6W0J4AXa+BFP4LXp4YwqsPQZOLSzQlFKQhQ/GdmZt5ZnsjkrL0gBGmNMcPiTbGYAt+I0Y2UCg32f67IU6CkiXUUkApgCvFtpm3eB63290kYCh1R1dx37vg2cBSAivYAIIJvG4sxZoF6Yf1+VohlndCcq3M2Dn1vtxhgTWupMNqqararXqmqKqiar6nWqmuPHfh7gNuATnJrQHFXNEJEZIjLDt9mHwBZgE/A0cEtt+/r2eQ7oJiJrcDoOTNXGNA9zy07OBGsrXoGsDRWK2sRFcsPoLry3chdrdtp8N8aY0CE1/U6LSH+gu6q+6/v8f0ALX/Gjqvr9yQnxxKWlpWl6evrJO2F+Djw0CLqNgymvVCjKLSpl3N+/pF/7BF6eNoJqbjcZY0yjICLLVDWtIY5VW83mfio2T50HfAB8CdzZECdvtmJbw2m3w/r3YceSCkUJUeHcdlZPFmzK4euNjaf1zxhjAqm2ZNNOVReW+5yrqv9V1X/TWG7IN2ajboHYZPjsL1Cp9njdyE50TIzm/o/W4/U2nhZAY4wJlNqSTXz5D6o6stzH5MCE04xExDqDdG5fCBs/rVAUGebmd+f2Zt3uXN5esbOGAxhjTPNRW7LZJSIjKq/09Rqr3IXZVOfUGyCxG3x+F5R5KhRdOLA9Azq04J+f/kBRqT+PLRljTNNVW7K5A3hdRP4iIhf6lruA14A/nJTomjp3OJx9F+xbC8tfqlDkcgmzJvZh58FCXvpuaxCCM8aYk6e2sdGWACMAN3CDb3EBI31lxh99L4LOp8EX90BRxe7Oo3u04YzeSTz6xSb255cEKUBjjAm82qYYmA2cBjygqpf5ljtVde/JC68ZEIHz/urM5Pn1A1WK/+f8vuSXlPGvzzZUs7MxxjQPtTWjPQcMAj4UkXkicoeIDDpJcTUv7QfDkGth0ZOQU3Fem14p8fxkZGdeXbyddbtzgxSgMcYEVm3NaItU9S5VPR24EtgO/FZElovIcyJy5UmLsjk4688QFgmfVX1E6ddn96JFdDh3v5dBYxoMwRhjGoq/oz7nqOprqnq9qg7BmdisZ2BDa2bi28Lpv3Ee9NzyVYWiFjHh/Obc3izasp+P1+wJUoDGGBM4fiWbagxU1XsbNJJQMPJWaNEJPp5VpSv01cM60qdtPPd+uM66Qhtjmp3jTTZ3N2gUoSI8Cs67F/ZlwNKnKxSFuV3cOakfmQcKeeabLUEK0BhjAiOspgIRWVVTEZASmHBCQN8LocfZ8MW90O9iSGh3tGh0jzZM6N+Wx77czMVDOpDaKiaIgRpjTMOprWaTAlwPXFjNUucUA6YGIjDx71BWAp/+qUrxnyb1BeDu99ae7MiMMSZgaks27wNxqrqt0rIVmH9SomuuWneHMb+GNW9W6SyQ2iqG28/uyWdr9/JphnUWMMY0D7V1fZ6mqt/WUHZN4EIKEWN+Ba26wIe/A0/F0QOmjelKr5Q47n5vLQUlnur3N8aYJqS2EQTi6trZn21MDcKjYeIDkP0DfPdoxSK3i3svGcDOg4U8NM+mkDbGNH21NaO9IyL/FJGxIhJ7ZKWIdBORaSLyCTAh8CE2Y73OhT6TnGFsDmyrUDSsSyJXpqXy7Dc/smHP4SAFaIwxDaO2ZrTxwDzg50CGiBwSkRzgZaAtMFVV3zw5YTZjE+4HBN7/dZVJ1mZO7Et8VBh/enu1TbJmjGnSan3ORlU/VNVrVbWLqrZQ1daqOlpV71VVu3vdEFp2dKYh2DwPVv2nQlFibASzJvZl6dYDvL50R1DCM8aYhnC8D3WahjTsZ9BxJHw8E/L2VSi6Ii2V0d1b89cP17HrYGGQAjTGmBNjyaYxcLngokegJB8+qjgvnYhw/6UDKfMq//PWahuo0xjTJFmyaSySesG4P0DGW7D+gwpFnVrH8IcJvZm/IYu53+8MUoDGGHP8/Eo2ItJKRAaKyNAjS6ADC0mn/QpSToEPfguFBysUTR3VhbTOrbj7vQz25RYFKUBjjDk+dSYbEflfYBXwMPBP3/KPAMcVmtzhMPlRyNsLn/6xQpHLJfz98oEUe7z88e011pxmjGlS/KnZXAl0V9UzVPVM33JWoAMLWe2HODWc5S/D+g8rFHVLiuM35/Tis7V7eW/V7iAFaIwx9edPslkDtAx0IKacM2ZB2wHw3i8hP7tC0c9O78bgji3589tr2HPImtOMMU2DP8nmPmC5iHwiIu8eWQIdWEgLi4BLZkPRIXjv9goPe7pdwr+uHESJx8vv3lhpD3saY5oEf5LNi8DfgPs5ds/mn4EMygAp/WD8nc400iterVDULSmOP0/qx7ebsnlh4dbgxGeMMfVQ4+Rp5WSr6sMBj8RUNfJW2PAxfHQHdBkDrTofLbp6eEfmrdvL/R+v57QebejdNj6IgRpjTO38qdksE5H7RGSUdX0+yVwuuOQJ5/3bN4O37GiRiPC3yweSEBXG7a8vp9hTVsNBjDEm+PxJNkOAkcBfsa7PJ1/LTnD+A7BtAXxTsfWyTVwkf7tsIOv3HOafn/4QpACNMaZudTajqeqZJyMQU4tBU2DLlzD/Puh8GnQ57WjR+L4pXDuiE09/s4XTerRhXK+kIAZqjDHVk5oeDhSR39S2o6r+KyARBUBaWpqmp6cHO4wTU3wYnhoHpYUw41uIbX20qLCkjIsfW0B2XjEf3n46KQlRQQzUGNNciMgyVU1riGPV1owW71vSgJuBDr5lBtCvIU5u6iEyHq54Hgqy4Z1bKnSHjo5w89i1QyksLeOXry3HU+YNYqDGGFNVbZOn3a2qdwNtgKGq+ltV/S1wKpDqz8FFZIKIbBCRTSIys5pyEZGHfeWrync88GPf34mIikgbf2JpFtoNgnPvgR8+hkWPVyjqkRzHPRefwuIf99tU0saYRsefDgKdgJJyn0uALnXtJCJu4DFgIk5N6GoRqVwjmgj09C3TgSf82VdEOgLnANv9iL95GT7dmUr6s79A5rIKRZcOTeXKtFQe/XIT32zMClKAxhhTlT/J5t/AEhG5S0T+AizGedCzLsOBTaq6RVVLgNeByZW2mQy8pI5FQEsRaefHvv8H/AEIvcfnRZy5bxLawZyfQF7FpHL3RafQMzmOX72+gr02OrQxppGoM9mo6r3AjcAB4CBwo6re58exOwDl5zLO9K3zZ5sa9xWRi4CdqrqytpOLyHQRSReR9KysZvZXfkwiXPUyFOTAmzdCmedoUXSEm8d9929ueeV7Sjx2/8YYE3w1JhvfD/VDIjIBWKuqD/mW5X4eW6pZV7kmUtM21a4XkRjgj8CddZ1cVWerapqqpiUlNcPuwO0GwaT/g63fwLy7KhT1SI7ngcsHsWzbAe5+LyM48RljTDm11WxGAm8BZwBficiHInK7iPTy89iZQMdyn1OBXX5uU9P67kBXYKWIbPWt/15E2voZU/My+BoY9jNY+Igzw2c5Fwxsx81ndOeVxdt5fUno3doyxjQutfVG86jqfFWdqaojgGnAYeAeEVkuIo/XtK/PUqCniHQVkQhgClB5tOh3get9vdJGAodUdXdN+6rqalVNVtUuqtoFJykNVdU9x/Hdm4fz7oPU4fD2rbBvXYWi353bm7G9krjznQy+334gSAEaY4yf00L75ANvqOqVON2fX6ltY1X1ALcBnwDrgDmqmiEiM0Rkhm+zD4EtwCbgaeCW2vatR6yhIywCrnwJImLhtauhYP/RIrdLeHjKYNq2iOLml5ex77B1GDDGBEeNIwgc3UBkGPAczgOegtNJ4KequqzWHRuRZjGCQF12LIEXJkHqMPjJW04S8lm/J5dLHltIn3bxvHbTSKLC3UEM1BjTVJysEQSOeBa4xdd01Rm4FXi+IU5uGlDH4TD5Udj2LXzwmwojDPRpm8D/XTWI5dsP8vs3V9mEa8aYk86fZHNYVb858kFVv8W5d2Mam4FXwtjfw/J/w3ePViiacEo7Zk7sw3srd/Hg5zZCtDHm5PJn8rQlIvIU8BpOt+SrgPlHhpZR1e8DGJ+przP+B7J/gE//DK17QO+JR4t+PrYbP2bl8/AXm+jcOpbLTvVr1CFjjDlh/iSbwb7Xv1RaPxon+ZzVoBGZE+NywcVPwsHt8OY0uPEDaD8EcCZc+9+LT2HHgQJmzl1FaqtoRnRrXccBjTHmxNXZQaA5CIkOApUd3gPPnAOeQpj2KSR2O1p0qKCUS55YwP78Et6cMYoeyTaltDGmqobsIGDz2TRn2Rvh2XMhqgVM+wzijo2ksD2ngEufWEiEW/jvLaNp1yI6iIEaYxojm8/G+KdNT7hmjlPLefUKKM47WtSpdQwv3DiM3CIP1z+7hIMFJbUcyBhjTkxA57MxjUDHYXDFC7B7Fcy5HspKjxad0qEFs68/lW05BUx7MZ3CkrLgxWmMadYCNp+NaUR6T4ALH4TN82DudPAeSyqju7fhoSmD+X77AW599XtKbZZPY0wABHI+G9OYDL0ezvl/kDEX3v0FeI8llYkD2nHPxafwxfp9/HbOSsrsoU9jTAOrs+uzqt4rIh8Bp/tW3ViPaQZMY3La7VBaCPPvg/BoOP8fzmRswLUjOnO4yMP9H60n3O3igcsH4nJVN9ODMcbUnz/P2Rx5cNMe3mwOxt0BJfmw8GEn4Zzzv0cTzoxx3Sku9fJ/n/9ARJhw78UDLOEYYxqEX8nGNCMiTnNaaaEzD054LJw562jxL8f3oKSsjMe+3Ey428XdF/VHxBKOMebEWLIJRSIw8e9OwvnqfkDhjFkggojwu3N7U+Lx8vQ3PxLhdvHHC/pawjHGnBBLNqHK5YKLHnbef/U38BTD2XcdTTj/c35fSsuUZ779EY9XuXNSP2tSM8YcN0s2oczlhosecea+WfAglJXAeX89mnD+cmE/XCI8t+BHikrLuPeSAbgt4RhjjoMlm1DncsEF/wJ3JCx63KnhnP8PcLkQEf48qS+xkW4e+WIThaVl/POKQYS56zPBqzHGWLIx4NzDmXCfr4bzkJNwLnwI3GGICL89tzdR4W4e+GQDRaVlPHz1ECLDbLZPY4z/7E9U4xCBs+92OgqseBn+cx2UFBwtvvXMHtw5qR+fZOzlppeWkV/sCWKwxpimxpKNOUYEzpzydDgAABtmSURBVJjpNKP98DH8+xIo2H+0+KdjuvL3ywby7cYsrn56Edl5xUEM1hjTlFiyMVUNv8kZvHPX9/D8+XBo59GiK4d1ZPZP0vhh72Eue2IhW7PzgxenMabJsGRjqtf/Yrjuv3Ao05kTZ9/6o0Vn90vh1ZtGkltYymVPLGTljoNBDNQY0xRYsjE16zrWmVbaWwrPngOb5h0tGtqpFW/ePJroCDdTZi9i3rq9QQzUGNPYWbIxtWs3CH42D1p2gleugCVPHy3qnhTH3FtG0z05lp+9lM7TX28hFKYZN8bUnyUbU7eWHeGnH0PPc+DD38GHf4AypzdacnwUc34+iomntOXeD9fxhzdXUeKxOXGMMRVZsjH+iYyHKa/CqNtgyVPw2lVQ6NyriYkI49Grh/LLs3rwxrJMrntmMTnWU80YU44lG+M/lxvOu9d54HPLfJh9BuxZ4xS5hN+c25tHrh7CysyDTH5sARm7DgU1XGNM42HJxtTfqTfADR84o0Y/czaseuNo0YWD2jPn56PwlCmXPr6QN9J3BC9OY0yjYcnGHJ9OI+HnX0P7ITD3Z/DRHVBWCsCgji15/5djOLVzK37/5ipmzV1FUWlZkAM2xgSTJRtz/OJTYOq7MOJmWPwkvDAJDjo1mTZxkfx72ghuPbM7ry3ZweVPLmTH/oI6DmiMaa4s2ZgT4w6HiffDZc/C3gx4cgyse88pcgm/P68Pz1yfxracAiY98i0fr9kT5ICNMcFgycY0jAGXw4yvIbGrM4jn+79x7ungjDjw/i/G0CkxhhkvL2PW3NUUlNhAnsaEEks2puEkdoOffgqjfwnpz8LTZ8HetQB0bh3Lf28ezc/HdeO1Jdu58JFvrbeaMSHEko1pWGERcO7/OuOq5WfB7HHwzb+gzENEmItZE/vy8rQRHC7ycMljC3nmmy14vTbqgDHNXUCTjYhMEJENIrJJRGZWUy4i8rCvfJWIDK1rXxF5QETW+7Z/S0RaBvI7mOPU42y4ZRH0ngjz7obnzoWsDQCM6dmGj381lrG9krjng3Vc88wituc0484DpUWweyVs/BzW/BfWvQ9bF0DubrDhfUyIkECNZSUibuAH4BwgE1gKXK2qa8ttcz7wC+B8YATwkKqOqG1fETkX+EJVPSLyNwBVvaO2WNLS0jQ9Pb3Bv6Px05r/wge/dSZjO+tPMPIWcIehqsxJ38E976+jTJWZE/tw3YjOuFwS7IhPjNcLmUucOYE2fg771oLW0PU7LsUZBqjvRdB9PLht8lzTeIjIMlVNa5BjBTDZjALuUtXzfJ9nAajqfeW2eQqYr6qv+T5vAM4AutS1r2/9JcDlqnptbbFYsmkEDu+F938NGz6AtgNh0oOQeioAuw4WMnPuar7+IYuR3RL5+2WD6NQ6JsgBH4f8bFjxCix7AfZvAVcYdBrlPJOU0h8SOjjD/pSVQEEOZG+EHYth42dQnAstO8OoW2HoVAiPCva3MaZBk00g/4zqAJR/fDwTp/ZS1zYd/NwX4KfAf6o7uYhMB6YDdOrUqT5xm0CIT4Epr0DGW/DxLHhmPKT9FMbfSfuWLXnxxmFHaznnPfg1vzmnFzec1oVwdxO4rViwHxY8BIufAk8hdBoN4+6AXhMgupZW3u5nwYifg6fYqQV99xh89AdY9Dicey/0ucCZPdWYZiCQ/5Or+19SuRpV0zZ17isifwQ8wCvVnVxVZ6tqmqqmJSUl+RGuCTgROOVSuG0pjJgBy56HR4fBytcRVa4a1olPfj2W0d1bc++H67jwkW9Ztu1AsKOuWXEefP0APDTYSTb9LnLuU/30Ixg0pfZEU15YJPSbDNM+hZ+8DWHR8J9r4Y0boLARf39j6iGQySYT6Fjucyqwy89tat1XRKYCk4Br1SZQaXqiEpwHQafPd6YveOvn8PSZ8OM3tG8ZzTNT03jqJ6dyyDcT6Ky5qzlYUBLsqI/xFDu1mIcHwxf3QJfT4OYFcOlsSO57YsfufibM+BbG3wnr34cnTnM6ExjTxAXynk0Yzk3+8cBOnJv816hqRrltLgBu41gHgYdVdXht+4rIBOBfwDhVzfInFrtn04h5vbB6Dsz7f5C7E3pNhHPuhqTe5Bd7ePDzH3huwVZaRodzx8Q+XD40NXgdCLxlsPJ1mH8/HNoOXU53kkLH4YE5387vYe5NcGArnP+A0+xozEnUJDoIwNHeZg8CbuA5Vb1XRGYAqOqTIiLAo8AEoAC4UVXTa9rXt34TEAnk+E6zSFVn1BaHJZsmoLQQFj3hPJNTWgADr4TTfwtterJ2Vy5/ens1328/yIAOLfjzpH4M75p48mJTdYbg+eIeyN4A7QbD2X+BbmcG/p5K0SF4cxps+gyGT4cJ9ztTPRhzEjSZZNNYWLJpQvKznYST/hyUFUP/S+D03+FN6st7q3Zx/0fr2X2oiPMHtGXWxL50TAxgrzVV58b9l3+FPaugTS+n63bfi07ujXtvGXx2J3z3qHM9LpntPDxrTIBZsqknSzZNUF6W8+O69BkoyYOuY2H4dAq7nsvTC7bzxPzNlHmVG8d04ZZxPWgRE95w5/aWwQ+fwNd/h13LoVVXp3fZgCuC+xzMgoecpNPjbLjy3xDRBLuHmybFkk09WbJpwgr2w/cvwtJn4dAOSEiFgVeS3fUi/rpMeGv5TuIiw5gxrjs3jO5CbOQJJIPCA7D8FVgyGw5ug5adYOwfnJ5l7gZMZidi2Yvw/q+g4wi49g3nuR1jAsSSTT1ZsmkGvGVOk1b6c7D5S+eJ/OR+ZLc/k+f3dOXprUkkxMVy65k9uGZEJyLD/LyvkZ8Dmz6HjLmwaR54S50HMUf8HPpMajxJpryMt5z7OB1HwHVvQkRssCMyzZQlm3qyZNPM5GXB2rdhzVxnWBivh7KwaDZJVxYXdmBvVDeGDzqF0YNPITymhXNDXb1OzSU/B3I2QdZ6yEyHfb7OkQkdnPshA6+EdoOC+/38sWYu/HcadD4NrpljTWomICzZ1JMlm2as+DBs/Ra2fIXuXkHZrtWEefLq3i860UkqXU+HruOg/VBwNYHRCspb9Qa8Nd3pgn3NfyA8OtgRmWamqQxXY0zgRcY7I0v3nogAYV4vengXyzPWM2/JSnbty6JFlIuxvZIZ2b87MS1TnHl3YtsEO/ITN/AK8Hrg7Zvh9Wthyqs2pppptKxmY5q1xVtyePTLTXyzMZuEqDCuG9mZ60d1oW2LZvSj/P2/4d3bnLHYrvy3dYs2Dcaa0erJko1ZueMgT8zfzKdr9+ASYdLAdkwb040BqS2CHVrDWPosfPAb6HshXP584+zYYJocSzb1ZMnGHLE9p4AXFm5lTvoO8oo9DO+SyE/HdOGcfm1xN/V5dBY9CR/fAadcBpc+bSMNmBNmyaaeLNmYynKLSpmzdAcvLNxK5oFC2rWI4qphHblqWEfatWjCN9qPPPg56GqY/HjT6/RgGhVLNvVkycbUxFPm5fN1+3hl8Ta+2ZiNS+CsPilcM6Ij43olN83azlcPwJf3wNDrYdJDlnDMcbPeaMY0kDC3iwmntGXCKW3ZnlPA60u3Myc9k8/X7aVDy2iuSEvl0iGpTWvm0HG/d8aV+/oBcEfA+f+wSdhM0FnNxphKSjxePl+3l1cXb2fB5mxUIa1zKy4Z2oFJA9o37DhsgaLqNKctfBhG3grn3WsJx9SbNaPVkyUbc7x2HSzk7RU7mfv9TjbtyyPC7WJ832QuHZrKuF5JRIQ14iYqVWcK7sVPwGm/grPvsoRj6sWSTT1ZsjEnSlVZszOXucszeXfFLnLyS4iPCuOcfimcf0o7Tu/Vxv/x2E4mVadLdPpzzlTc591n93CM3yzZ1JMlG9OQSsu8fLsxm/dX7eaztXvILfIQHxnG+L7JnD+gHWN7JREV3ogSj9cLn/4JFj0GA6+CyY/ZczjGL9ZBwJggCne7OLNPMmf2SabEM4AFm7P5aPVuPl27l7dX7CI2ws3pPZM4q28yZ/ZOJik+MrgBu1zOPZuYRPjif6EoF6543sZSMyeV1WyMaSClZV6+25zDxxl7+GLdPvbkFiECg1JbMr5PMmf1TaZfuwQkmPdNlj4DH/wOOgx1xlKLbxu8WEyjZ81o9WTJxpxsqsra3bnMW7ePeev3sXLHQQDaJkRxes82jOnZhtHd2wSn1rP+A/jvTRDVAq5+DdoPPvkxmCbBkk09WbIxwbbvcBHzN2Tx5fp9LNycw6HCUgD6tI1nTA8n+QzvmkhMxElq2d6zGl67GvKz4cIHndlIjanEkk09WbIxjUmZV8nYdYhvNmazYFM26VsPUFLmJdwtDOnUiuFdEknr0opTO7ciPiqAN/Lz9sEbN8C2Bc7wNuc/YNNMmwos2dSTJRvTmBWWlJG+bT/fbspm0eYc1uzKpcyruAT6tktgWJdEZ+naiuT4Bp4aocwD3/wDvvobtOoKl86G1Ab5bTHNgCWberJkY5qS/GIPy7cfZMnW/aRv3c/32w9QVOoFoHPrGAaltmRQx5YMSm1B//YtiI5ogG7WWxfA3OmQuxOGT4fxf7ZajrFkU1+WbExTVlrmZc3OQyzdup9l2w6wKvMQuw8VAeB2Cb1T4o8mn4GpLemZEke4+zge3CzKdbpGL3kaEtrD+DthwJX2EGgIs2RTT5ZsTHOzL7eIlZmHWLnjICszD7Jyx0FyizwARLhd9EyJo2+7BN8ST792CbSM8XMGzx1L4cPfwe4V0HYAjP8L9DjbhroJQZZs6smSjWnuVJWtOQWsyjzI2l25rN2dy7rdh8nOKz66TfsWUUcTUM+UOHokx9GtTVz1zXBeL2TMhXl3w8HtkHIKnHY79L/ERh8IIZZs6smSjQlVWYeLWbf7SPJxls1Z+ZR5nf/3ItChZTQ9kuPokeQkoCNLy5gI8JTAmjedSdmy1kNcW6eb9JDroE3PIH87E2iWbOrJko0xxxR7ytiaXcCmfXnOkuW8bsnKo9jjPbpd69gIOreOoUvrWDolRjHcs4x+u+bSIvNLRMug/VDoc4GzJPWxZrZmyJJNPVmyMaZuZV5l54FCNmUdZtO+PDbvy2fb/ny25xSwO7eIIz8VSRzkqsiFXBi2hN5lPwBwODqVw+1GIV1OI6H3OGKTu1ryaQYs2dSTJRtjTkxRaRmZBwrYlnNkyWfb/gLys3bQ5/ACxrKCYa71tJR8APbRiq1h3dkX24u8Vn3xJp9CfLvutE+MJyUhiqT4yMY5JYOpwJJNPVmyMSZwyrxK1uFidh7I5/D2VYTtWEhc9ira5G2gXek23DhNc6XqJlPbsE3b8qO2ZV94ewqi2+GNb4+rRXuiW7YjOSGalIQokhMiSY6PJDk+qmGeIzLHxZJNPVmyMSZISosgaz1FO1dRsPsHPNmbCTv4I3H524koy6+waYm62Ucr9mgie7QVWdqSHE0gz92C0qhEvNFtkLgk3HFJxCQkkhgXRWJsBK3jIkiMjaR1bASJsRHERLiDO7J2M2Lz2RhjmobwKGg/mKj2g6kw0I6qMwjooR1weDfk7iL80C6SDmTS6tBO+uXuIqxwLRGew872Jb7lkPPRoy4OEE+OJnCQOPZqLD9oLIeIJV/i8ES2oCwiAW9kC4huiSumFeGxiUTEJRIfG0OL6HBaRIeT4Hs9sjTqab6bOEs2xpiTTwTikpzlyCog0rcc5SmGghwnMRVkO6/52bjzs2mZl0Xsob2kFhyAooO4incQXppLeFkheHCWAuBAxVMXagSHiCVPo8kjmj0axWac90USgycshpKwOLzhsXjD49DIeCQyDldkPK7oBMJjEgiPTiAyNoG4qAhiI8OIO7JEhREbGUZMuJuw4xnFoRmzZGOMabzCIp2hcxLaV1gtQLhvqcJTAkWHoOggFB4s9/4AnoIDaN5+YvIPEFl4mFZFuVCSh6vkEO7SnYR58oksy8ddUubUpOpQrOEUEkEBkRRqJPuJYCeRFGgkxRJJqSuSUlcUpa5oysKi8Lqj8IZF4w2LhfBoJDwaiYhBwqMJi4wiLCKasMhoIiKiCI+KIdK3LjwyiojwSCIj3ES4XUSGu3yvbiLDXIS5pNE3HQY02YjIBOAhwA08o6r3VyoXX/n5OH+D3KCq39e2r4gkAv8BugBbgStVtdLfLsaYkBUWUaXWdLQIP370VJ0aVUkeFB92lpI8KM6DksOUFhyitOAwJQW5eIoL8BTnQ3E+kaUFRJYU0KK0ELenEFfZIcLKipzFW0REURHhlB731ypToZgISgijmHAKNJwDhFNCOMWEUyrheCSCUomgzOW897oiUFc46gpD3eHgigB3BLjLfw5Hwnyv7gjEHYErLAIJa9iRIgKWbETEDTwGnANkAktF5F1VXVtus4lAT98yAngCGFHHvjOBeap6v4jM9H2+I1DfwxgTYkSce03hURDbpkrxkRpVzPEcu8wDnkIoLYSSfOe1tAA8RZSVFFFSUkhJUSGlxYV4igvxlBbhLSnCW1qEt7QYPEWox3mlrBjxlBBVVkRMWQlSVoy7rASXtwC3txi3twR3WSluTxluLSVMPbjxEEbZCV6g4xPIms1wYJOqbgEQkdeByUD5ZDMZeEmdLnGLRKSliLTDqbXUtO9k4Azf/i8C87FkY4xpCtxh4I6vdvoGNxDtWwLK6wVvKZSVQlmJ8+p13ntLS/B4SvCUFuMpKYa7T2+w0wYy2XQAdpT7nIlTe6lrmw517JuiqrsBVHW3iCRXd3IRmQ5M930sFpE1x/MlmqE2QHawg2gk7FocY9fiGLsWx/RuqAMFMtlUd7eq8kM9NW3jz761UtXZwGwAEUlvqL7iTZ1di2PsWhxj1+IYuxbHiEiDPaAYyL55mUDHcp9TgV1+blPbvnt9TW34Xvc1YMzGGGMCIJDJZinQU0S6ikgEMAV4t9I27wLXi2MkcMjXRFbbvu8CU33vpwLvBPA7GGOMaQABa0ZTVY+I3AZ8gnPv6zlVzRCRGb7yJ4EPcbo9b8Lp+nxjbfv6Dn0/MEdEpgHbgSv8CGd2w32zJs+uxTF2LY6xa3GMXYtjGuxahMTYaMYYY4LLxlMwxhgTcJZsjDHGBFyzSDYislVEVovIiiNd9UQkUUQ+E5GNvtdW5bafJSKbRGSDiJwXvMgbnu/B2DdFZL2IrBORUaF4LUSkt+/fw5ElV0R+FYrXAkBEfi0iGSKyRkReE5GoEL4Wt/uuQ4aI/Mq3LiSuhYg8JyL7yj93eDzfXURO9f3mbhKRh8WfgdlUtckvOGOktam07u/ATN/7mcDffO/7AStxBpftCmwG3MH+Dg14LV4EfuZ7HwG0DNVrUe6auIE9QOdQvBY4D0n/CET7Ps8BbgjRa3EKsAZntJkw4HOc4bJC4loAY4GhwJpy6+r93YElwCicZyI/AibWde5mUbOpwWScH158rxeXW/+6qhar6o84PeGGByG+BiciCTj/mJ4FUNUSVT1ICF6LSsYDm1V1G6F7LcKAaBEJw/mh3UVoXou+wCJVLVBVD/AVcAkhci1U9Wtgf6XV9fruvucbE1T1O3Uyz0vl9qlRc0k2CnwqIst8w9RApWFtgCPD2tQ0RE5z0A3IAp4XkeUi8oyIxBKa16K8KcBrvvchdy1UdSfwD5xHBXbjPM/2KSF4LXBqNWNFpLWIxOA8etGR0LwWR9T3u3fwva+8vlbNJdmcpqpDcUaRvlVExtay7QkPhdOIheFUkZ9Q1SFAPk61uCbN+VoA4Hso+CLgjbo2rWZds7gWvjb4yThNIe2BWBG5rrZdqlnXLK6Fqq4D/gZ8BnyM00zkqWWXZnst/NCgw4k1i2Sjqrt8r/uAt3CquTUNa+PPMDpNVSaQqaqLfZ/fxEk+oXgtjpgIfK+qe32fQ/FanA38qKpZqloKzAVGE5rXAlV9VlWHqupYnCaljYTotfCp73fP9L2vvL5WTT7ZiEisiMQfeQ+ci1NVrmlYm3eBKSISKSJdcW4OLjm5UQeGqu4BdojIkZFax+NMyxBy16KcqznWhAaheS22AyNFJMbXa2g8sI7QvBaIb6R4EekEXIrz7yMkr4VPvb67r6ntsIiM9P17uh5/hg0Ldu+IBuhd0Q2nKrwSyAD+6FvfGpiH81fLPCCx3D5/xOlZsQE/elE0pQUYDKQDq4C3gVYhfC1igBygRbl1oXot7gbW4/wh9m+cHkahei2+wfkjbCUwPpT+XeAk1t1AKU4NZdrxfHcgzfdvaTPwKL7RaGpbbLgaY4wxAdfkm9GMMcY0fpZsjDHGBJwlG2OMMQFnycYYY0zAWbIxxhgTcJZsjDHGBJwlG2OMMQFnycaYehCRaBH5SkTcx7FvhIh87Rt5uc5ji0iKiLwqIlt8g8x+JyKX1HGO+ZXnXPHN4/N4Xec3JpAs2RhTPz8F5qpqWX13VNUSnCe0r6rr2L5hQN4GvlbVbqp6Ks7o1ak17HvEa77typsCvObH+Y0JGEs2xviIyBsi8qiIfCsi20RkjIi8JCI/iMizvs2uxTcOlIj0FGeW2B6+z+EislJETq9hfSpOArm2hhCOHhs4CyhR1SePFKrqNlV9pFy814nIEnFmIn3KVyN6E5gkIpG+bbrgjPT8rW+32s5vTMBYsjHmmAHAFlUdgzOJ1LPAHTizO17q+wHvpqpbAVR1IzAbONJsdRvwjqp+U8P6TJzxpIZVPrFvKoSjxwb6A9/XFKiI9MWpoZymqoOBMuBaVc3BGShygm/TKcB/9Ni4VNWe35hAs7ZbYwARicKZQvtB36pC4Fn1TSolIgVAG+BgpV3XAGeLSCLOoIYjalvvayIrEZF4VT1c7jjVHbt8fI8BY3BqO8NwRm4+FVjqm/49mmNDwx9pSnvH9/rTI8ep5fzGBJTVbIxx9MeZ98br+zwIWAzga/7aBRQAUZX2+wHoDdwF/ENV8+tYD86Iy0WVjlNY6dgZOHMRAaCqt+IkmCTfKgFeVNXBvqW3qt7lK3sbGC8iQ4FoVa1cQ6ru/MYElCUbYxwDcIacP2IgzjQN4CSeVap6AHD7akFHbMZJCsNx5mKvdb2ItAaOTGJ2VDXH/gKIEpGby20WU+79PODycnOzJIpIZ9+x8oD5wHNUnMunxvMbE2iWbIxxDABWwNEmtWhfAoCKiedTnOYsAHw/2rnAzHK1ohrXA2cCH9YQw9Fj++6xXAyME5EfRWQJzn2kO3zla4E/AZ+KyCqcaY7blTvWazhJ8vVK56jt/MYEjM1nY0w9iMgQ4Deq+pNy67YDnbXSf6bq1ovIXGCWqm7w59gBiL/G8xsTSFazMaYeVHU58GW5By+7ANuqSTRV1vt6nL1d0w995WM3tLrOb0wgWc3GGGNMwFnNxhhjTMBZsjHGGBNwlmyMMcYEnCUbY4wxAWfJxhhjTMBZsjHGGBNwlmyMMcYE3P8H1L5Kgzo7yfIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from scipy.stats import expon, norm\n", "bkg = expon.freeze(loc=500, scale=100)\n", "sig = norm.freeze(loc=750, scale=30)\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "m = np.linspace(500, 1000, 1000)\n", "plt.plot(m, bkg.pdf(m), label='$H_0$')\n", "\n", "def model(m, fsig):\n", " \"\"\"normalized pdf for the mass distribution, fsig = signal fraction (0 <= fsig <= 1)\"\"\"\n", " return (1-fsig)*bkg.pdf(m) + fsig*sig.pdf(m)\n", "plt.plot(m, model(m, 0.1), label='$H_1$')\n", "\n", "plt.ylim(0, 0.010)\n", "plt.xlim(500, 1000)\n", "plt.legend()\n", "plt.xlabel('$m(\\gamma\\gamma)$ (GeV)')\n", "plt.ylabel('dp/dm (1/GeV)');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step 2: We create a function that returns the test statistic $\\ln Q$ for a given dataset." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "def LL(evts, nsig):\n", " \"\"\"\n", " returns log-likelihood for a model with nsig signal events,\n", " i.e., for a model with a signal fraction of nsig / #events\n", " \"\"\"\n", " fsig = nsig / len(evts)\n", " return np.sum(np.log(model(evts, fsig)))\n", "\n", "def lnQ(evts, nsig):\n", " \"\"\"returns test statistic\"\"\"\n", " return LL(evts, nsig) - LL(evts, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step 3: We generate 1000 toy Monte Carlo (MC) simulations for $H_0$ and for $H_1$. For each simulation, we generate 1000 events. We plot one simulation for $H_0$ and one for $H_1$ in the same plot." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# number of pseudo data sets\n", "ntoys = 1000\n", "\n", "# list for the value of the test statistic for each pseudo data set\n", "lnQ_B = []\n", "\n", "# generate pseudo data sets (background only model)\n", "for i in range(ntoys):\n", "\n", " # draw data points from background distribution\n", " toy = bkg.rvs(1000)\n", " \n", " # test statistic (assuming 30 signal events in H1)\n", " lnQ_B.append(lnQ(toy, 30))\n", "\n", "# list for the value of the test statistic for each pseudo data set \n", "lnQ_SB = []\n", "\n", "# generate pseudo data sets (signal + background model)\n", "for i in range(ntoys):\n", " \n", " # generate data points (970 background events, 30 signal events)\n", " toyB = bkg.rvs(1000-30)\n", " toyS = sig.rvs(30)\n", " toy = np.concatenate((toyS, toyB))\n", " \n", " # test statistic (assuming 30 signal events in H1)\n", " lnQ_SB.append(lnQ(toy, 30))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAd60lEQVR4nO3de5hcVZnv8e/PQNIJBxJIAhPShE40wBDDBIg4CiICKigMl0EFUZkDx+BIzigzz3MIlzOTcQ4jjiLKAYQwAXEeTLhfnKCCcAaEB4EOl1wAJUATCiKJAcIlF0h4zx97d1Jpqruuuy5dv8/z1NO1V+2199ubUG/vtfZaSxGBmZlZpT7Q6ADMzKy1OZGYmVlVnEjMzKwqTiRmZlYVJxIzM6vKNo0OoBpjxoyJrq6uRodhZtZSFi5c+KeIGFur47V0Iunq6qK7u7vRYZiZtRRJL9TyeG7aMjOzqjiRmJlZVZxIzMysKi3dR2JmlqV3332XXC7H+vXrGx1KRTo6Oujs7GTbbbfN9DxOJGZm/cjlcmy//fZ0dXUhqdHhlCUiWL16NblcjokTJ2Z6LjdtmZn1Y/369YwePbrlkgiAJEaPHl2XuyknEjOzAbRiEulVr9idSMzMrCpOJGZmVpWW7mx/+o9v0jVrQdn1xo8azgOzDs0gIjOz9tPSieTdTe/Rc8Hny65XSfIxMzvwgnt46fV1NTteqX/UXnHFFTzxxBNcdtllm8umTJnCTTfdxF577VWzeCrV0onEzKyeXnp9XUV/vPan1D9qFy1axL777rt5e/369SxfvpzJkyfXLJZquI/EzKzJLV68mP3222+r7T322IMhQ4Y0MKotfEdiZtbkli5dyvHHH7/5cd633nqLo446qsFRbeFEYmbWxF588UXGjh3L008/vbls5syZTJo0ibfffptvfvObDB06lEMOOYSTTz65ITG6acvMrIktWrSIKVOmbFX25JNPMnXqVG6++WZOOOEErrzySm6//fYGRehEYmbW1BYvXszee++9VdnSpUvZZ599yOVy7LbbbgAN7S9x05aZWYnGjxpe0+ED40cNL7rP4sWLt+oPefXVV4kIdtllFzo7O8nlckybNo333nuvZnGVy4nEzKxEjRjIfO211261vdNOO7Fy5UoAjj/+eGbOnMmCBQs4+uij6x5bLycSM7MWtd1223H11Vc3Ogz3kZiZWXUySySSrpK0UtKSvLLrJD2evnokPZ6Wd0lal/fZ5VnFZWZmtZVl09ZPgUuAn/UWRMSXet9LuhBYk7f/sxExLcN4zMwsA5klkoi4T1JXoc+UDM/8IuApeM3MWlyj+kg+AbwSEc/klU2U9JikeyV9or+KkmZI6pbUvWntmv52MzOzOmnUU1snAfPytlcAEyJitaT9gVslTYmIN/pWjIg5wByAYeMmR12iNTOzftX9jkTSNsDxwHW9ZRGxISJWp+8XAs8Ce9Q7NjMzK18jmrYOB56OiFxvgaSxkoak7ycBk4HnGhCbmZmVKcvHf+cBDwJ7SspJOi396ES2btYCOBhYJOkJ4EbgGxHxalaxmZlZ7WT51NZJ/ZT/TYGym4CbsorFzKwmLpoKa5bX7ngjJ8CZi4vu5qV2zcwGizXLYXYNnxadPbKk3bzUrpmZVcVL7ZqZWVW81K6ZmVVsoKV2n3vuOc4//3zWrFnDjTfe2LAY3bRlZtbEBlpqd9KkScydO7dBkW3hRGJm1sQGWmq3Wbhpy8ysVCMnlPykVcnHK2KgpXabRUsnkj21vKL/qPcPGwN8vvYBmdngVsKYj1obaKnd1atXc+655/LYY4/x3e9+l7PPPrvu8UGLJ5KhbKzome7OWv5FYWbWIKNHj+byyxu/DqD7SMzMrCpOJGZmVhUnEjMzq4oTiZnZACJad/28esXuRGJm1o+Ojg5Wr17dkskkIli9ejUdHR2Zn6uln9oyM8tSZ2cnuVyOVatWNTqUinR0dNDZ2Zn5eZxIzMz6se222zJx4sRGh9H03LRlZmZVcSIxM7OqZLlm+1WSVkpaklc2W9JLkh5PX5/L++xsScsk/V7SZ7OKy8zMaivLO5KfAkcUKL8oIqalrzsAJO0NnAhMSetcJqk5lv4yM7MBZdbZHhH3SeoqcfdjgPkRsQF4XtIy4ADgwSxiy8WY8ufbGjmhIRO2mZk1u0Y8tTVT0teAbuAfIuI1YDzwu7x9cmnZ+0iaAcwAmDBSFQVw0IaL6bmgzNl/PdGjmVlB9e5s/wnwQWAasAK4MC0vlBEKjgCKiDkRMT0ipo8dUVkiMTOz2qlrIomIVyJiU0S8B1xJ0nwFyR3Ibnm7dgIv1zM2MzOrTF0TiaRxeZvHAb1PdN0OnChpmKSJwGTg4XrGZmZmlcmsj0TSPOAQYIykHPBPwCGSppE0W/UApwNExFJJ1wNPAhuBMyJiU1axmZlZ7WT51NZJBYrnDrD/+cD5WcVjZmbZ8Mh2MzOrihOJmZlVxYnEzMyq4kRiZmZVcSIxM7OqOJGYmVlVnEjMzKwqTiRmZlaVtlyzffyo4XTNWlBWnZ6OjIIxM2txRROJpE6SRac+AewKrCOZI2sB8Mt0AsaW8sCsQ8uvNLvmYZiZDQoDJhJJV5OsC/KfwPeAlUAHsAfJSobnSpoVEfdlHaiZmTWnYnckF0bEkgLlS4CbJQ0FJtQ+LDMzaxXFEskRkl6PiFyhDyPiHWBZ7cNqTuX2q0DSH1NRU5qZWYsolkjGAw9Keh6YB9wQEX/KPqwmNHICPXy57Gq5dWOAZ2sfj5lZk1BEwRVtt+wgCTiYpMP9GOAJkqRyS0S8mXmEA5i+65DofrnJly2ZPRJmr2l0FGZmm0laGBHTa3W8ouNIInFvRPwtyXK4PwLOBF6pVRBmZta6Sh5HImkqyV3Jl4DVwDlZBWVmZq2j2OO/k0mSx0nAJmA+8JmIeK4OsZmZWQsodkfya5L+kC9FxOI6xGNmZi1mwD6SiJgUEedGxGJJu0s6HEDScEnbD1RX0lWSVkpaklf2fUlPS1ok6RZJo9LyLknrJD2evi6vxS9nZmbZK2nSRklfB24ErkiLOoFbi1T7Kcno93x3AR+OiH2APwBn5332bERMS1/fKCUuMzNrvFJn/z0DOBB4AyAingF2HqhCOm3Kq33K7oyIjenm70gSkpmZtbBSE8mGdBQ7AJK2AQYegFLcqcAv87YnSnpM0r2SPtFfJUkzJHVL6l61ttoQzMysWqUmknslnQMMl/Rp4AbgF5WeVNK5wEbg2rRoBTAhIvYF/h74uaQdCtWNiDkRMT0ipo8doUpDMDOzGik1kcwCVgGLgdOBO4DzKjmhpFOAo4CTIx1WHxEbImJ1+n4hyZwie1RyfDMzq6+SBiSma45cmb4qJukI4CzgkxGxNq98LPBqRGySNAmYDHisiplZCxjwjkTSMZLOyNt+SNJz6esLRerOAx4E9pSUk3QacAmwPXBXn8d8DwYWSXqC5Omwb0TEqwUPbGZmTaXYHcn/IhnZ3msY8BFgO+Bqkr6SgiLipALFc/vZ9ybgpiKxmJlZEyqWSIZGxIt52/enfRmrJW2XYVxmZtYiinW275i/EREz8zbH1j4cMzNrNcUSyUPpqPatSDodeDibkMzMrJUUa9o6E7hV0peBR9Oy/Un6So7NMjAzM2sNAyaSiFgJfFzSocCUtHhBRNyTeWRmZtYSSh1Hcg/g5GFmZu9T6sh2MzOzgpxIzMysKiUnEkk7Sdqx+J5mZtZOik2RMkHSfEmrgIeAR9JVD+dL6qpHgGZm1tyK3ZFcB9wC/FlETI6IDwHjSFZHnJ91cGZm1vyKJZIxEXFdRGzqLYiITRExHxidbWhmZtYKij3+u1DSZcA1QO+cW7sBpwCPZRmYmZm1hmKJ5GvAacA/A+MBkSSUX9DPTL5mZtZeio1sfwf4SfoyMzN7n6KP/0r6rKTTJO3ep/zU7MIyM7NWUezx338FzgWmAvdI+p95H88sXMvMzNpJsTuSo4FDI+LbJLP+HinpovQzZRqZmZm1hGKd7dtExEaAiHhd0tHAHEk3AEMHqijpKuAoYGVEfDgt24lkbEoX0AN8MSJeSz87m6RjfxPwdxHx60p/qWbTNWtBWfuPHzWcB2YdmlE0Zma1VSyRPCvpkxFxLyRjSIDTJP0f4K+L1P0pcAnws7yyWcDdEXGBpFnp9lmS9iZZG34KsCvwG0l75I9faWU9F3y+rP3LTTxmZo1UrGnrCxRYCTEiziMZT9KviLgPeLVP8TEkY1JIfx6bVz4/IjZExPPAMuCAIrGZmVkTKJZIdomIdYU+iIiXlOgs43y7RMSKtP4KYOe0fDxbBjwC5NKy95E0Q1K3pO5Va6OMU5uZWRaKNW19X9IHgNuAhcAqoAP4EPAp4DDgn0i++KtRqOO+YJaIiDnAHIDpuw5xJjEza7BiAxK/kPZfnAycSjJh41rgKeAO4PyIWF/G+V6RNC4iVkgaB6xMy3Ns3VTWCbxcxnHNzKxBii61GxFPkowlqYXbSebpuiD9eVte+c8l/ZCks30yBfpmzMys+ZS0ZnslJM0DDgHGSMqRNIFdAFwv6TRgOUlnPhGxVNL1wJPARuCMwfLElpnZYJdZIomIk/r56LB+9j8fOD+reMzMLBvFpkjJLNGYmdngUCxR/C5tlvoV8KuI6Mk+JDMzayXFntqans76eyTwI0njgfuBXwL3RsSGOsRoZmZNrOg08hHxQkRcHhHHAh8nWdTqcOC3kjyXh5lZmyurDyQi3gXuSV+kdyhmZtbGit6RDCQiXqpVIGZm1pqqSiRmZmZlJxJJH5C0QxbBmJlZ6ympj0TSz4FvkCw6tRAYKemHEfH9LIMbFEZOgNkjy6py/7AxQHlrmJiZNUqpne17R8Qbkk4mmazxLJKE4kRSzJmLy67SWWbiMTNrpFKbtraVtC3JQlS3pU9vmZmZlZxIriBZY3074L50kOKarIIyM7PWUWoi+UVEjI+Iz0VEkMzce2qGcZmZWYsoNZHclL+RJpP5tQ/HzMxazYCd7ZL2AqaQPKV1fN5HO5AsuWtmZm2u2FNbewJHAaOAo/PK3wS+nlVQZmbWOorN/nsbcJukj0XEg3WKyczMWkip40iWSToH6MqvExHucDcza3OlJpLbgN8CvyEZ3V4xSXsC1+UVTQL+kaT57OvAqrT8nIi4o5pztapcjCl/UOLICRUNfjQzq1apiWRERJxVixNGxO+BaQCShgAvAbcA/x24KCJ+UIvztLKDNlxMzwVlTpHi0fBm1iClPv77n5I+l8H5DwOejYgXMji2mZnVQamJ5FskyWS9pDckvSnpjRqc/0RgXt72TEmLJF0laccaHN/MzDJWUiKJiO0j4gMR0RERO6TbVU0lL2ko8FfADWnRT4APkjR7rQAu7KfeDEndkrpXrY1qQjAzsxooKZEo8RVJ/zvd3k3SAVWe+0jg0Yh4BSAiXomITRHxHnAlUPD4ETEnIqZHxPSxI1RlCGZmVq1SO9svA94DDgX+BXgLuBT4SBXnPom8Zi1J4yJiRbp5HLCkimO3pa5ZC8quM37UcB6YdWgG0ZhZuyg1kXw0IvaT9BhARLyWNk1VRNII4NPA6XnF/yZpGhAkMw2fXqBqWxg/anjZSaGng/Kf9KKy5GNmlq/URPJu+qhuAEgaS3KHUpGIWAuM7lP21UqPN9hUdIcwu+ZhmJmVpNSnti4mGeuxs6TzgfuBf80sKjMzaxkl3ZFExLWSFpKM+xBwbEQ8lWlkZmbWEkpKJJJ+DFwXEZdmHI+ZmbWYUpu2HgXOk7RM0vclTc8yKDMzax2lDki8JiI+RzK24w/A9yQ9k2lkZmbWEkq9I+n1IWAvkunkn655NGZm1nJKHdneewfyHZKBgvtHxNFFqpmZWRsodRzJ88DHIuJPWQZjZmatp9SmrTnAEZL+EUDShBrMtWVmZoNAqXckl7Jlrq3vAG8CN1HdXFtWSyMnVLS41f3DxgDlT61iZtarIXNtWQYqXGa37CV9zcz6KLVpq6ZzbZmZ2eBR6h1J37m2TgDOyywqq6tyZwD21PNmls9zbVnZ08976nkzy1fqHQkR8TQehGhmZn2UO7LdzMxsK04kZmZWFScSMzOrihOJmZlVpeTO9lqS1EMyOn4TsDEipkvaCbiOZGbhHuCLEfFaI+IzM7PSNfKO5FMRMS0iehfJmgXcHRGTgbvTbTMza3LN1LR1DHBN+v4a4NgGxmJmZiVqVCIJ4E5JCyXNSMt2iYgVAOnPnQtVlDRDUrek7lVro07hmplZfxrSRwIcGBEvS9oZuEtSyQMdI2IOybT2TN91iDOJmVmDNeSOJCJeTn+uJJnD6wDgFUnjANKfKxsRm5mZlafuiUTSdpK2730PfIZk+d7bgVPS3U4Bbqt3bGZmVr5GNG3tAtwiqff8P4+IX0l6BLhe0mnAcuALDYjNSvBgx7dg9pfLqrOCsYybvSyjiMyskeqeSCLiOeAvCpSvJpld2JrcOFbB7DXl1fECWmaDVjM9/mtmZi3IicTMzKriRGJmZlVxIjEzs6o0akCiNYuRE6DcjvCRE8o+TS7G0FlJh/vICXDm4vLrmVndOJG0uzp9SR+04eKy14YHyk9yZlZ3btoyM7Oq+I7EmlulTW9uDjOrGycSa26VJAQ3h5nVlZu2zMysKk4kZmZWFScSMzOrihOJmZlVxZ3tVhfjRw2na9aCiuo9MOvQDCIys1pxIrG6qDQZVJJ8KnbRVFizvLw6ftTYzInEbLM1y8teZ8WPGpu5j8TMzKrkRGJmZlWpe9OWpN2AnwF/BrwHzImIH0uaDXwdWJXuek5E3FHv+Ky5VNJJ39ORUTCFeAoXs4b0kWwE/iEiHpW0PbBQ0l3pZxdFxA8aEJM1qYo66WfXPIz+eQoXs/onkohYAaxI378p6SlgfL3jMDOz2mhoH4mkLmBf4KG0aKakRZKukrRjP3VmSOqW1L1qbdQpUjMz60/DEomk/wbcBHw7It4AfgJ8EJhGcsdyYaF6ETEnIqZHxPSxI1S3eM3MrLCGJBJJ25IkkWsj4maAiHglIjZFxHvAlcABjYjNzMzK04intgTMBZ6KiB/mlY9L+08AjgOW1Ds2s3KfEPMULmaNeWrrQOCrwGJJj6dl5wAnSZoGBNADnN6A2KzNlbuufF2ncDFrUo14aut+oFDnhseMWE3kYgydFTxim4sxdGYQj9lg57m2bNA5aMPFZd9ZABw0awE9tQ/HbNBzIjGrt0pGw/fW84h4a0JOJGZVqGQKl/GjfswDsysZse8R8dacnEjMqlDJE1vuoLfBxonEBp1qVmOsh0rjq+tklGZlUETrTjMyfdch0f3ypkaHYVYfs0eWv/CWWQGSFkbE9Fodz+uRmJlZVZxIzMysKk4kZmZWFXe2m9n7HHjBPbz0+rqy6njesfblRGJm7/PS6+s875iVzInErEWsYCzjyhyUmNRZVva57h/2dzD7y2XWGQOUPzWNtT4nErMWUUlCKDfx9OrUn8p+1LiSiTJtcHAiMbPWc9FUWLO8vDqeqywzTiRm1nrWLC9/cKbvmDLjRGI2iFXSr7K5Xpl1KloHxncJg4ITidkgVkm/ClB2EoEK14Gp512Cp+/PjBOJmTVONV/u5ao0GVQQX7uNw3EiMbPGGaR/6bfbOJymSySSjgB+DAwB/j0iLmhwSGbWZCr5i79SPR3lf8lXsiRBJWN3gKZoemuqRCJpCHAp8GkgBzwi6faIeLKxkZlZMZWus1Lpucruj6nURRPoocwv+GETgPK+3CsZuwM0xdNoTZVIgAOAZRHxHICk+cAxgBOJWZNr1fb9oir5a/+iqWV/wediDJ3lnwlofLNYsyWS8cCLeds54KP5O0iaAcxINzdIWlKn2JrdGOBPjQ6iSfhabOFrsUWTX4s34DuqsO5R5VbYs8ITFdRsiaTQVdxqCceImAPMAZDUXctVvlqZr8UWvhZb+Fps4WuxhaTuWh6v2dYjyQG75W13Ai83KBYzMytBsyWSR4DJkiZKGgqcCNze4JjMzGwATdW0FREbJc0Efk3y+O9VEbF0gCpz6hNZS/C12MLXYgtfiy18Lbao6bVQRBTfy8zMrB/N1rRlZmYtxonEzMyq0vSJRFKPpMWSHu99ZE3STpLukvRM+nPHvP3PlrRM0u8lfbZxkdeWpFGSbpT0tKSnJH2sTa/Dnum/hd7XG5K+3Y7XAkDSmZKWSloiaZ6kjja+Ft9Kr8NSSd9Oy9rmWki6StLK/LF1lfz+kvZPv3OXSbpYUvHBLRHR1C+gBxjTp+zfgFnp+1nA99L3ewNPAMOAicCzwJBG/w41ug7XAP8jfT8UGNWO16HPNRkC/BHYvR2vBckA3ueB4en29cDftOm1+DCwBBhB8hDRb4DJ7XQtgIOB/YAleWVl//7Aw8DHSMb1/RI4sti5m/6OpB/HkHyxkv48Nq98fkRsiIjngWUk0660NEk7kPwjmQsQEe9ExOu02XUo4DDg2Yh4gfa9FtsAwyVtQ/Il+jLteS3+HPhdRKyNiI3AvcBxtNG1iIj7gFf7FJf1+0saB+wQEQ9GklV+llenX62QSAK4U9LCdHoUgF0iYgVA+nPntLzQFCvj6xZpdiYBq4CrJT0m6d8lbUf7XYe+TgTmpe/b7lpExEvAD4DlwApgTUTcSRteC5K7kYMljZY0AvgcyeDmdrwW+cr9/cen7/uWD6gVEsmBEbEfcCRwhqSDB9i36BQrLWobklvWn0TEvsDbJLep/Rms12GzdMDqXwE3FNu1QNmguBZpe/cxJE0TuwLbSfrKQFUKlA2KaxERTwHfA+4CfkXSbLNxgCqD9lqUqL/fv6Lr0vSJJCJeTn+uBG4huf18Jb0FI/25Mt19sE6xkgNyEfFQun0jSWJpt+uQ70jg0Yh4Jd1ux2txOPB8RKyKiHeBm4GP057XgoiYGxH7RcTBJE08z9Cm1yJPub9/Ln3ft3xATZ1IJG0nafve98BnSG5hbwdOSXc7BbgtfX87cKKkYZImknS2PVzfqGsvIv4IvCipd8bOw0im1m+r69DHSWxp1oL2vBbLgb+UNCJ9suYw4Cna81ogaef05wTgeJJ/H215LfKU9funzV9vSvrL9N/U1/Lq9K/RTxoUeQphEskt6hPAUuDctHw0cDfJXxx3Azvl1TmX5AmE31PC0wat8gKmAd3AIuBWYMd2vA7p7zYCWA2MzCtr12vxz8DTJH9g/QfJUzjtei1+S/IH1hPAYe3274Ikca4A3iW5szitkt8fmJ7+e3oWuIR0BpSBXp4ixczMqtLUTVtmZtb8nEjMzKwqTiRmZlYVJxIzM6uKE4mZmVXFicTMzKriRGJmZlVxIjFLSRou6V5JQyqoO1TSfeksvEWPLWkXST+X9Fw6IemDko4rco7/6rtuRroWy2XFzm+WJScSsy1OBW6OiE3lVoyId0hGDn+p2LHTqSduBe6LiEkRsT/JTMad/dTtNS/dL9+JwLwSzm+WGScSawuSbpB0iaT7Jb0g6SBJP5P0B0lz091OJp1XSNJkJatzfijd3lbSE5I+0U95J0lyOLmfEDYfGzgUeCciLu/9MCJeiIj/mxfvVyQ9rGQVyCvSO5kbgaMkDUv36SKZ9ff+tNpA5zfLjBOJtYupwHMRcRDJAj9zgbNIVtY7Pv1ynhQRPQAR8QwwB+htSpoJ3BYRv+2nPEcyP9FH+p44nfJ+87GBKcCj/QUq6c9J7iwOjIhpwCbg5IhYTTKx4BHpricC18WWeY4Knt8sa25PtUFPUgfJ0sQ/SovWAXMjXfBH0lpgDPB6n6pLgMMl7UQyAd5HBypPm63ekbR9RLyZd5xCx86P71LgIJK7lI+QzOK7P/BIulz2cLZM/93bvHVb+vPU3uMMcH6zTPmOxNrBFJK1S95Lt/8CeAggbZJ6GVgLdPSp9wdgT2A28IOIeLtIOSSz767vc5x1fY69lGQ9GQAi4gyS5DE2LRJwTURMS197RsTs9LNbgcMk7UeyVnvfO5tC5zfLlBOJtYOpJFOL99qHZDp+SJLKooh4DRiS3r30epbkC/8AkrWrByyXNBroXWRqswLHvgfokPS3ebuNyHt/N3BC3voaO0naPT3WW8B/AVex9Xos/Z7fLGtOJNYOpgKPw+ZmruHplztsnVTuJGliAiD9Qn4DmJV3N9NvOfAp4I5+Yth87LRP41jgk5Kel/QwSb/NWennTwLnAXdKWkSyfOy4vGPNI0mA8/ucY6Dzm2XG65GYpSTtC/x9RHw1r2w5sHv0+R+lULmkm4GzI+L3pRw7g/j7Pb9ZlnxHYpaKiMeA/5c3aLALeKFAEnlfefpk1q39fYn3PXatFTu/WZZ8R2JmZlXxHYmZmVXFicTMzKriRGJmZlVxIjEzs6o4kZiZWVWcSMzMrCpOJGZmVpX/D4ixZE/cDdNvAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# generate an example of a background and a signal+background \n", "# pseudo data set and plot them together\n", "toyBonly = bkg.rvs(1000)\n", "\n", "toyB = bkg.rvs(1000-30)\n", "toyS = sig.rvs(30)\n", "toySB = np.concatenate((toyS, toyB))\n", "\n", "plt.hist(toyBonly, bins=25, range=(500, 1000), histtype='step', label='$H_0$');\n", "plt.hist(toySB, bins=25, range=(500, 1000), histtype='step', label='$H_1$');\n", "plt.legend()\n", "plt.xlabel('$m(\\gamma\\gamma)$ (GeV)')\n", "plt.ylabel('events / (20 GeV)')\n", "plt.xlim(500, 1000);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step 4: For each toy MC, we calculate the test statistic and fill it in two histograms, one for $H_0$ and one for $H_1$. Then we plot the histograms on the same figure to compare them." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAUBElEQVR4nO3df5DcdX3H8eebAxMQHCE5aMwJFzS0KqlRblSwdUIRpFYMdopBawnITOqIoNixJnY6dDqTkaltqUJtPSvmsKLJoA6RVjANiENVMMG08hvEFM6c5AdSxBIg4d0/7hs8wuXudr+7t3efPB8zmf3u98fu+3Obfd3nPvv9fjYyE0lSWQ7odAGSpNYz3CWpQIa7JBXIcJekAhnuklSgAztdAMDs2bOzt7e302VI0rSycePG7ZnZPdq2KRHuvb29bNiwodNlSNK0EhH/s69tDstIUoEMd0kqkOEuSQWaEmPuklTHM888w+DgIDt37ux0KW0xc+ZMenp6OOiggyZ8jOEuadobHBzksMMOo7e3l4jodDktlZns2LGDwcFB5s2bN+HjHJaRNO3t3LmTWbNmFRfsABHBrFmzGv6rxHCXVIQSg32PZtpmuEtSgRxzl1Scy9bd19LHu/jU48bdp6uriwULFpCZdHV1ccUVV3DSSSe1tI5GGO6aNOO94SbyBpKmqoMPPphNmzYBcMMNN7BixQpuvvnmjtXjsIwktdjjjz/O4Ycf3tEaxu25R8SVwDuArZl5fLXuU8AZwNPAT4DzMvOxatsK4HxgN3BRZt7Qptolacp48sknWbhwITt37mRoaIgbb7yxo/VMpOe+Cjh9r3XrgOMz87eB+4AVABHxauBs4DXVMZ+NiK6WVStJU9SeYZl77rmH66+/nnPOOYdOfkf1uOGemd8FHt1r3bczc1d19wdAT7W8GPhqZj6VmT8FHgDe0MJ6JWnKO/HEE9m+fTvbtm3rWA2tGHN/P/Ctanku8PCIbYPVuheIiGURsSEiNnTyByBJrXbPPfewe/duZs2a1bEaap0tExF/AewCvrxn1Si7jfp3SWb2A/0AfX19nfvbRVJxOnHm1Z4xdxieMmBgYICurs6NSjcd7hGxlOEPWk/JXw8sDQIvH7FbD7Cl+fIkaXrYvXt3p0t4nqaGZSLidODjwDsz8/9GbFoLnB0RMyJiHjAfuK1+mZKkRkzkVMivAIuA2RExCFzC8NkxM4B11ZwHP8jMD2TmnRGxBriL4eGaCzJzav06k6T9wLjhnpnvGWX1F8bYfyWwsk5RkqR6nH5AKtFNnxx7+8krJqcOdYzTD0hSgey5qyFjTf7lxF/S1GG4SyrPeMNSjZrAMNbKlSu5+uqr6erq4oADDuBzn/scb3zjG8c8ZtWqVXzsYx9j7ty5PPPMM7zqVa/iqquu4pBDDqldsuGulmn1HNrSdPH973+f6667jttvv50ZM2awfft2nn766eft09vby+bNm19w7JIlS7jiiisAeO9738vq1as577zzatdkuEtSTUNDQ8yePZsZM2YAMHv27IYfY9euXfzqV79q2VTBfqAqSTWddtppPPzwwxx33HF88IMfbOhLOlavXs3ChQuZO3cujz76KGeccUZLajLcJammQw89lI0bN9Lf3093dzdLlixh1apVXHDBBSxcuJCFCxeyZcuW55ZXrvz1pUBLlixh06ZN/PznP2fBggV86lOfaklNDstIUgt0dXWxaNEiFi1axIIFCxgYGOCb3/zmc9t7e3uf+xq+0UQEZ5xxBpdffjnLly+vXY89d0mq6d577+X+++9/7v6mTZs45phjGn6cW265hVe84hUtqcmeu6TyTPIVuE888QQXXnghjz32GAceeCCvfOUr6e/vn9Cxq1ev5pZbbuHZZ5+lp6eHVatWtaQmw12SajrhhBP43ve+N+Y+o50Gee6553Luuee2pSaHZSSpQIa7JBXIcJdUhF9/IVx5mmmb4S5p2ps5cyY7duwoMuAzkx07djBz5syGjvMDVUnTXk9PD4ODg2zbtq3TpbTFzJkz6enpaegYw13StHfQQQcxb968TpcxpTgsI0kFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgo07qmQEXEl8A5ga2YeX607AlgN9AKbgXdn5i+qbSuA84HdwEWZeUNbKpf2d63+EmgVZSI991XA6XutWw6sz8z5wPrqPhHxauBs4DXVMZ+NiK6WVStJmpBxwz0zvws8utfqxcBAtTwAnDli/Vcz86nM/CnwAPCGFtUqSZqgZsfcj8rMIYDq9shq/Vzg4RH7DVbrXiAilkXEhojYUOolw5LUKa3+QDVGWTfqTD6Z2Z+ZfZnZ193d3eIyJGn/1my4PxIRcwCq263V+kHg5SP26wG2NF+eJKkZzYb7WmBptbwUuHbE+rMjYkZEzAPmA7fVK1GS1KiJnAr5FWARMDsiBoFLgEuBNRFxPvAQcBZAZt4ZEWuAu4BdwAWZubtNtUuS9mHccM/M9+xj0yn72H8lsLJOUZKkerxCVZIKZLhLUoEMd0kqkOEuSQUy3CWpQH5BtqaFy9bdN+b2i089bpIqKcR4M0qevGJy6lDbGO56nvFCVNL04LCMJBXInrs0VfllHKrBnrskFchwl6QCOSwjdZJDL2oTe+6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQJ4KqSnDeW2k1rHnLkkFMtwlqUAOy0h6obGunHWu92nBnrskFchwl6QCGe6SVKBa4R4RF0fEnRFxR0R8JSJmRsQREbEuIu6vbg9vVbGSpIlpOtwjYi5wEdCXmccDXcDZwHJgfWbOB9ZX9yVJk6jusMyBwMERcSBwCLAFWAwMVNsHgDNrPockqUFNh3tm/gz4W+AhYAj438z8NnBUZg5V+wwBR452fEQsi4gNEbFh27ZtzZYhSRpFnWGZwxnupc8DXga8OCLeN9HjM7M/M/sys6+7u7vZMiRJo6gzLPNW4KeZuS0znwG+DpwEPBIRcwCq2631y5QkNaJOuD8EvCkiDomIAE4B7gbWAkurfZYC19YrUZLUqKanH8jMWyPiGuB2YBfwI6AfOBRYExHnM/wL4KxWFCpJmrhac8tk5iXAJXutforhXrykseZokdrIK1QlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoFqTfkrTRWXrbtvzO0Xn3rcJFUiTQ323CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVqFa4R8RLI+KaiLgnIu6OiBMj4oiIWBcR91e3h7eqWEnSxNTtuX8auD4zfwt4LXA3sBxYn5nzgfXVfUnSJGo63CPiJcBbgC8AZObTmfkYsBgYqHYbAM6sW6QkqTF1ZoU8FtgGfDEiXgtsBD4MHJWZQwCZORQRR452cEQsA5YBHH300TXKkMY31qyRzhjZoJs+Ofb2k1dMTh0aU51wPxB4PXBhZt4aEZ+mgSGYzOwH+gH6+vqyRh1S54wXdFKH1BlzHwQGM/PW6v41DIf9IxExB6C63VqvRElSo5oO98z8OfBwRPxmteoU4C5gLbC0WrcUuLZWhZKkhtX9JqYLgS9HxIuAB4HzGP6FsSYizgceAs6q+RySpAbVCvfM3AT0jbLplDqPK0mqxytUJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUB1v2ZP09Bl6+7rdAmS2syeuyQVyHCXpAIZ7pJUIMfcJbXWTZ/c97aTV0xeHfs5w13S5Bkr+MHwbyGHZSSpQLXDPSK6IuJHEXFddf+IiFgXEfdXt4fXL1OS1IhW9Nw/DNw94v5yYH1mzgfWV/clSZOo1ph7RPQAfwCsBD5arV4MLKqWB4DvAB+v8zxSO413UdfFpx43SZVIrVO35/4PwJ8Dz45Yd1RmDgFUt0eOdmBELIuIDRGxYdu2bTXLkCSN1HS4R8Q7gK2ZubGZ4zOzPzP7MrOvu7u72TIkSaOoMyzzZuCdEfF2YCbwkoj4V+CRiJiTmUMRMQfY2opCJUkT13TPPTNXZGZPZvYCZwM3Zub7gLXA0mq3pcC1tauUJDWkHRcxXQqsiYjzgYeAs9rwHBqDsz5Kakm4Z+Z3GD4rhszcAZzSiseVJDXHK1QlqUDOLSONZ7z5UKQpyJ67JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQU/5KNXz/wR1jbj/x2FmTVIn0fPbcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQVq+iKmiHg5cBXwG8CzQH9mfjoijgBWA73AZuDdmfmL+qVK7fGmh/rH3sELkTQN1blCdRfwZ5l5e0QcBmyMiHXAucD6zLw0IpYDy4GP1y9Ve1y27r5OlyBpimt6WCYzhzLz9mr5l8DdwFxgMTBQ7TYAnFm3SElSY1oyt0xE9AKvA24FjsrMIRj+BRARR+7jmGXAMoCjjz66FWVIbTHe/DFqoZs+Ofb2k1dMTh0FqP2BakQcCnwN+EhmPj7R4zKzPzP7MrOvu7u7bhmSpBFqhXtEHMRwsH85M79erX4kIuZU2+cAW+uVKElqVNPhHhEBfAG4OzP/fsSmtcDSankpcG3z5UmSmlFnzP3NwJ8AP46ITdW6TwCXAmsi4nzgIeCseiVKkhrVdLhn5i1A7GPzKc0+riSpPq9QlaQCGe6SVCDDXZIKZLhLUoFacoWqWs/5YyTVYc9dkgpkz137hXGn9W2T8ealOdHphBsz1twzzjvzPPbcJalAhrskFchhGRWhU8Mu0lRlz12SCmS4S1KBDHdJKpBj7h3iRUrSJNvPTqO05y5JBTLcJalADstIKsNYwy77IXvuklQge+5SB40194zzzqgOw72NPCNGUqc4LCNJBTLcJalAhrskFcgxd00LzvqothrvNMppeAWr4V6DH5i2lgH+fON9i9NYPNNGbRuWiYjTI+LeiHggIpa363kkSS/Ulp57RHQB/wicCgwCP4yItZl5VzueT9Lz+d2talfP/Q3AA5n5YGY+DXwVWNym55Ik7aVdY+5zgYdH3B8E3jhyh4hYBiyr7j4REffWeL7ZwPYax08VpbQDbMtUVEo7YNLb8ol2Pnidthyzrw3tCvcYZV0+705mP9CST9AiYkNm9rXisTqplHaAbZmKSmkH2JaJaNewzCDw8hH3e4AtbXouSdJe2hXuPwTmR8S8iHgRcDawtk3PJUnaS1uGZTJzV0R8CLgB6AKuzMw72/FclVJOkC6lHWBbpqJS2gG2ZVyRmePvJUmaVpxbRpIKZLhLUoGmbbhHxFkRcWdEPBsRfSPW90bEkxGxqfr3z52scyL21ZZq24pqCod7I+JtnaqxGRHxVxHxsxGvxds7XVMjSppCIyI2R8SPq9dhQ6fraUREXBkRWyPijhHrjoiIdRFxf3V7eCdrnKh9tKUt75NpG+7AHcAfAt8dZdtPMnNh9e8Dk1xXM0ZtS0S8muEzjV4DnA58tpraYTq5bMRr8e+dLmaiRkyh8fvAq4H3VK/HdHZy9TpMt/PDVzH8/3+k5cD6zJwPrK/uTwereGFboA3vk2kb7pl5d2bWuap1yhijLYuBr2bmU5n5U+ABhqd2UPs5hcYUkZnfBR7da/ViYKBaHgDOnNSimrSPtrTFtA33ccyLiB9FxM0R8budLqaG0aZxmNuhWpr1oYj47+rP0Wnxp3OlhJ/9SAl8OyI2VlN/THdHZeYQQHV7ZIfrqavl75MpHe4R8R8Rccco/8bqQQ0BR2fm64CPAldHxEsmp+J9a7It407j0GnjtOufgFcACxl+Xf6uo8U2Zsr/7Bv05sx8PcPDTBdExFs6XZCe05b3yZT+so7MfGsTxzwFPFUtb4yInwDHAR39EKmZtjANpnGYaLsi4vPAdW0up5Wm/M++EZm5pbrdGhHfYHjYabTPq6aLRyJiTmYORcQcYGunC2pWZj6yZ7mV75Mp3XNvRkR07/nQMSKOBeYDD3a2qqatBc6OiBkRMY/httzW4ZomrHrT7fEuhj84ni6KmUIjIl4cEYftWQZOY3q9FqNZCyytlpcC13awllra9T6Z0j33sUTEu4DLgW7g3yJiU2a+DXgL8NcRsQvYDXwgMyflA4xm7astmXlnRKwB7gJ2ARdk5u5O1tqgv4mIhQwPZ2wG/rSz5UxcB6bQaKejgG9EBAy/56/OzOs7W9LERcRXgEXA7IgYBC4BLgXWRMT5wEPAWZ2rcOL20ZZF7XifOP2AJBWouGEZSZLhLklFMtwlqUCGuyQVyHCXpAIZ7lJNEfGRiDik03VII3kqpFRTRGwG+jJze6drkfaw5679QkScU03M9F8R8aWIOCYi1lfr1kfE0dV+qyLij0Yc90R1uygivhMR10TEPRHx5Rh2EfAy4KaIuCkiuqrHuKOaP/3izrRY+7tpe4WqNFER8RrgLxiePGt7RBzB8DSxV2XmQES8H/gM408b+zqG59bfAvxn9XifiYiPMjxX+vaIOAGYm5nHV8/90jY1SxqTPXftD34PuGbPsEk1HcWJwNXV9i8BvzOBx7ktMwcz81lgE9A7yj4PAsdGxOURcTrweN3ipWYY7tofBONP17tn+y6q90UMT8byohH7PDVieTej/OWbmb8AXgt8B7gA+JemKpZqMty1P1gPvDsiZsHw928C32N4pkeAPwZuqZY3AydUy4uBgybw+L8E9sy6OBs4IDO/Bvwl8PoW1C81zDF3Fa+aXXMlcHNE7AZ+BFwEXBkRHwO2AedVu38euDYibmP4l8KvJvAU/cC3ImII+AjwxYjY03Fa0cKmSBPmqZCSVCCHZSSpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKtD/A9Z3enc4qoNSAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pltopts = dict(bins=40, range=(-15, 15), alpha=0.5)\n", "plt.hist(lnQ_B, **pltopts, label='B');\n", "plt.hist(lnQ_SB, **pltopts, label='S+B');\n", "plt.legend()\n", "plt.xlabel(\"test statistic lnQ\");\n", "plt.xlabel(\"counts\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step 5. For this exercise, we will use a test size $\\alpha = 5\\%$. We determine the critical test statistic and draw a line for where it is in the histogram above." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The critical test statistic is tc = 0.7538\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbQElEQVR4nO3df3RdZZ3v8feH8KNUikAbsDaUVH5UELxROzjgGlcZfohIhUGhuKbajjhVQbSMODZcUe9AL9wLXGChOETRFhFsh8pQEIS2tnqRXoFCxBYoFMSSaaS0WFBpkZbv/ePsHEI5JzlJztnPSfJ5rZWVc/bZe+ezm57zzfM8ez9bEYGZmRnATqkDmJlZ/XBRMDOzIhcFMzMrclEwM7MiFwUzMyvaOXWAgRgzZkw0NzenjmFmNqisXLlyY0Q0lnptUBeF5uZmHnzwwdQxzMwGFUm/L/eau4/MzKzIRcHMzIpcFMzMrGhQjymYmQ3Eq6++SkdHB1u3bk0dpSZGjBhBU1MTu+yyS8XbuCiY2bDV0dHBqFGjaG5uRlLqOFUVEWzatImOjg4mTJhQ8XbuPjKzYWvr1q2MHj16yBUEAEmMHj26z60gFwUzG9aGYkHo0p9jc1EwM7MijymYmWWuXPxEVfd33vGH9LpOQ0MDRxxxBBFBQ0MD3/rWtzj66KOrmqMvXBTMcnbfffcBJH3jW/3YfffdaW9vB+Duu++mtbWVX/ziF8nyuCiY5czFwMp56aWX2HvvvZNmcFEwy5lbCtbdli1baGlpYevWrXR2dvLzn/88aR4XBbOcXXDBBQAsX748bRCrC927j1asWMGnPvUpVq1aleysKJ99ZGZWJ4466ig2btzI888/nyyDi4KZWZ14/PHH2b59O6NHj06Wwd1HZmaZSk4hrbauMQUoTE0xb948Ghoacs/RxUXBzCyh7du3p47wBi4KZjm76qqrUkcwK8tFwSxnXV0FZvXIRcHqXm9TD6ToBx6IJUuWAHDcccclTmL2ZjU7+0jS9yVtkLSq27LLJD0u6RFJt0raq9trrZLWSloj6UO1ymWW2sUXX8zFF1+cOoZZSbU8JXUucOIOyxYDh0fEu4EngFYASYcBZwLvyra5VlK64Xczs2GqZkUhIn4JvLDDsnsiYlv29P8BTdnjU4AfR8QrEfE7YC1wZK2ymZlZaSnHFD4NzM8ej6NQJLp0ZMveRNJMYCbA+PHja5nPzIabZZdUd3/HtPa6ypw5c7jppptoaGhgp5124rrrruP9739/j9vMnTuXr3zlK4wbN45XX32VQw89lBtuuIGRI0cOOHKSK5ol/XdgG/CjrkUlVotS20ZEW0RMiohJjY2NtYpoZlZzK1as4I477uChhx7ikUceYcmSJey///5vWKe5ubnktlOnTqW9vZ3Vq1ez6667Mn/+/JLr9VXuLQVJ04GTgWMjouuDvwPo/i/RBKzPO5tZHq677rrUEaxOdHZ2MmbMGHbbbTcAxowZ0+d9bNu2jb/85S9Vm3I715aCpBOBrwIfjYiXu720CDhT0m6SJgAHA/fnmc0sLxMnTmTixImpY1gdOOGEE3j22Wc55JBDOPvss/t0c5358+fT0tLCuHHjeOGFF5gyZUpVMtXylNSbgRXAREkdks4CvgWMAhZLapf07wARsRpYADwK/Aw4JyLq69pvsyq5/fbbuf3221PHsDqwxx57sHLlStra2mhsbGTq1KnMnTuXc845h5aWFlpaWli/fn3x8Zw5c4rbdnUf/eEPf+CII47gsssuq0qmmnUfRcQnSiy+vof15wBzyr1uNlRcccUVAFX7y84Gt4aGBiZPnszkyZM54ogjmDdv3hv+aGhubi7eb6EUSUyZMoVrrrmG2bNnDziPp842M0tkzZo1PPnkk8Xn7e3tHHDAAX3ez7333suBBx5YlUye5sLMXtfbKZkVnGI5qOV8fH/+858599xz2bx5MzvvvDMHHXQQbW1tFW07f/587r33Xl577TWampqYO3duVTK5KFguepq/aLDNXWRWLe973/uK9+wu55lnnnnTshkzZjBjxoyaZHL3kZmZFbmlYMn1NgvqUPPDH/4wdQSzslwUzHK24xWrllZEIJWaVGHwe/364Mq5+8gsZ/Pnz6/alAQ2MCNGjGDTpk39+vCsdxHBpk2bGDFiRJ+2c0vBLGff+c53gMLFR5ZWU1MTHR0dPP/886mj1MSIESNoamrqfcVuXBTMbNjaZZddmDBhQuoYdcXdR2ZmVuSiYGZmRS4KZmZW5DEFs5zdcsstqSOYleWiYJaz/txIxSwv7j4yy9ncuXOrNnmZWbW5KJjlzEXB6pmLgpmZFbkomJlZkYuCmZkVuSiYmVmRT0k1y9mdd96ZOoJZWS4KZjkbOXJk6ghmZbn7yCxn1157Lddee23qGGYl1aylIOn7wMnAhog4PFu2DzAfaAaeAc6IiD9mr7UCZwHbgS9GxN21ymaW0oIFCwA4++yz0wRYdkman2uDQi1bCnOBE3dYNhtYGhEHA0uz50g6DDgTeFe2zbWSGmqYzczMSqhZUYiIXwIv7LD4FGBe9ngecGq35T+OiFci4nfAWuDIWmUzM7PS8h5T2C8iOgGy7/tmy8cBz3ZbryNb9iaSZkp6UNKDQ/UWemZmqdTLQLNKLCt5J+2IaIuISRExqbGxscaxzMyGl7xPSX1O0tiI6JQ0FtiQLe8A9u+2XhOwPudsZrlYvnx56ghmZeXdUlgETM8eTwdu67b8TEm7SZoAHAzcn3M2M7Nhr5anpN4MTAbGSOoAvgFcCiyQdBawDjgdICJWS1oAPApsA86JiO21ymaW0uWXXw7A+eefnziJ2ZvVrChExCfKvHRsmfXnAHNqlcesXtxxxx2Ai4LVp3oZaDYzszrgomBmZkUuCmZmVuRZUs1ytvvuu6eOYFaWi4INaVcufqLH1887/pCckrzurrvuyv1nVk1vk+kd05pPDqsZdx+ZmVmRWwpWFb39RW6vu+iiiwC48MILEycxezMXBbOcLV26FKhhUfD9EmwA3H1kZmZFLgpmZlbk7iOzwchdRFYjLgpmORs9enTqCGZluSiY5WzhwoWpI5iV5TEFMzMrclEwy1lrayutrb7y1+qTu4/McrZixYrUEczKckvBzMyK3FKwQc9TbJhVj1sKZmZW5JaCWc6amppSRzAry0XBLGc33nhj6gi109OV1r7XwqDg7iMzMytyUTDL2axZs5g1a1bqGGYlJek+knQe8BkggN8C/wSMBOYDzcAzwBkR8ccU+cxqqb29PXUEs7JybylIGgd8EZgUEYcDDcCZwGxgaUQcDCzNnpuZWY5SdR/tDOwuaWcKLYT1wCnAvOz1ecCpibKZmQ1buReFiPgv4HJgHdAJvBgR9wD7RURntk4nsG/e2czMhrvcxxQk7U2hVTAB2Az8h6Rpfdh+JjATYPz48TXJaFZLhxxySOoIZmWlGGg+DvhdRDwPIOknwNHAc5LGRkSnpLHAhlIbR0Qb0AYwadKkyCmzWdW0tbWljmBWVooxhXXA30oaKUnAscBjwCJgerbOdOC2BNnMzIa13FsKEfFrSbcADwHbgIcp/OW/B7BA0lkUCsfpeWczy8PMmTMBtxisPiW5TiEivgF8Y4fFr1BoNZgNaU88UcGsrj1NF2FWQ76i2czMiioqCpK+JGlPFVwv6SFJJ9Q6nJmZ5avSlsKnI+Il4ASgkcK0FJfWLJWZmSVR6ZiCsu8nAT+IiN9kZw6ZWR+1tLSkjmBWVqVFYaWkeyhccNYqaRTwWu1imQ1dV111VeoIZmVVWhTOAlqApyPiZUmjKXQhmZnZEFLpmMLiiHgoIjYDRMQm4MraxTIbuqZNm8a0aRXP7GKWqx5bCpJGUJjFdEw2Z1HXOMKewNtrnM1sSOro6Egdways3rqPPgvMolAAVvJ6UXgJ+HYNc5mZWQI9FoWIuBq4WtK5EXFNTpnMzCyRigaaI+IaSUdTuFXmzt2W31CjXGZmlkBFRUHSD4EDgXZge7Y4ABcFsz466qijUkcwK6vSU1InAYdFhO9fYDZAl1ziye6sflV6Suoq4G21DGJmZulV2lIYAzwq6X4KU1wDEBEfrUkqs5xcubjnaazPO776t8782Mc+BsDChQurvm+zgaq0KHyzliHMhpNNmzaljmBWVqVnH/2i1kHMzCy9Ss8++hOFs40AdgV2Af4SEXvWKpiZmeWv0pbCqO7PJZ0KHFmTRGZmlky/7tEcEf8paXa1w5gNB8ce61uRW/2qtPvotG5Pd6Jw3YKvWTDrhwsvvDB1BLOyKm0pTOn2eBvwDHBK1dOYmVlSlY4p+IY6ZlXy4Q9/GIC77rorcRKzN6voimZJTZJulbRB0nOSFkpq6u8PlbSXpFskPS7pMUlHSdpH0mJJT2bf9+7v/s3q2ZYtW9iyZUvqGGYlVTrNxQ+ARRTuqzAOuD1b1l9XAz+LiHcC/w14DJgNLI2Ig4Gl2XMzM8tRpUWhMSJ+EBHbsq+5QGN/fqCkPYEPAtcDRMRfs9t8ngLMy1abB5zan/2bmVn/VVoUNkqaJqkh+5oG9Pda/XcAzwM/kPSwpO9JeguwX0R0AmTf9+3n/s3MrJ8qPfvo08C3gCspnIp6H9DfweedgfcC50bEryVdTR+6iiTNBGYCjB8/vp8RzCrT04R5/Z0s7+STT+5vnMFtWS9Thh/Tmk8O61GlReEiYHpE/BFA0j7A5RSKRV91AB0R8evs+S0UisJzksZGRKekscCGUhtHRBvQBjBp0iRfK2GDzvnnn1/4gOztQ9IsgUq7j97dVRAAIuIF4D39+YER8QfgWUkTs0XHAo9SGMieni2bDtzWn/2bmVn/VdpS2EnS3ju0FPo1RUbmXOBHknYFnqbQFbUTsEDSWcA64PQB7N+sbk2ePBk2r2P5lf+cOorZm1T6wX4FcJ+kWyiMKZwBzOnvD42IdgpTZezIk8KYmSVU6RXNN0h6EPh7QMBpEfFoTZOZmVnuKu4CyoqAC4GZ2RBW6UCzmZkNAwMZLDazfjjjjDPgibtTxzAryUXBLGdnn302LHsxdQyzktx9ZJazl19+mZe3/jV1DLOSXBTMcnbSSSdxUuu83lc0S8BFwczMilwUzMysyEXBzMyKXBTMzKzIp6Sa5WzGjBnw2B2pY5iV5KJglrMZM2bAss7UMcxKcveRWc42btzIxhf/kjqGWUkuCmY5+/jHP87Hv3lT6hhmJbkomJlZkYuCmZkVeaDZKnbl4idSRzCzGnNLwczMitxSMMvZ5z//eVj9n6ljmJXkomCWs6lTp8Kyp1PHqD/LLin/2jGt+eUY5tx9ZJazZ599lmc3bE4dw6wktxTMcvbJT34SNq9j+ZX/nDrK4NFTKwLckqgitxTMzKwoWVGQ1CDpYUl3ZM/3kbRY0pPZ971TZTMzG65Sdh99CXgM2DN7PhtYGhGXSpqdPf9qqnBmventuo3zjj8kpyRm1ZOkpSCpCfgI8L1ui08Bum5cOw84Ne9cZmbDXaqWwlXAvwKjui3bLyI6ASKiU9K+pTaUNBOYCTB+/Pha5zSrui9/+cvwyH+kjmFWUu4tBUknAxsiYmV/to+ItoiYFBGTGhsbq5zOrPamTJnClKMPTR3DrKQULYUPAB+VdBIwAthT0o3Ac5LGZq2EscCGBNnMam7NmjWw7nkmjvcfNVZ/ci8KEdEKtAJImgycHxHTJF0GTAcuzb7flne24c4T3uXjs5/9rK9TsLpVT9cpXAocL+lJ4PjsuZmZ5SjpFc0RsRxYnj3eBBybMo+Z2XDnaS7MaqXc1Ayb1+Wbw6wP6qn7yMzMEnNLwSxnX5t2TOoIZmW5KJjl7Lj3HZQ6gllZ7j4yy1n72vW0r12fOoZZSW4pmOVs1rd/CuDrFKwuuaVgZmZFLgpmZlbkomBmZkUuCmZmVuSBZrOc/c+zTkgdwawsFwWznB19+AGpI5iV5e4js5zdt+r33Lfq96ljmJXkloJZzi64/h7A1ylYfXJLwczMitxSMEvgpS2vsuLpTWVfP+odo3NMY/Y6txTMzKzIRcHMzIrcfWSWs6vO+Qi/6XgxdQyzklwUzHLWctDb2bLTbqljmJXk7iOznC1ZuZYHHvV1Claf3FIwy9nFNy7jpS2v8jeH+cpmqz9uKZiZWVHuLQVJ+wM3AG8DXgPaIuJqSfsA84Fm4BngjIj4Y975zCr1t+vael7B1xrYIJSi+2gb8OWIeEjSKGClpMXADGBpRFwqaTYwG/hqgnxD1pWLn0gdwczqXO7dRxHRGREPZY//BDwGjANOAeZlq80DTs07m5nZcJd0oFlSM/Ae4NfAfhHRCYXCIWnfMtvMBGYCjB8/Pp+gZv1QbhqLz50+Od8gw8GyS3p+/ZjWfHIMAckGmiXtASwEZkXES5VuFxFtETEpIiY1NjbWLqBZjRzwtn044G37pI5hVlKSoiBpFwoF4UcR8ZNs8XOSxmavjwU2pMhmVmv3tj/Fve1PpY5hVlLuRUGSgOuBxyLi/3R7aREwPXs8Hbgt72xmebj5nge5+Z4HU8cwKynFmMIHgE8Cv5XUni27ALgUWCDpLGAdcHqCbGZmw1ruRSEi7gVU5uVj88xiZmZv5CuazcysyEXBzMyKPCGeWc6+/pkPp45gVpaLwhDjqSzq33777Jk6gllZ7j4yy9mS+x9nyf2Pp45hVpJbCmY96HUm1H64dflvADjuyHeWXafcFBldjvIMrH3T0zQYngLjDdxSMDOzIhcFMzMrcveRDWu16B4yG8zcUjAzsyK3FMxyNufzU1JHMCvLRWGQ8XUIg99eo0amjmB9MczOXHL3kVnOfvqrVfz0V6tSxzAryUXBLGd3/mo1d/5qdeoYZiW5KJiZWZHHFMxseOtpzGAYclEwG4R6mgbDU2DYQLj7yMzMitxSqEM+7XRou+JLp6WOYFaWi4JZzkbstkvqCGZlufvILGcLl7WzcFl76hhmJbmlYENaPU549/MH1gDwsWNaEiexAevtzKVBeMWzWwpmZlZUdy0FSScCVwMNwPci4tLEkarOA8nVVY+tgZR6u2tbT3w6q9VVUZDUAHwbOB7oAB6QtCgiHk2bzGx48G1Ard66j44E1kbE0xHxV+DHwCmJM5mZDRt11VIAxgHPdnveAby/+wqSZgIzs6d/lrRmAD9vDLBxANvXi6FyHDCMjuXoz1yRY5QBGTa/k+q7oJY7H8ixHFDuhXorCiqxLN7wJKINqEonsqQHI2JSNfaV0lA5DvCx1KOhchzgY6lEvXUfdQD7d3veBKxPlMXMbNipt6LwAHCwpAmSdgXOBBYlzmRmNmzUVfdRRGyT9AXgbgqnpH4/Imp5N5Khci7jUDkO8LHUo6FyHOBj6ZUiove1zMxsWKi37iMzM0vIRcHMzIqGXVGQdLqk1ZJekzSp2/JmSVsktWdf/54yZyXKHUv2WquktZLWSPpQqoz9Iembkv6r2+/ipNSZ+kLSidm/+1pJs1PnGQhJz0j6bfZ7eDB1nr6Q9H1JGySt6rZsH0mLJT2Zfd87ZcZKlTmWmrxPhl1RAFYBpwG/LPHaUxHRkn19Ludc/VHyWCQdRuHMrXcBJwLXZlOIDCZXdvtd3Jk6TKW6TdXyYeAw4BPZ72MwOyb7PQy28/vnUvj/391sYGlEHAwszZ4PBnN587FADd4nw64oRMRjETGQq6DrRg/Hcgrw44h4JSJ+B6ylMIWI1Z6naqkTEfFL4IUdFp8CzMsezwNOzTVUP5U5lpoYdkWhFxMkPSzpF5L+LnWYASg1Xci4RFn66wuSHsmazYOiiZ8ZCv/23QVwj6SV2RQzg91+EdEJkH3fN3Gegar6+2RIFgVJSyStKvHV019sncD4iHgP8C/ATZL2zCdxef08ll6nC0mtl+P6DnAg0ELh9zJoJgliEPzb99EHIuK9FLrDzpH0wdSBrKgm75O6unitWiLiuH5s8wrwSvZ4paSngEOApINr/TkWBsF0IZUel6TvAnfUOE411f2/fV9ExPrs+wZJt1LoHis1HjdYPCdpbER0ShoLbEgdqL8i4rmux9V8nwzJlkJ/SGrsGoyV9A7gYODptKn6bRFwpqTdJE2gcCz3J85UsezN2uUfKAyoDxZDZqoWSW+RNKrrMXACg+t3UcoiYHr2eDpwW8IsA1Kr98mQbCn0RNI/ANcAjcBPJbVHxIeADwL/JmkbsB34XETkMrDTX+WOJSJWS1oAPApsA86JiO0ps/bR/5bUQqHb5Rngs2njVC7BVC21tB9wqyQofFbcFBE/SxupcpJuBiYDYyR1AN8ALgUWSDoLWAecni5h5cocy+RavE88zYWZmRW5+8jMzIpcFMzMrMhFwczMilwUzMysyEXBzMyKXBRsUJK0l6SzB7D9LEkj+7D+qZVMbLfjepL+TVLZC/X6uv4O2zZ3nzWzh/XeKukGSU9lXz8aZFOHWI5cFGyw2gvod1EAZgEVFwUKE6dVMtvpG9aLiK9HxJIqrt8f1wNPR8SBEXEghQkS51b5Z9gQ4esUbFCS1DX76BpgcUR8RdJXgDOA3YBbI+Ib2ZW4CyhMN9EAXEThoqzLs203RsQxO+z7UuCjFC78uwf4CYUpBF7Mvj4G/D0wE9iVwofsJynMQbPjehcCd0TELRXut/v6fwNcDbyFwhQsx0bEn7rlbM7WPVzSjGzfIynMh3NrRPyrpIOAxcBBXRcwZlfuPwV8aKjMGGzVM+yuaLYhYzZweES0AEg6gcJ0HkdSmJRuUTZ5WyOwPiI+kq331oh4UdK/ULhPwMbuO5W0D4UpA94ZESFpr4jYLGkR2Yd1tt7miPhu9vhi4KyIuKbEen3db9f6uwLzgakR8UA2OeOWXv5NWoD3UCggayRdQ6EV0t79ivaI2C7pYeBQCoXRrMjdRzZUnJB9PQw8BLyTQpH4LXCcpP8l6e8i4sVe9vMSsBX4nqTTgJfLrHe4pP8r6bfAP1K4oVE19ttlItAZEQ8ARMRLEbGtl22WRsSLEbGVwhQnB1AokKW6A0rN5mrmomBDhoBLut2F6qCIuD4ingDeR6E4XCLp6z3tJPvgPRJYSKG/v9xcP3OBL0TEEcD/AEZUab/dj6evfbuvdHu8nUJPwGrgPZKK7/Xs8bspFE+zN3BRsMHqT8Cobs/vBj4taQ8ASeMk7Svp7cDLEXEjhXGE95bZnmy7PYC3Zrc2nEWhS6bU+qOATkm7UGgplMvV1/12eRx4ezaugKRRkvrc3RsRaym0nr7WbfHXKLQq1vV1fzb0eUzBBqWI2CTpV9kpmXdlA82HAiuyfvk/A9OAg4DLJL0GvAp8PttFG3CXpM4dBppHAbdJGkHhr/XzsuU/Br4r6YvAxykMCP8a+D2FVsioMuv1db9dx/dXSVOBayTtTmE84bjsuPrq09l+1gJvpTC995R+7MeGAZ99ZDaMSJoI3AmcW60bvdvQ4qJgZmZFHlMwM7MiFwUzMytyUTAzsyIXBTMzK3JRMDOzIhcFMzMr+v9SfwLwRNzAewAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# test rule: the null hypothesis will be rejected if we observe t in the critical region t > tc\n", "# alpha is the probability to reject H0 although it is actually true\n", "\n", "tc = np.percentile(lnQ_B, 95)\n", "print(f'The critical test statistic is tc = {tc:.4f}')\n", "\n", "fig, ax = plt.subplots()\n", "plt.hist(lnQ_B, **pltopts, label='B');\n", "plt.hist(lnQ_SB, **pltopts, label='S+B');\n", "plt.legend()\n", "plt.plot([tc, tc], [0, ax.get_ylim()[1]], color='black', linestyle='dashed');\n", "plt.xlabel(\"test statistic lnQ\");\n", "plt.ylabel(\"counts\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step 6. Determine the power $1-\\beta$ of the test." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The power of the test is 1-beta = 0.913\n" ] } ], "source": [ "# beta is the probability to accept H0 although it is actually false (i.e. H1 is true)\n", "# a test has good performance if 1-beta is large\n", "\n", "# The Neyman-Pearson lemma tells us that the likelihood-ratio Q \n", "# is the test sttistic that for a given size alpha gives us the highest power\n", "oneminusbeta = np.sum(lnQ_SB>tc) / len(lnQ_SB)\n", "print(f'The power of the test is 1-beta = {oneminusbeta}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Now it is your turn!\n", "#### Step 7.: Read the data and plot them as a histogram. Evaluate the test-statistic for the data. What are the $p$-values of the data for both $H_0$ and $H_1$? To which significance do they correspond?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "d = pd.read_csv('https://www.physi.uni-heidelberg.de/~reygers/lectures/2020/smipp/two_photon_inv_masses.csv');\n", "data = d['m']" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Your solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Step 8: What is your conclusion on the existence of the new particle given the observed data?" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Your solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7.4 Best fit (10 points)\n", "\n", "Given is a set of 6 data points with an $x$ and $y$ value and an uncertainty in $y$ (see below). Implement three least square fits (a constant fit, a linear fit and a quadratic fit). Calculate the $\\chi^2$/dof for each of the fits as well as the p-value. Plot the data together with the fit results and comment on the quality of the fits" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from iminuit import Minuit\n", "from pprint import pprint\n", "\n", "xd = np.array([ 1., 2., 3., 4., 5., 6.])\n", "yd = np.array([1.1, 1.2, 1.3, 1.4, 2.0, 1.9])\n", "yd_err = np.array([0.1, 0.1, 0.3, 0.1, 0.5, 0.2])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }