CP2_Muon_decay_pseudo-experiemnt.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "39f739db",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Making general imports\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import random as rn\n",
    "import sklearn\n",
    "import tensorflow as tf\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.inspection import PartialDependenceDisplay\n",
    "from sklearn import metrics\n",
    "%matplotlib inline\n",
    "\n",
    "from sklearn import datasets\n",
    "from sklearn.tree import DecisionTreeClassifier, export_graphviz , plot_tree\n",
    "from sklearn.metrics import accuracy_score, confusion_matrix, r2_score\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "\n",
    "from sklearn import tree\n",
    "from sklearn.tree import DecisionTreeClassifier, plot_tree\n",
    "from sklearn.tree import export_graphviz\n",
    "from sklearn.model_selection import train_test_split # Import train_test_split function\n",
    "from sklearn import metrics # Import scikit-learn metrics module for accuracy calculation\n",
    "from sklearn import preprocessing # Import preprocessing for String-Int conversion\n",
    "from sklearn.compose import ColumnTransformer\n",
    "from sklearn.preprocessing import OneHotEncoder, StandardScaler, MinMaxScaler\n",
    "\n",
    "\n",
    "from sklearn.neural_network import MLPClassifier\n",
    "\n",
    "from tensorflow.python.keras.models import Model\n",
    "from tensorflow.python.keras.layers import Input, Dense\n",
    "from tensorflow.python.keras.regularizers import l1_l2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7336c28d",
   "metadata": {},
   "source": [
    "## Single experiment"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f761cbe7",
   "metadata": {},
   "source": [
    "<b>1. Simulate the decay times of 1000 muons<b/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "e41b8908",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Generating a random sample of 1500 muon decay times\n",
    "decay_times = np.random.exponential(scale=2.2e-6,size = 1500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "c9875d69",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Selecting muon decays with lifetimes between 0 and 10 microseconds\n",
    "decays_cut = np.where(decay_times<10e-6)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "11b815f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Selecting 1000 muon decays from the decays between 0 and 10 microseconds\n",
    "decay_sample = decay_times[decays_cut[:1000]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "079c7c80",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAHFCAYAAAAUpjivAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA13UlEQVR4nO3de1xVVf7/8fcREAEB8cZFUbDQ8G7eJrRBS/GW1piTo46jZY2Nppmpyc9UdArSJofStNExtBLtO2Vlo6V4Q02bUdJJ0bEyMTWJNAS8YcL+/dGDM55Qk+s5LF/Px2M/Hp6119n7c9YDOW/WXmcfm2VZlgAAAAxVzdkFAAAAVCTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOcItZtmyZbDabbDabtm7dWmy/ZVm6/fbbZbPZ1K1bt0qvryRGjhypsLAwh7awsDCNHDnSoW3v3r2Kjo6Wv7+/bDabEhMTb/ocFy5cUFxc3DXHqmgsMzIySlw7gMrj7uwCADiHr6+vli5dWizQpKam6siRI/L19XVOYWX03nvvyc/Pz6HtkUce0fnz57Vq1SoFBAQUC0g3cuHCBc2aNUuSio1Vv379tGvXLgUHB5e1bAAViLAD3KIGDx6sFStW6NVXX3UIB0uXLtVdd92l3NxcJ1ZXeu3atSvWduDAAT322GPq06dPuZ6rXr16qlevXrkeE0D54zIWcIsaMmSIJGnlypX2tpycHL377rt65JFHivXfunXrNS99ZWRkyGazadmyZQ7ta9as0V133SVvb2/5+vqqZ8+e2rVrl0OfuLg42Ww2paena8iQIfL391dgYKAeeeQR5eTklOp1XX0Zq+gy05UrV7Ro0SL75bsimZmZGj16tBo2bKjq1asrPDxcs2bN0pUrV+yvrSjMzJo1y/78nx//6stY3bp1U8uWLbVr1y5FRUXJy8tLYWFhSkpKkiStXbtWd955p7y9vdWqVSt9/PHHxV7Dl19+qaFDh6p+/fry9PRUZGSkXn31VYc+hYWFeu6559SsWTN5eXmpVq1aat26tV5++eVSjRtgMsIOcIvy8/PToEGD9Prrr9vbVq5cqWrVqmnw4MFlOnZycrLuv/9++fn5aeXKlVq6dKmys7PVrVs37dixo1j/Bx98UE2bNtW7776rqVOnKjk5WU899VSZapD+d5lJkgYNGqRdu3bZH2dmZqpTp05av369ZsyYoY8++kijRo1SQkKCHnvsMUlScHCwPYyMGjXK/vzp06ff8LyZmZl6+OGH9eijj+qDDz5Qq1at9Mgjj2j27NmKjY3VlClT9O6776pmzZp64IEH9O2339qfe/DgQXXs2FEHDhzQSy+9pH/+85/q16+fxo8fb7+cJklz585VXFychgwZorVr1+rtt9/WqFGjdPbs2TKPG2AcC8AtJSkpyZJk7d6929qyZYslyTpw4IBlWZbVsWNHa+TIkZZlWVaLFi2s6Oho+/OK+m7ZssXheEePHrUkWUlJSZZlWVZBQYEVEhJitWrVyiooKLD3y8vLs+rXr29FRUXZ22bOnGlJsubOnetwzDFjxlg1atSwCgsLb/haRowYYTVu3NihrXHjxtaIESMc2iRZY8eOdWgbPXq0VbNmTevYsWMO7X/5y18sSVZ6erplWZb1/fffW5KsmTNnFjt/0VgePXrU3hYdHW1Jsvbs2WNvO3PmjOXm5mZ5eXlZJ0+etLfv27fPkmS98sor9rZevXpZDRs2tHJychzO9cQTT1g1atSwfvjhB8uyLOu+++6z2rZte+2BAeCAmR3gFhYdHa3bbrtNr7/+uvbv36/du3df8xJWSRw+fFjffvuthg8frmrV/vcrpmbNmnrwwQf16aef6sKFCw7PGTBggMPj1q1b69KlS8rKyipTLTfyz3/+U927d1dISIiuXLli34rW9aSmppb62MHBwWrfvr39ce3atVW/fn21bdtWISEh9vbIyEhJ0rFjxyRJly5d0qZNm/Sb3/xG3t7eDnX17dtXly5d0qeffipJ6tSpk/7zn/9ozJgxWr9+fZVdYwVUBhYoA7cwm82mhx9+WK+88oouXbqkpk2b6u677y7TMc+cOSNJ1/yEUkhIiAoLC5WdnS1vb297e506dRz6eXp6SpIuXrxYplpu5LvvvtOHH34oDw+Pa+4/ffp0qY9du3btYm3Vq1cv1l69enVJP4Uc6aexu3LliubPn6/58+ffsK7Y2Fj5+Pjorbfe0muvvSY3Nzf9+te/1pw5c9ShQ4dS1w6YiLAD3OJGjhypGTNm6LXXXtPzzz9/3X41atSQJOXn5zu0/zwUFAWXU6dOFTvGt99+q2rVqikgIKCsZZdZ3bp11bp16+u+5qtnYCpLQECA3NzcNHz4cI0dO/aafcLDwyVJ7u7umjhxoiZOnKizZ89q48aN+n//7/+pV69eOn78uEOYBG51hB3gFtegQQNNnjxZ//3vfzVixIjr9iu6N83nn3+uXr162dvXrFnj0K9Zs2Zq0KCBkpOTNWnSJPunn86fP693333X/gktZ7vvvvu0bt063XbbbTcMX5Uxy1TE29tb3bt31969e9W6dWv7zM8vqVWrlgYNGqSTJ09qwoQJysjIUPPmzSu4WqDqIOwA0AsvvPCLfYKCgtSjRw8lJCQoICBAjRs31qZNm7R69WqHftWqVdPcuXM1bNgw3XfffRo9erTy8/P14osv6uzZszd1rsowe/ZspaSkKCoqSuPHj1ezZs106dIlZWRkaN26dXrttdfUsGFD+fr6qnHjxvrggw907733qnbt2qpbt26JbkxYEi+//LK6du2qu+++W3/6058UFhamvLw8ffXVV/rwww+1efNmSVL//v3VsmVLdejQQfXq1dOxY8eUmJioxo0bKyIiokJqA6oqwg6Am/bmm29q3LhxeuaZZ1RQUKD+/ftr5cqVxdaIDB06VD4+PkpISNDgwYPl5uamX/3qV9qyZYuioqKcVL2j4OBg7dmzR3/+85/14osv6sSJE/L19VV4eLh69+7tMNuzdOlSTZ48WQMGDFB+fr5GjBhR7L5C5aV58+b67LPP9Oc//1nPPvussrKyVKtWLUVERKhv3772ft27d9e7776rv//978rNzVVQUJB69uyp6dOnX3cdEnCrslmWZTm7CAAAgIrCR88BAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIzGfXYkFRYW6ttvv5Wvr6/9bq8AAMC1WZalvLw8hYSEOHzx8M8RdvTT9/WEhoY6uwwAAFAKx48fV8OGDa+7n7AjydfXV9JPg+Xn5+fkagAAwM3Izc1VaGio/X38egg7kv3SlZ+fH2EHAIAq5peWoLBAGQAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0d2cXAEdhU9fedN+MF/pVYCUAAJiBmR0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEZzatjZtm2b+vfvr5CQENlsNr3//vsO+y3LUlxcnEJCQuTl5aVu3bopPT3doU9+fr7GjRununXrysfHRwMGDNCJEycq8VUAAABX5tSwc/78ebVp00YLFiy45v65c+dq3rx5WrBggXbv3q2goCD17NlTeXl59j4TJkzQe++9p1WrVmnHjh06d+6c7rvvPhUUFFTWywAAAC7MqffZ6dOnj/r06XPNfZZlKTExUdOmTdPAgQMlScuXL1dgYKCSk5M1evRo5eTkaOnSpXrzzTfVo0cPSdJbb72l0NBQbdy4Ub169aq01wIAAFyTy67ZOXr0qDIzMxUTE2Nv8/T0VHR0tHbu3ClJSktL048//ujQJyQkRC1btrT3uZb8/Hzl5uY6bAAAwEwuG3YyMzMlSYGBgQ7tgYGB9n2ZmZmqXr26AgICrtvnWhISEuTv72/fQkNDy7l6AADgKlw27BSx2WwOjy3LKtb2c7/UJzY2Vjk5Ofbt+PHj5VIrAABwPS4bdoKCgiSp2AxNVlaWfbYnKChIly9fVnZ29nX7XIunp6f8/PwcNgAAYCaXDTvh4eEKCgpSSkqKve3y5ctKTU1VVFSUJKl9+/by8PBw6HPq1CkdOHDA3gcAANzanPpprHPnzumrr76yPz569Kj27dun2rVrq1GjRpowYYLi4+MVERGhiIgIxcfHy9vbW0OHDpUk+fv7a9SoUXr66adVp04d1a5dW5MmTVKrVq3sn84CAAC3NqeGnT179qh79+72xxMnTpQkjRgxQsuWLdOUKVN08eJFjRkzRtnZ2ercubM2bNggX19f+3P++te/yt3dXQ899JAuXryoe++9V8uWLZObm1ulvx4AAOB6bJZlWc4uwtlyc3Pl7++vnJwcp6/fCZu69qb7ZrzQrwIrAQDAtd3s+7fLrtkBAAAoD4QdAABgNMIOAAAwmlMXKN8KSrIGBwAAlD9mdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0biDchXGN6QDAPDLmNkBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Vw67Fy5ckXPPvuswsPD5eXlpSZNmmj27NkqLCy097EsS3FxcQoJCZGXl5e6deum9PR0J1YNAABciUuHnTlz5ui1117TggULdOjQIc2dO1cvvvii5s+fb+8zd+5czZs3TwsWLNDu3bsVFBSknj17Ki8vz4mVAwAAV+HSYWfXrl26//771a9fP4WFhWnQoEGKiYnRnj17JP00q5OYmKhp06Zp4MCBatmypZYvX64LFy4oOTnZydUDAABX4NJhp2vXrtq0aZO++OILSdJ//vMf7dixQ3379pUkHT16VJmZmYqJibE/x9PTU9HR0dq5c+d1j5ufn6/c3FyHDQAAmMnd2QXcyDPPPKOcnBzdcccdcnNzU0FBgZ5//nkNGTJEkpSZmSlJCgwMdHheYGCgjh07dt3jJiQkaNasWRVXOAAAcBkuPbPz9ttv66233lJycrI+++wzLV++XH/5y1+0fPlyh342m83hsWVZxdquFhsbq5ycHPt2/PjxCqkfAAA4n0vP7EyePFlTp07V7373O0lSq1atdOzYMSUkJGjEiBEKCgqS9NMMT3BwsP15WVlZxWZ7rubp6SlPT8+KLR4AALgEl57ZuXDhgqpVcyzRzc3N/tHz8PBwBQUFKSUlxb7/8uXLSk1NVVRUVKXWCgAAXJNLz+z0799fzz//vBo1aqQWLVpo7969mjdvnh555BFJP12+mjBhguLj4xUREaGIiAjFx8fL29tbQ4cOdXL1AADAFbh02Jk/f76mT5+uMWPGKCsrSyEhIRo9erRmzJhh7zNlyhRdvHhRY8aMUXZ2tjp37qwNGzbI19fXiZUDAABXYbMsy3J2Ec6Wm5srf39/5eTkyM/Pr1yPHTZ1bbker7QyXujn7BIAAChXN/v+7dJrdgAAAMqKsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADCau7MLQOUIm7q2RP0zXuhXQZUAAFC5mNkBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjObu7AJQ9YVNXXvTfTNe6FeBlQAAUBwzOwAAwGiEHQAAYDTCDgAAMBprdlCpWN8DAKhszOwAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGguH3ZOnjyp3//+96pTp468vb3Vtm1bpaWl2fdblqW4uDiFhITIy8tL3bp1U3p6uhMrBgAArsSlvwg0OztbXbp0Uffu3fXRRx+pfv36OnLkiGrVqmXvM3fuXM2bN0/Lli1T06ZN9dxzz6lnz546fPiwfH19nVd8FVeSL+wEAMCVuXTYmTNnjkJDQ5WUlGRvCwsLs//bsiwlJiZq2rRpGjhwoCRp+fLlCgwMVHJyskaPHl3ZJQMAABfj0pex1qxZow4dOui3v/2t6tevr3bt2mnJkiX2/UePHlVmZqZiYmLsbZ6enoqOjtbOnTudUTIAAHAxLh12vv76ay1atEgRERFav369Hn/8cY0fP15vvPGGJCkzM1OSFBgY6PC8wMBA+75ryc/PV25ursMGAADM5NKXsQoLC9WhQwfFx8dLktq1a6f09HQtWrRIf/jDH+z9bDabw/MsyyrWdrWEhATNmjWrYooGAAAuxaVndoKDg9W8eXOHtsjISH3zzTeSpKCgIEkqNouTlZVVbLbnarGxscrJybFvx48fL+fKAQCAqyhV2GnSpInOnDlTrP3s2bNq0qRJmYsq0qVLFx0+fNih7YsvvlDjxo0lSeHh4QoKClJKSop9/+XLl5WamqqoqKjrHtfT01N+fn4OGwAAMFOpLmNlZGSooKCgWHt+fr5OnjxZ5qKKPPXUU4qKilJ8fLweeugh/fvf/9bixYu1ePFiST9dvpowYYLi4+MVERGhiIgIxcfHy9vbW0OHDi23OgAAQNVVorCzZs0a+7/Xr18vf39/++OCggJt2rTJ4aPhZdWxY0e99957io2N1ezZsxUeHq7ExEQNGzbM3mfKlCm6ePGixowZo+zsbHXu3FkbNmzgHjsAAECSZLMsy7rZztWq/XTVy2az6edP8/DwUFhYmF566SXdd9995VtlBcvNzZW/v79ycnLK/ZIWN+crvYwX+jm7BACAC7vZ9+8SzewUFhZK+mmtzO7du1W3bt2yVQkAAFDBSrVm5+jRo+VdBwAAQIUo9X12Nm3apE2bNikrK8s+41Pk9ddfL3NhAAAA5aFUYWfWrFmaPXu2OnTooODg4BvewA8AAMCZShV2XnvtNS1btkzDhw8v73oAAADKValuKnj58uUb3rQPAADAVZQq7Dz66KNKTk4u71oAAADKXakuY126dEmLFy/Wxo0b1bp1a3l4eDjsnzdvXrkUBwAAUFalCjuff/652rZtK0k6cOCAwz4WKwMAAFdSqrCzZcuW8q4DAACgQpRqzQ4AAEBVUaqZne7du9/wctXmzZtLXRAAAEB5KlXYKVqvU+THH3/Uvn37dODAAY0YMaI86gIAACgXpQo7f/3rX6/ZHhcXp3PnzpWpIAAAgPJUrmt2fv/73/O9WAAAwKWUa9jZtWuXatSoUZ6HBAAAKJNSXcYaOHCgw2PLsnTq1Cnt2bNH06dPL5fCAAAAykOpwo6/v7/D42rVqqlZs2aaPXu2YmJiyqUwAACA8lCqsJOUlFTedQDFhE1dW6L+GS/0q6BKAABVWanCTpG0tDQdOnRINptNzZs3V7t27cqrLgAAgHJRqrCTlZWl3/3ud9q6datq1aoly7KUk5Oj7t27a9WqVapXr1551wkAAFAqpfo01rhx45Sbm6v09HT98MMPys7O1oEDB5Sbm6vx48eXd40AAAClVqqZnY8//lgbN25UZGSkva158+Z69dVXWaAMAABcSqlmdgoLC+Xh4VGs3cPDQ4WFhWUuCgAAoLyUKuzcc889evLJJ/Xtt9/a206ePKmnnnpK9957b7kVBwAAUFalCjsLFixQXl6ewsLCdNttt+n2229XeHi48vLyNH/+/PKuEQAAoNRKtWYnNDRUn332mVJSUvTf//5XlmWpefPm6tGjR3nXBwAAUCYlmtnZvHmzmjdvrtzcXElSz549NW7cOI0fP14dO3ZUixYttH379gopFAAAoDRKFHYSExP12GOPyc/Pr9g+f39/jR49WvPmzSu34gAAAMqqRGHnP//5j3r37n3d/TExMUpLSytzUQAAAOWlRGHnu+++u+ZHzou4u7vr+++/L3NRAAAA5aVEYadBgwbav3//dfd//vnnCg4OLnNRAAAA5aVEYadv376aMWOGLl26VGzfxYsXNXPmTN13333lVhwAAEBZleij588++6xWr16tpk2b6oknnlCzZs1ks9l06NAhvfrqqyooKNC0adMqqlYAAIASK1HYCQwM1M6dO/WnP/1JsbGxsixLkmSz2dSrVy8tXLhQgYGBFVIoAABAaZT4poKNGzfWunXrlJ2dra+++kqWZSkiIkIBAQEVUR8AAECZlOoOypIUEBCgjh07lmctAAAA5a5U340FAABQVRB2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIzm7uwCAGcIm7r2pvtmvNCvAisBAFQ0ZnYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABitSoWdhIQE2Ww2TZgwwd5mWZbi4uIUEhIiLy8vdevWTenp6c4rEgAAuJQqE3Z2796txYsXq3Xr1g7tc+fO1bx587RgwQLt3r1bQUFB6tmzp/Ly8pxUKQAAcCVVIuycO3dOw4YN05IlSxQQEGBvtyxLiYmJmjZtmgYOHKiWLVtq+fLlunDhgpKTk51YMQAAcBVVIuyMHTtW/fr1U48ePRzajx49qszMTMXExNjbPD09FR0drZ07d173ePn5+crNzXXYAACAmVz+DsqrVq3SZ599pt27dxfbl5mZKUkKDAx0aA8MDNSxY8eue8yEhATNmjWrfAsFAAAuyaVndo4fP64nn3xSb731lmrUqHHdfjabzeGxZVnF2q4WGxurnJwc+3b8+PFyqxkAALgWl57ZSUtLU1ZWltq3b29vKygo0LZt27RgwQIdPnxY0k8zPMHBwfY+WVlZxWZ7rubp6SlPT8+KKxwAALgMl57Zuffee7V//37t27fPvnXo0EHDhg3Tvn371KRJEwUFBSklJcX+nMuXLys1NVVRUVFOrBwAALgKl57Z8fX1VcuWLR3afHx8VKdOHXv7hAkTFB8fr4iICEVERCg+Pl7e3t4aOnSoM0oGAAAuxqXDzs2YMmWKLl68qDFjxig7O1udO3fWhg0b5Ovr6+zSAACAC6hyYWfr1q0Oj202m+Li4hQXF+eUegAAgGtz6TU7AAAAZUXYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIzm7uwCgPISNnWt04+b8UK/CqkBAFB6zOwAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBrfeg6Uo5J+8zrfkg4AFY+ZHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNO6gDFQRJbk7M3dmBoD/YWYHAAAYjbADAACMxmUswIlK+sWhAICSY2YHAAAYjbADAACMRtgBAABGY80OAD7WDsBozOwAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEZz6bCTkJCgjh07ytfXV/Xr19cDDzygw4cPO/SxLEtxcXEKCQmRl5eXunXrpvT0dCdVDAAAXI1L30E5NTVVY8eOVceOHXXlyhVNmzZNMTExOnjwoHx8fCRJc+fO1bx587Rs2TI1bdpUzz33nHr27KnDhw/L19fXya8AuLVxZ2YArsClw87HH3/s8DgpKUn169dXWlqafv3rX8uyLCUmJmratGkaOHCgJGn58uUKDAxUcnKyRo8e7YyyAQCAC3Hpy1g/l5OTI0mqXbu2JOno0aPKzMxUTEyMvY+np6eio6O1c+fO6x4nPz9fubm5DhsAADBTlQk7lmVp4sSJ6tq1q1q2bClJyszMlCQFBgY69A0MDLTvu5aEhAT5+/vbt9DQ0IorHAAAOFWVCTtPPPGEPv/8c61cubLYPpvN5vDYsqxibVeLjY1VTk6OfTt+/Hi51wsAAFyDS6/ZKTJu3DitWbNG27ZtU8OGDe3tQUFBkn6a4QkODra3Z2VlFZvtuZqnp6c8PT0rrmAAAOAyXHpmx7IsPfHEE1q9erU2b96s8PBwh/3h4eEKCgpSSkqKve3y5ctKTU1VVFRUZZcLAABckEvP7IwdO1bJycn64IMP5Ovra1+H4+/vLy8vL9lsNk2YMEHx8fGKiIhQRESE4uPj5e3traFDhzq5egAA4ApcOuwsWrRIktStWzeH9qSkJI0cOVKSNGXKFF28eFFjxoxRdna2OnfurA0bNnCPHQAAIMnFw45lWb/Yx2azKS4uTnFxcRVfEAAAqHJcOuwAcD0luSsyALgCl16gDAAAUFaEHQAAYDTCDgAAMBprdgADsa4GAP6HmR0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTuoAwAVynJ3aczXuhXgZUAKC/M7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2PngOokkryEXHTlXQs+Mg8bjXM7AAAAKMRdgAAgNG4jAXAJVTFy1JcPgKqBmZ2AACA0Qg7AADAaIQdAABgNNbsAEAl4RvVAedgZgcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI07KAPALYY7OeNWw8wOAAAwGmEHAAAYjctYAACjcJkOP8fMDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0fjoOQCgXPCRb7gqZnYAAIDRCDsAAMBoXMYCABdUkktCFclV6qgoXHq7NTCzAwAAjEbYAQAARiPsAAAAo7FmBwBQ6UxfCyS5zmusqLVGVWm9EzM7AADAaMaEnYULFyo8PFw1atRQ+/bttX37dmeXBAAAXIARl7HefvttTZgwQQsXLlSXLl30t7/9TX369NHBgwfVqFEjZ5cHAIDTVKXLTRXFiJmdefPmadSoUXr00UcVGRmpxMREhYaGatGiRc4uDQAAOFmVDzuXL19WWlqaYmJiHNpjYmK0c+dOJ1UFAABcRZW/jHX69GkVFBQoMDDQoT0wMFCZmZnXfE5+fr7y8/Ptj3NyciRJubm55V5fYf6Fcj8mAKDylfQ9oir+/i/JayzJ66uI99erj2tZ1g37VfmwU8Rmszk8tiyrWFuRhIQEzZo1q1h7aGhohdQGAKj6/BOdXUHFq6jXWNFjl5eXJ39//+vur/Jhp27dunJzcys2i5OVlVVstqdIbGysJk6caH9cWFioH374QXXq1LluQCqN3NxchYaG6vjx4/Lz8yu34+LaGO/Kw1hXHsa68jDWlae8xtqyLOXl5SkkJOSG/ap82Klevbrat2+vlJQU/eY3v7G3p6Sk6P7777/mczw9PeXp6enQVqtWrQqr0c/Pj/84lYjxrjyMdeVhrCsPY115ymOsbzSjU6TKhx1JmjhxooYPH64OHTrorrvu0uLFi/XNN9/o8ccfd3ZpAADAyYwIO4MHD9aZM2c0e/ZsnTp1Si1bttS6devUuHFjZ5cGAACczIiwI0ljxozRmDFjnF2GA09PT82cObPYJTNUDMa78jDWlYexrjyMdeWp7LG2Wb/0eS0AAIAqrMrfVBAAAOBGCDsAAMBohB0AAGA0wg4AADAaYaeMFi5cqPDwcNWoUUPt27fX9u3bb9g/NTVV7du3V40aNdSkSRO99tprlVRp1VeSsV69erV69uypevXqyc/PT3fddZfWr19fidVWfSX92S7yySefyN3dXW3btq3YAg1S0rHOz8/XtGnT1LhxY3l6euq2227T66+/XknVVm0lHesVK1aoTZs28vb2VnBwsB5++GGdOXOmkqqturZt26b+/fsrJCRENptN77///i8+p0LfHy2U2qpVqywPDw9ryZIl1sGDB60nn3zS8vHxsY4dO3bN/l9//bXl7e1tPfnkk9bBgwetJUuWWB4eHtY777xTyZVXPSUd6yeffNKaM2eO9e9//9v64osvrNjYWMvDw8P67LPPKrnyqqmk413k7NmzVpMmTayYmBirTZs2lVNsFVeasR4wYIDVuXNnKyUlxTp69Kj1r3/9y/rkk08qseqqqaRjvX37dqtatWrWyy+/bH399dfW9u3brRYtWlgPPPBAJVde9axbt86aNm2a9e6771qSrPfee++G/Sv6/ZGwUwadOnWyHn/8cYe2O+64w5o6deo1+0+ZMsW64447HNpGjx5t/epXv6qwGk1R0rG+lubNm1uzZs0q79KMVNrxHjx4sPXss89aM2fOJOzcpJKO9UcffWT5+/tbZ86cqYzyjFLSsX7xxRetJk2aOLS98sorVsOGDSusRhPdTNip6PdHLmOV0uXLl5WWlqaYmBiH9piYGO3cufOaz9m1a1ex/r169dKePXv0448/VlitVV1pxvrnCgsLlZeXp9q1a1dEiUYp7XgnJSXpyJEjmjlzZkWXaIzSjPWaNWvUoUMHzZ07Vw0aNFDTpk01adIkXbx4sTJKrrJKM9ZRUVE6ceKE1q1bJ8uy9N133+mdd95Rv379KqPkW0pFvz8acwflynb69GkVFBQU+2b1wMDAYt/AXiQzM/Oa/a9cuaLTp08rODi4wuqtykoz1j/30ksv6fz583rooYcqokSjlGa8v/zyS02dOlXbt2+Xuzu/Vm5Wacb666+/1o4dO1SjRg299957On36tMaMGaMffviBdTs3UJqxjoqK0ooVKzR48GBdunRJV65c0YABAzR//vzKKPmWUtHvj8zslJHNZnN4bFlWsbZf6n+tdhRX0rEusnLlSsXFxentt99W/fr1K6o849zseBcUFGjo0KGaNWuWmjZtWlnlGaUkP9uFhYWy2WxasWKFOnXqpL59+2revHlatmwZszs3oSRjffDgQY0fP14zZsxQWlqaPv74Yx09epQvma4gFfn+yJ9gpVS3bl25ubkV+4sgKyurWDotEhQUdM3+7u7uqlOnToXVWtWVZqyLvP322xo1apT+8Y9/qEePHhVZpjFKOt55eXnas2eP9u7dqyeeeELST2/IlmXJ3d1dGzZs0D333FMptVc1pfnZDg4OVoMGDeTv729vi4yMlGVZOnHihCIiIiq05qqqNGOdkJCgLl26aPLkyZKk1q1by8fHR3fffbeee+45ZuPLUUW/PzKzU0rVq1dX+/btlZKS4tCekpKiqKioaz7nrrvuKtZ/w4YN6tChgzw8PCqs1qquNGMt/TSjM3LkSCUnJ3ONvQRKOt5+fn7av3+/9u3bZ98ef/xxNWvWTPv27VPnzp0rq/QqpzQ/2126dNG3336rc+fO2du++OILVatWTQ0bNqzQequy0oz1hQsXVK2a49ukm5ubpP/NOqB8VPj7Y7ksc75FFX2McenSpdbBgwetCRMmWD4+PlZGRoZlWZY1depUa/jw4fb+RR+te+qpp6yDBw9aS5cu5aPnN6mkY52cnGy5u7tbr776qnXq1Cn7dvbsWWe9hCqlpOP9c3wa6+aVdKzz8vKshg0bWoMGDbLS09Ot1NRUKyIiwnr00Ued9RKqjJKOdVJSkuXu7m4tXLjQOnLkiLVjxw6rQ4cOVqdOnZz1EqqMvLw8a+/evdbevXstSda8efOsvXv32j/mX9nvj4SdMnr11Vetxo0bW9WrV7fuvPNOKzU11b5vxIgRVnR0tEP/rVu3Wu3atbOqV69uhYWFWYsWLarkiquukox1dHS0JanYNmLEiMovvIoq6c/21Qg7JVPSsT506JDVo0cPy8vLy2rYsKE1ceJE68KFC5VcddVU0rF+5ZVXrObNm1teXl5WcHCwNWzYMOvEiROVXHXVs2XLlhv+Dq7s90ebZTEXBwAAzMWaHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADgpm3btk39+/dXSEiIbDab3n///Qo9X1xcnGw2m8MWFBRUomMQdgBIUrFfWv/973/1q1/9SjVq1FDbtm1LfdyMjAzZbDbt27evzDWWxuHDhxUUFKS8vDynnL+8hIWFKTEx8ab6duzYUatXr67YgnDLOn/+vNq0aaMFCxZU2jlbtGihU6dO2bf9+/eX6Pl86zlwixg5cqTOnj173b/CTp06pYCAAPvjmTNnysfHR4cPH1bNmjVLfY7Q0FCdOnVKdevWLUv5pTZt2jSNHTtWvr6+Tjm/M0yfPl2TJk3SAw88UOyLLIGy6tOnj/r06XPd/ZcvX9azzz6rFStW6OzZs2rZsqXmzJmjbt26lfqc7u7uJZ7NuRr/CwBIkoKCguTp6Wl/fOTIEXXt2lWNGzdWnTp1Sn1cNzc3BQUFyd298v+2OnHihNasWaOHH3640s/tTP369VNOTo7Wr1/v7FJwC3r44Yf1ySefaNWqVfr888/129/+Vr1799aXX35Z6mN++eWXCgkJUXh4uH73u9/p66+/LtHzCTsAJDlexrLZbEpLS9Ps2bNls9kUFxcnSTp58qQGDx6sgIAA1alTR/fff78yMjIk/XRdffny5frggw/s19W3bt1a7DLW1q1bZbPZtH79erVr105eXl665557lJWVpY8++kiRkZHy8/PTkCFDdOHCBXt9lmVp7ty5atKkiby8vNSmTRu98847N3xN//d//6c2bdqoYcOG9rZjx46pf//+CggIkI+Pj1q0aKF169bZ9x88eFB9+/ZVzZo1FRgYqOHDh+v06dP2/YWFhZozZ45uv/12eXp6qlGjRnr++eft+/fv36977rlHXl5eqlOnjv74xz/q3Llz9v0jR47UAw88oL/85S8KDg5WnTp1NHbsWP3444/2PllZWerfv7+8vLwUHh6uFStWFHttcXFxatSokTw9PRUSEqLx48fb97m5ualv375auXLlDccHKG9HjhzRypUr9Y9//EN33323brvtNk2aNEldu3ZVUlJSqY7ZuXNnvfHGG1q/fr2WLFmizMxMRUVF6cyZMzd9DMIOgGJOnTqlFi1a6Omnn9apU6c0adIkXbhwQd27d1fNmjW1bds27dixQzVr1lTv3r11+fJlTZo0SQ899JB69+5tv64eFRV13XPExcVpwYIF2rlzp44fP66HHnpIiYmJSk5O1tq1a5WSkqL58+fb+z/77LNKSkrSokWLlJ6erqeeekq///3vlZqaet1zbNu2TR06dHBoGzt2rPLz87Vt2zbt379fc+bMsV+mO3XqlKKjo9W2bVvt2bNHH3/8sb777js99NBD9ufHxsZqzpw5mj59ug4ePKjk5GQFBgZKki5cuKDevXsrICBAu3fv1j/+8Q9t3LhRTzzxhEMNW7Zs0ZEjR7RlyxYtX75cy5Yt07Jly+z7R44cqYyMDG3evFnvvPOOFi5cqKysLPv+d955R3/961/1t7/9TV9++aXef/99tWrVyuEcnTp10vbt2687NkBF+Oyzz2RZlpo2baqaNWvat9TUVB05ckTS/9bx3Wi7+v9Mnz599OCDD6pVq1bq0aOH1q5dK0lavnz5TdfFmh0AxRRddqpZs6b9Ovnrr7+uatWq6e9//7tsNpskKSkpSbVq1dLWrVsVExMjLy8v5efn39S19eeee05dunSRJI0aNUqxsbE6cuSImjRpIkkaNGiQtmzZomeeeUbnz5/XvHnztHnzZt11112SpCZNmmjHjh3629/+pujo6GueIyMjQ+3bt3do++abb+y/OIuOU2TRokW68847FR8fb297/fXXFRoaqi+++ELBwcF6+eWXtWDBAo0YMUKSdNttt6lr166SpBUrVujixYt644035OPjI0lasGCB+vfvrzlz5thDUUBAgBYsWCA3Nzfdcccd6tevnzZt2qTHHntMX3zxhT766CN9+umn6ty5syRp6dKlioyMdHgNQUFB6tGjhzw8PNSoUSN16tTJ4XU2aNBA33zzjQoLC1m3g0pTWFgoNzc3paWlyc3NzWFf0R8VDRo00KFDh254nKvXD/6cj4+PWrVqVaLLYoQdADclLS1NX331VbGFvpcuXbL/xVYSrVu3tv87MDBQ3t7eDsEjMDBQ//73vyX9dGnp0qVL6tmzp8MxLl++rHbt2l33HBcvXlSNGjUc2saPH68//elP2rBhg3r06KEHH3zQXktaWpq2bNlyzQXZR44c0dmzZ5Wfn6977733muc7dOiQ2rRpYw86ktSlSxcVFhbq8OHD9rDTokULhzeC4OBg+6dLDh06JHd3d4cZqTvuuEO1atWyP/7tb3+rxMRENWnSRL1791bfvn3Vv39/h3VRXl5eKiwsVH5+vry8vK47RkB5ateunQoKCpSVlaW77777mn08PDx0xx13lPoc+fn5OnTo0HWPfy2EHQA3pbCwUO3bt7/m+pF69eqV+HgeHh72f9tsNofHRW2FhYX2c0vS2rVr1aBBA4d+Vy+q/rm6desqOzvboe3RRx9Vr169tHbtWm3YsEEJCQl66aWXNG7cOBUWFtpnYX4uODj4FxdFWpZln/X6uavbb/RaLcsq1v/nQkNDdfjwYaWkpGjjxo0aM2aMXnzxRaWmptqP/cMPP8jb25ugg3J37tw5ffXVV/bHR48e1b59+1S7dm01bdpUw4YN0x/+8Ae99NJLateunU6fPq3NmzerVatW6tu3b4nPN2nSJPXv31+NGjVSVlaWnnvuOeXm5tpnV28Gc5sAbsqdd96pL7/8UvXr19ftt9/usPn7+0uSqlevroKCgnI/d/PmzeXp6alvvvmm2LlDQ0Ov+7x27drp4MGDxdpDQ0P1+OOPa/Xq1Xr66ae1ZMkS+2tMT09XWFhYsfP4+PgoIiJCXl5e2rRp03Xr3Ldvn86fP29v++STT1StWjU1bdr0pl5rZGSkrly5oj179tjbDh8+rLNnzzr08/Ly0oABA/TKK69o69at2rVrl8O9Rw4cOKA777zzps4JlMSePXvUrl07+6zqxIkT1a5dO82YMUPST5e3//CHP+jpp59Ws2bNNGDAAP3rX/+64f/VGzlx4oSGDBmiZs2aaeDAgapevbo+/fRTNW7c+KaPwcwOcAvJyckpdnO/2rVrq1GjRr/43GHDhunFF1/U/fffr9mzZ6thw4b65ptvtHr1ak2ePFkNGzZUWFiY1q9fr8OHD6tOnTr2EFRWvr6+mjRpkp566ikVFhaqa9euys3N1c6dO1WzZs3r/oXXq1cvPfrooyooKLBfNpowYYL69Omjpk2bKjs7W5s3b7avhxk7dqyWLFmiIUOGaPLkyapbt66++uorrVq1SkuWLFGNGjX0zDPPaMqUKapevbq6dOmi77//Xunp6Ro1apSGDRummTNnasSIEYqLi9P333+vcePGafjw4fZLWL+kWbNm6t27tx577DEtXrxY7u7umjBhgsMMzbJly1RQUKDOnTvL29tbb775pry8vBx++W/fvl0xMTGlHXLgurp162afgbwWDw8PzZo1S7NmzSqX861atarMx2BmB7iFbN261f4XWdFW9NfYL/H29ta2bdvUqFEjDRw4UJGRkXrkkUd08eJF+fn5SZIee+wxNWvWTB06dFC9evX0ySeflFvtf/7znzVjxgwlJCQoMjJSvXr10ocffqjw8PDrPqdv377y8PDQxo0b7W0FBQUaO3asIiMj1bt3bzVr1kwLFy6UJIWEhOiTTz5RQUGBevXqpZYtW+rJJ5+Uv7+/fZHv9OnT9fTTT2vGjBmKjIzU4MGD7Z+U8vb21vr16/XDDz+oY8eOGjRokO69994S32k2KSlJoaGhio6O1sCBA/XHP/5R9evXt++vVauWlixZoi5duqh169batGmTPvzwQ/v9kE6ePKmdO3fecvcXAq7HZt0ongFAFbdw4UJ98MEHt9QN9iZPnqycnBwtXrzY2aUALoHLWACM9sc//lHZ2dnKy8u7Zb4yon79+po0aZKzywBcBjM7AADAaKzZAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABG+/92DNCakkpFkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(decay_sample, bins=40)\n",
    "#plt.savefig('outputHistogram.pdf')\n",
    "plt.title(\"Muon lifetimes\")\n",
    "plt.xlabel(\"Lifetime (seconds)\")\n",
    "plt.ylabel(\"Count\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf7b9994",
   "metadata": {},
   "source": [
    "<b>2. Estimate the lifetime of the muon<b/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "8a1673a6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean muon lifetime from dataset = 2.079653901251733e-06 seconds\n"
     ]
    }
   ],
   "source": [
    "#Estimate of the muon lifetime from this dataset\n",
    "#Taking the mean of all the datapoints\n",
    "mean_lifetime = decay_sample.mean()\n",
    "print(\"Mean muon lifetime from dataset =\",mean_lifetime,\"seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7c47f4f",
   "metadata": {},
   "source": [
    "Mean muon lifetime found through this dataset is lower than the actual lifetime of the muon which is 2.2e-6 microseconds."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d777471",
   "metadata": {},
   "source": [
    "The value that we have obtained is lower than the actual value because we have cut off the tail of the exponential distribution at 10 microseconds. The exponential distribution can otherwise produce decay times between 0ms and infinity. In this way we are disregarding rare decay times that happen to be longer than 10 microseconds. As a result of the missing longer lifetimes, we get an lower average value."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7abed8e2",
   "metadata": {},
   "source": [
    "## Repeat 500 times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d362420",
   "metadata": {},
   "source": [
    "<b>3. Repeat the experiment 500 times and obtain distribution of mean lifetimes<b/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "aa3b0345",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function for the experiment\n",
    "def experiment(lifetime,sample_size):\n",
    "    decay_times = np.random.exponential(scale=lifetime,size = int(1.5*sample_size))\n",
    "    decays_cut = np.where(decay_times<10e-6)[0]\n",
    "    decay_sample = decay_times[decays_cut[:sample_size]]\n",
    "    return decay_sample.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "0a08d7b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Repeat experiment 500 times\n",
    "mean_lifetimes = np.array([])\n",
    "for i in range(0,500):\n",
    "    mean_lifetime = experiment(2.2e-6,1000)\n",
    "    mean_lifetimes = np.append(mean_lifetimes,mean_lifetime)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "id": "2863c997",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAHFCAYAAAAHcXhbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4oElEQVR4nO3de1hVZf7//9cWcAsKqCgnRdDygHmiNM+DlqKojIdKyzRwyg6iZmaamYlOidpHsxkPHaY0+3pq8pCTppIKVmopo2lq5gHSCqVMAUkxZf3+6MeetoACAnsvfD6ua19X6173Xut9cxu8uNdiL4thGIYAAABMqpKjCwAAALgZhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBmgHCxevFgWi0UWi0WJiYn59huGodtvv10Wi0VdunQp9/puZampqbJYLFq8eLGtLW++UlNT7fq++OKLqlevnlxdXVW9evVinWfDhg2Ki4srcF9ISIhiYmKKdTwA/0OYAcqRp6en3nnnnXztSUlJOn78uDw9PR1QFa7Vu3dv7dy5UwEBAba2jz76SK+88ooeeeQRJSUl6dNPPy3WMTds2KCpU6cWuG/NmjWaPHnyTdUM3MpcHV0AcCsZNGiQli5dqvnz58vLy8vW/s4776h9+/bKzMx0YHXIU7t2bdWuXduu7ZtvvpEkjR49Wr6+vqV6vrCwsFI9HnCrYWUGKEcPPfSQJGn58uW2toyMDK1atUp/+9vfCnzP5cuX9fLLL6tJkyayWq2qXbu2hg0bpp9//tmu38qVKxUREaGAgAC5u7srNDRUzz//vLKzs+36xcTEqFq1ajp27Jh69eqlatWqKSgoSM8++6xycnJuOIaQkBD16dNHH3/8scLCwmzn+vjjjyX9cYkmNDRUVatW1d133609e/bYvb9Lly4FXkqLiYlRSEiIXduvv/6qESNGqE6dOqpcubIaNGigSZMm5avTYrFo5MiRev/99xUaGioPDw+1bNnSVlNxXXuZKSQkRC+++KIkyc/PTxaLxe6S0cqVK9W+fXtVrVpV1apVU48ePbR37167sc2fP99Wa97rz8f/82WmxMREWSwWLVu2TBMmTFBAQICqVaumqKgonTlzRllZWXr88cdVq1Yt1apVS8OGDdOFCxfsxmAYhhYsWKBWrVrJ3d1dNWrU0P33368TJ07Y9du7d6/69OkjX19fWa1WBQYGqnfv3vrhhx9K9LUDHIEwA5QjLy8v3X///Xr33XdtbcuXL1elSpU0aNCgfP1zc3PVt29fzZgxQ4MHD9b69es1Y8YMJSQkqEuXLrp48aKt79GjR9WrVy+988472rhxo8aMGaMPPvhAUVFR+Y77+++/669//avuvfdeffTRR/rb3/6m1157TTNnzizSOL7++mtNnDhREyZM0OrVq+Xt7a0BAwZoypQp+te//qXp06dr6dKlysjIUJ8+fezqLKpLly6pa9euWrJkicaOHav169dryJAhmjVrlgYMGJCv//r16zVv3jxNmzZNq1atUs2aNdW/f/98P7xLYs2aNXr00UclSRs3btTOnTv12GOPSZKmT5+uhx56SE2bNtUHH3yg999/X1lZWercubMOHTokSZo8ebLuv/9+SdLOnTttrz9fxirICy+8oPT0dC1evFizZ89WYmKiHnroId13333y9vbW8uXLNX78eL3//vt64YUX7N77xBNPaMyYMerWrZvWrl2rBQsW6ODBg+rQoYPOnDkjScrOzlb37t115swZzZ8/XwkJCZo7d67q1aunrKysm/66AeXGAFDmFi1aZEgydu/ebWzbts2QZHzzzTeGYRhGmzZtjJiYGMMwDOOOO+4wwsPDbe9bvny5IclYtWqV3fF2795tSDIWLFhQ4Plyc3ON33//3UhKSjIkGV9//bVtX3R0tCHJ+OCDD+ze06tXL6Nx48Y3HEtwcLDh7u5u/PDDD7a2ffv2GZKMgIAAIzs729a+du1aQ5Kxbt06W1t4eLjdGP9cV3BwsG37jTfeKLDOmTNnGpKMzZs329okGX5+fkZmZqat7fTp00alSpWM+Pj4644nJSXFkGQsWrTI1pY3XykpKba2KVOmGJKMn3/+2dZ28uRJw9XV1Rg1apTdMbOysgx/f39j4MCBtrbY2FijsG+5wcHBRnR0tG07799IVFSUXb8xY8YYkozRo0fbtffr18+oWbOmbXvnzp2GJGP27Nl2/U6dOmW4u7sb48ePNwzDMPbs2WNIMtauXVtgXYBZsDIDlLPw8HDddtttevfdd3XgwAHt3r270EtMH3/8sapXr66oqChduXLF9mrVqpX8/f3t/jLqxIkTGjx4sPz9/eXi4iI3NzeFh4dLkg4fPmx3XIvFkm/FpkWLFvr++++LNIZWrVqpTp06tu3Q0FBJf1xC8vDwyNde1OP+2datW1W1alXbikaevMsxW7ZssWvv2rWr3Q3Ufn5+8vX1LdG5i2rTpk26cuWKHnnkEbv5qVKlisLDwwv8y7Xi6NOnj9123tezd+/e+dp//fVX26Wmjz/+WBaLRUOGDLGry9/fXy1btrTVdfvtt6tGjRqaMGGC3njjDdtKEmA23AAMlDOLxaJhw4bpH//4hy5duqRGjRqpc+fOBfY9c+aMzp8/r8qVKxe4/5dffpEkXbhwQZ07d1aVKlX08ssvq1GjRvLw8NCpU6c0YMCAfJd5PDw8VKVKFbs2q9WqS5cuFWkMNWvWtNvOq6+w9qIe98/Onj0rf39/WSwWu3ZfX1+5urrq7Nmzdu0+Pj75jmG1Wkt0iauo8i7XtGnTpsD9lSrd3O+LJfk6V6tWTWfOnJFhGPLz8yvwuA0aNJAkeXt7KykpSa+88opeeOEFnTt3TgEBARo+fLhefPFFubm53VT9QHkhzAAOEBMTo5deeklvvPGGXnnllUL71apVSz4+Ptq4cWOB+/NWIrZu3aqffvpJiYmJttUYSTp//nyp1l0aqlSpooyMjHztecEsj4+Pj7788ksZhmEXaNLT03XlyhXVqlWrzGu9kbwaPvzwQwUHBzu4mv+pVauWLBaLPvvsM1mt1nz7/9zWvHlzrVixQoZhaP/+/Vq8eLGmTZsmd3d3Pf/88+VZNlBihBnAAerUqaPnnntO3377raKjowvt16dPH61YsUJXr15V27ZtC+2X98P+2h9cb775ZukUXIpCQkL073//Wzk5ObZ6z549qx07dtj9ufq9996rDz74QGvXrlX//v1t7UuWLLHtd7QePXrI1dVVx48f13333XfdvnljvXjxotzd3cu0rj59+mjGjBn68ccfNXDgwCK9x2KxqGXLlnrttde0ePFi/fe//y3TGoHSRJgBHGTGjBk37PPggw9q6dKl6tWrl55++mndfffdcnNz0w8//KBt27apb9++6t+/vzp06KAaNWroySef1JQpU+Tm5qalS5fq66+/LoeRFM/QoUP15ptvasiQIRo+fLjOnj2rWbNm2QUZSXrkkUc0f/58RUdHKzU1Vc2bN9fnn3+u6dOnq1evXurWrZuDRvA/ISEhmjZtmiZNmqQTJ06oZ8+eqlGjhs6cOaOvvvpKVatWtX1QXvPmzSVJM2fOVGRkpFxcXNSiRYtCLyHejI4dO+rxxx/XsGHDtGfPHv3lL39R1apVlZaWps8//1zNmzfXU089pY8//lgLFixQv3791KBBAxmGodWrV+v8+fPq3r17qdcFlBXCDODEXFxctG7dOr3++ut6//33FR8fL1dXV9WtW1fh4eG2H5A+Pj5av369nn32WQ0ZMkRVq1ZV3759tXLlSt15550OHoW9jh076r333tOMGTPUt29fNWjQQFOmTNGGDRvsbpitUqWKtm3bpkmTJunVV1/Vzz//rDp16mjcuHGaMmWK4wZwjYkTJ6pp06Z6/fXXtXz5cuXk5Mjf319t2rTRk08+aes3ePBgffHFF1qwYIGmTZsmwzCUkpKS77N1Ssubb76pdu3a6c0339SCBQuUm5urwMBAdezYUXfffbckqWHDhqpevbpmzZqln376SZUrV1bjxo21ePHi664YAs7GYhiG4egiAAAASoo/zQYAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZW4T9nJjc3Vz/99JM8PT3zPeMFAAA4J8MwlJWVpcDAwBs+56zCh5mffvpJQUFBji4DAACUwKlTp1S3bt3r9qnwYSbvQXynTp3K93HpAADAOWVmZiooKMj2c/x6KnyYybu05OXlRZgBAMBkinKLCDcAAwAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAUyPMAAAAU3N1dAEAcDNCnl9f5L6pM3qXYSUAHIWVGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGoODTMLFy5UixYt5OXlJS8vL7Vv316ffPKJbb9hGIqLi1NgYKDc3d3VpUsXHTx40IEVAwAAZ+PQMFO3bl3NmDFDe/bs0Z49e3TPPfeob9++tsAya9YszZkzR/PmzdPu3bvl7++v7t27Kysry5FlAwAAJ+LQMBMVFaVevXqpUaNGatSokV555RVVq1ZNu3btkmEYmjt3riZNmqQBAwaoWbNmeu+99/Tbb79p2bJljiwbAAA4Eae5Z+bq1atasWKFsrOz1b59e6WkpOj06dOKiIiw9bFarQoPD9eOHTsKPU5OTo4yMzPtXgAAoOJy+OMMDhw4oPbt2+vSpUuqVq2a1qxZo6ZNm9oCi5+fn11/Pz8/ff/994UeLz4+XlOnTi3TmgGUreI8ogAAHL4y07hxY+3bt0+7du3SU089pejoaB06dMi232Kx2PU3DCNf259NnDhRGRkZttepU6fKrHYAAOB4Dl+ZqVy5sm6//XZJUuvWrbV79269/vrrmjBhgiTp9OnTCggIsPVPT0/Pt1rzZ1arVVartWyLBgAATsPhKzPXMgxDOTk5ql+/vvz9/ZWQkGDbd/nyZSUlJalDhw4OrBAAADgTh67MvPDCC4qMjFRQUJCysrK0YsUKJSYmauPGjbJYLBozZoymT5+uhg0bqmHDhpo+fbo8PDw0ePBgR5YNAACciEPDzJkzZzR06FClpaXJ29tbLVq00MaNG9W9e3dJ0vjx43Xx4kWNGDFC586dU9u2bbV582Z5eno6smwAAOBELIZhGI4uoixlZmbK29tbGRkZ8vLycnQ5AIqgrP6aKXVG7zI5LoDSV5yf3053zwwAAEBxEGYAAICpEWYAAICpEWYAAICpOfxD8wCUreLcTMsNsiVXVl/n4t4MzRziVsTKDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDWezQQAhSjuc5EcfVzgVsXKDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDUeZwCgRIrzkfypM3qXYSUAbnWszAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFMjzAAAAFNzaJiJj49XmzZt5OnpKV9fX/Xr109Hjhyx6xMTEyOLxWL3ateunYMqBgAAzsahYSYpKUmxsbHatWuXEhISdOXKFUVERCg7O9uuX8+ePZWWlmZ7bdiwwUEVAwAAZ+PqyJNv3LjRbnvRokXy9fVVcnKy/vKXv9jarVar/P39y7s8AABgAk51z0xGRoYkqWbNmnbtiYmJ8vX1VaNGjTR8+HClp6cXeoycnBxlZmbavQAAQMXlNGHGMAyNHTtWnTp1UrNmzWztkZGRWrp0qbZu3arZs2dr9+7duueee5STk1PgceLj4+Xt7W17BQUFldcQAACAAzj0MtOfjRw5Uvv379fnn39u1z5o0CDbfzdr1kytW7dWcHCw1q9frwEDBuQ7zsSJEzV27FjbdmZmJoEGAIAKzCnCzKhRo7Ru3Tpt375ddevWvW7fgIAABQcH6+jRowXut1qtslqtZVEmAABwQg4NM4ZhaNSoUVqzZo0SExNVv379G77n7NmzOnXqlAICAsqhQgAA4Owces9MbGys/t//+39atmyZPD09dfr0aZ0+fVoXL16UJF24cEHjxo3Tzp07lZqaqsTEREVFRalWrVrq37+/I0sHAABOwqErMwsXLpQkdenSxa590aJFiomJkYuLiw4cOKAlS5bo/PnzCggIUNeuXbVy5Up5eno6oGIAAOBsHH6Z6Xrc3d21adOmcqoGAACYkdP8aTYAAEBJEGYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpuTq6AAAoLyHPr3d0CQDKACszAADA1AgzAADA1AgzAADA1AgzAADA1AgzAADA1AgzAADA1AgzAADA1AgzAADA1AgzAADA1AgzAADA1HicAQDghorzKIjUGb3LsBIgP1ZmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqRFmAACAqTk0zMTHx6tNmzby9PSUr6+v+vXrpyNHjtj1MQxDcXFxCgwMlLu7u7p06aKDBw86qGIAAOBsHBpmkpKSFBsbq127dikhIUFXrlxRRESEsrOzbX1mzZqlOXPmaN68edq9e7f8/f3VvXt3ZWVlObByAADgLFwdefKNGzfabS9atEi+vr5KTk7WX/7yFxmGoblz52rSpEkaMGCAJOm9996Tn5+fli1bpieeeMIRZQMAACfiVPfMZGRkSJJq1qwpSUpJSdHp06cVERFh62O1WhUeHq4dO3Y4pEYAAOBcHLoy82eGYWjs2LHq1KmTmjVrJkk6ffq0JMnPz8+ur5+fn77//vsCj5OTk6OcnBzbdmZmZhlVDAAAnIHThJmRI0dq//79+vzzz/Pts1gsdtuGYeRryxMfH6+pU6eWSY0ASibk+fWOLgEFYF5QUTjFZaZRo0Zp3bp12rZtm+rWrWtr9/f3l/S/FZo86enp+VZr8kycOFEZGRm216lTp8qucAAA4HAODTOGYWjkyJFavXq1tm7dqvr169vtr1+/vvz9/ZWQkGBru3z5spKSktShQ4cCj2m1WuXl5WX3AgAAFZdDLzPFxsZq2bJl+uijj+Tp6WlbgfH29pa7u7ssFovGjBmj6dOnq2HDhmrYsKGmT58uDw8PDR482JGlAwAAJ+HQMLNw4UJJUpcuXezaFy1apJiYGEnS+PHjdfHiRY0YMULnzp1T27ZttXnzZnl6epZztQAAwBk5NMwYhnHDPhaLRXFxcYqLiyv7ggAAgOk4xQ3AAAAAJUWYAQAApkaYAQAApkaYAQAApkaYAQAApuY0jzMA4Hh8vD0AM2JlBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBqPMwCACoRHUuBWxMoMAAAwtRKFmQYNGujs2bP52s+fP68GDRrcdFEAAABFVaIwk5qaqqtXr+Zrz8nJ0Y8//njTRQEAABRVse6ZWbdune2/N23aJG9vb9v21atXtWXLFoWEhJRacQAAADdSrDDTr18/SZLFYlF0dLTdPjc3N4WEhGj27NmlVhwAAMCNFCvM5ObmSpLq16+v3bt3q1atWmVSFAAAQFGV6E+zU1JSSrsOAACAEinx58xs2bJFW7ZsUXp6um3FJs+7775704UBAAAURYnCzNSpUzVt2jS1bt1aAQEBslgspV0XAABAkZQozLzxxhtavHixhg4dWtr1AAAAFEuJPmfm8uXL6tChQ2nXAgAAUGwlWpl57LHHtGzZMk2ePLm06wEqjOI8Iyd1Ru8yrAQoX8V9PhT//nGzShRmLl26pLfeekuffvqpWrRoITc3N7v9c+bMKZXiAAAAbqREYWb//v1q1aqVJOmbb76x28fNwAAAoDyVKMxs27attOsAAAAokRLdAAwAAOAsSrQy07Vr1+teTtq6dWuJCwIAACiOEoWZvPtl8vz+++/at2+fvvnmm3wPoAQAAChLJQozr732WoHtcXFxunDhwk0VBAAAUByles/MkCFDeC4TAAAoV6UaZnbu3KkqVaqU5iEBAACuq0SXmQYMGGC3bRiG0tLStGfPHj4VGAAAlKsShRlvb2+77UqVKqlx48aaNm2aIiIiSqUwAACAoihRmFm0aFFp1wEAAFAiJQozeZKTk3X48GFZLBY1bdpUYWFhpVUXAABAkZQozKSnp+vBBx9UYmKiqlevLsMwlJGRoa5du2rFihWqXbt2adcJAABQoBL9NdOoUaOUmZmpgwcP6tdff9W5c+f0zTffKDMzU6NHjy7tGgEAAApVopWZjRs36tNPP1VoaKitrWnTppo/fz43AAMAgHJVopWZ3Nxcubm55Wt3c3NTbm7uTRcFAABQVCUKM/fcc4+efvpp/fTTT7a2H3/8Uc8884zuvffeIh9n+/btioqKUmBgoCwWi9auXWu3PyYmRhaLxe7Vrl27kpQMAAAqqBKFmXnz5ikrK0shISG67bbbdPvtt6t+/frKysrSP//5zyIfJzs7Wy1bttS8efMK7dOzZ0+lpaXZXhs2bChJyQAAoIIq0T0zQUFB+u9//6uEhAR9++23MgxDTZs2Vbdu3Yp1nMjISEVGRl63j9Vqlb+/f0nKBAAAt4Bircxs3bpVTZs2VWZmpiSpe/fuGjVqlEaPHq02bdrojjvu0GeffVaqBSYmJsrX11eNGjXS8OHDlZ6eft3+OTk5yszMtHsBAICKq1grM3PnztXw4cPl5eWVb5+3t7eeeOIJzZkzR507dy6V4iIjI/XAAw8oODhYKSkpmjx5su655x4lJyfLarUW+J74+HhNnTq1VM4POKOQ59c7ugQAcCrFWpn5+uuv1bNnz0L3R0REKDk5+aaLyjNo0CD17t1bzZo1U1RUlD755BN99913Wr++8G/mEydOVEZGhu116tSpUqsHAAA4n2KtzJw5c6bAP8m2HczVVT///PNNF1WYgIAABQcH6+jRo4X2sVqtha7aAACAiqdYKzN16tTRgQMHCt2/f/9+BQQE3HRRhTl79qxOnTpVpucAAADmUqww06tXL7300ku6dOlSvn0XL17UlClT1KdPnyIf78KFC9q3b5/27dsnSUpJSdG+fft08uRJXbhwQePGjdPOnTuVmpqqxMRERUVFqVatWurfv39xygYAABVYsS4zvfjii1q9erUaNWqkkSNHqnHjxrJYLDp8+LDmz5+vq1evatKkSUU+3p49e9S1a1fb9tixYyVJ0dHRWrhwoQ4cOKAlS5bo/PnzCggIUNeuXbVy5Up5enoWp2wAAFCBFSvM+Pn5aceOHXrqqac0ceJEGYYhSbJYLOrRo4cWLFggPz+/Ih+vS5cutmMUZNOmTcUpDwAA3IKK/aF5wcHB2rBhg86dO6djx47JMAw1bNhQNWrUKIv6AAAArqtEnwAsSTVq1FCbNm1KsxYAAIBiK9GzmQAAAJwFYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJiaq6MLAEpbyPPri9w3dUbvMqwEAFAeWJkBAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmRpgBAACmxrOZAAAVUnGe0ybxrDYzY2UGAACYGmEGAACYGmEGAACYGmEGAACYGmEGAACYGmEGAACYmkPDzPbt2xUVFaXAwEBZLBatXbvWbr9hGIqLi1NgYKDc3d3VpUsXHTx40DHFAgAAp+TQMJOdna2WLVtq3rx5Be6fNWuW5syZo3nz5mn37t3y9/dX9+7dlZWVVc6VAgAAZ+XQD82LjIxUZGRkgfsMw9DcuXM1adIkDRgwQJL03nvvyc/PT8uWLdMTTzxRnqUCAAAn5bT3zKSkpOj06dOKiIiwtVmtVoWHh2vHjh2Fvi8nJ0eZmZl2LwAAUHE57eMMTp8+LUny8/Oza/fz89P3339f6Pvi4+M1derUMq0NFYezfNx5cesAAPyP067M5LFYLHbbhmHka/uziRMnKiMjw/Y6depUWZcIAAAcyGlXZvz9/SX9sUITEBBga09PT8+3WvNnVqtVVqu1zOsDAADOwWlXZurXry9/f38lJCTY2i5fvqykpCR16NDBgZUBAABn4tCVmQsXLujYsWO27ZSUFO3bt081a9ZUvXr1NGbMGE2fPl0NGzZUw4YNNX36dHl4eGjw4MEOrBoAADgTh4aZPXv2qGvXrrbtsWPHSpKio6O1ePFijR8/XhcvXtSIESN07tw5tW3bVps3b5anp6ejSgYAAE7GoWGmS5cuMgyj0P0Wi0VxcXGKi4srv6IAAICpOO09MwAAAEVBmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKZGmAEAAKbm6ugCADMJeX69o0sAAFyDlRkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqPM4ApsBjBABIfC9AwViZAQAApkaYAQAApkaYAQAApkaYAQAApkaYAQAApkaYAQAApkaYAQAApkaYAQAApkaYAQAApkaYAQAApkaYAQAApsazmQAAKKbiPCMqdUbvMqwEEiszAADA5AgzAADA1AgzAADA1AgzAADA1AgzAADA1AgzAADA1Jw6zMTFxclisdi9/P39HV0WAABwIk7/OTN33HGHPv30U9u2i4uLA6sBAADOxunDjKurK6sxAACgUE59mUmSjh49qsDAQNWvX18PPvigTpw4cd3+OTk5yszMtHsBAICKy6lXZtq2baslS5aoUaNGOnPmjF5++WV16NBBBw8elI+PT4HviY+P19SpU8u5UgBASRXn0QBAQZx6ZSYyMlL33Xefmjdvrm7dumn9+j/+wb/33nuFvmfixInKyMiwvU6dOlVe5QIAAAdw6pWZa1WtWlXNmzfX0aNHC+1jtVpltVrLsSoAAOBITr0yc62cnBwdPnxYAQEBji4FAAA4CacOM+PGjVNSUpJSUlL05Zdf6v7771dmZqaio6MdXRoAAHASTn2Z6YcfftBDDz2kX375RbVr11a7du20a9cuBQcHO7o0AADgJJw6zKxYscLRJQAAACfn1JeZAAAAboQwAwAATI0wAwAATI0wAwAATM2pbwBGxcXHlwMASgsrMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNR4NhNKDc9bAmBmZvweVpyaU2f0LsNKHIuVGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGqEGQAAYGo8zuAWZMaP7AYAs+KRA2WPlRkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqhBkAAGBqPM7gJvEx1QCA0uIsj5sx2882VmYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpmSLMLFiwQPXr11eVKlV011136bPPPnN0SQAAwEk4fZhZuXKlxowZo0mTJmnv3r3q3LmzIiMjdfLkSUeXBgAAnIDTh5k5c+bo0Ucf1WOPPabQ0FDNnTtXQUFBWrhwoaNLAwAATsCpw8zly5eVnJysiIgIu/aIiAjt2LHDQVUBAABn4tSPM/jll1909epV+fn52bX7+fnp9OnTBb4nJydHOTk5tu2MjAxJUmZmZpnUmJvzW5H7llUNxVWcmgEAFUNxfgY5w8+2vOMahnHDvk4dZvJYLBa7bcMw8rXliY+P19SpU/O1BwUFlUltxeE919EVAABuVWX1M6isf7ZlZWXJ29v7un2cOszUqlVLLi4u+VZh0tPT863W5Jk4caLGjh1r287NzdWvv/4qHx+fQgNQWcrMzFRQUJBOnTolLy+vcj9/ebqVxirdWuNlrBUTY62YKspYDcNQVlaWAgMDb9jXqcNM5cqVdddddykhIUH9+/e3tSckJKhv374Fvsdqtcpqtdq1Va9evSzLLBIvLy9T/6MqjltprNKtNV7GWjEx1oqpIoz1RisyeZw6zEjS2LFjNXToULVu3Vrt27fXW2+9pZMnT+rJJ590dGkAAMAJOH2YGTRokM6ePatp06YpLS1NzZo104YNGxQcHOzo0gAAgBNw+jAjSSNGjNCIESMcXUaJWK1WTZkyJd+lr4roVhqrdGuNl7FWTIy1YrqVxprHYhTlb54AAACclFN/aB4AAMCNEGYAAICpEWYAAICpEWYAAICpEWaKYfv27YqKilJgYKAsFovWrl17w/fMnz9foaGhcnd3V+PGjbVkyZJ8fVatWqWmTZvKarWqadOmWrNmTRlUXzxlMdbFixfLYrHke126dKmMRlE08fHxatOmjTw9PeXr66t+/frpyJEjN3xfUlKS7rrrLlWpUkUNGjTQG2+8ka+Ps81tWY3VGee2JGNNS0vT4MGD1bhxY1WqVEljxowpsJ+zzatUduOtKHO7evVqde/eXbVr15aXl5fat2+vTZs25evnbHNbVmN1xnm9GYSZYsjOzlbLli01b968IvVfuHChJk6cqLi4OB08eFBTp05VbGys/vOf/9j67Ny5U4MGDdLQoUP19ddfa+jQoRo4cKC+/PLLshpGkZTFWKU/PpEyLS3N7lWlSpWyGEKRJSUlKTY2Vrt27VJCQoKuXLmiiIgIZWdnF/qelJQU9erVS507d9bevXv1wgsvaPTo0Vq1apWtjzPObVmNVXK+uS3JWHNyclS7dm1NmjRJLVu2LLCPM86rVHbjlSrG3G7fvl3du3fXhg0blJycrK5duyoqKkp79+619XHGuS2rsUrON683xUCJSDLWrFlz3T7t27c3xo0bZ9f29NNPGx07drRtDxw40OjZs6ddnx49ehgPPvhgqdV6s0prrIsWLTK8vb3LoMLSlZ6ebkgykpKSCu0zfvx4o0mTJnZtTzzxhNGuXTvbthnmtrTGaoa5LcpY/yw8PNx4+umn87WbYV4No/TGWxHnNk/Tpk2NqVOn2rbNMLelNVYzzGtxsDJThnJycvKlXHd3d3311Vf6/fffJf3xm0BERIRdnx49emjHjh3lVmdpKMpYJenChQsKDg5W3bp11adPn3y/KTiDjIwMSVLNmjUL7VPYvO3Zs8dUc1taY5Wcf26LMtaiMMO8SqU3Xqlizm1ubq6ysrLs3mOGuS2tsUrOP6/FQZgpQz169NC//vUvJScnyzAM7dmzR++++65+//13/fLLL5Kk06dP53sCuJ+fX74nhTu7ooy1SZMmWrx4sdatW6fly5erSpUq6tixo44ePerg6v/HMAyNHTtWnTp1UrNmzQrtV9i8XblyxTRzW5pjdfa5LepYi8LZ51Uq3fFW1LmdPXu2srOzNXDgQFubs89taY7V2ee1uEzxOAOzmjx5sk6fPq127drJMAz5+fkpJiZGs2bNkouLi62fxWKxe59hGPnanF1RxtquXTu1a9fO9p6OHTvqzjvv1D//+U/94x//cFTpdkaOHKn9+/fr888/v2Hfgubt2nZnntvSHKuzz21xxloUzjyvUumOtyLO7fLlyxUXF6ePPvpIvr6+dvuceW5Lc6zOPq/FxcpMGXJ3d9e7776r3377TampqTp58qRCQkLk6empWrVqSZL8/f3zpf709PR8vx04u6KM9VqVKlVSmzZtnOY3gVGjRmndunXatm2b6tate92+hc2bq6urfHx8rtvHGea2tMd6LWea2+KMtSiceV6l0h/vtcw+tytXrtSjjz6qDz74QN26dbPb58xzW9pjvZYzzWtJEGbKgZubm+rWrSsXFxetWLFCffr0UaVKf3zp27dvr4SEBLv+mzdvVocOHRxR6k273livZRiG9u3bp4CAgHKuMn8dI0eO1OrVq7V161bVr1//hu8pbN5at24tNze36/Zx5NyW1VgLOo+j57YkYy0KZ5xXqezGW9B5zDq3y5cvV0xMjJYtW6bevXvn2++Mc1tWYy3oPI6e15tSXncaVwRZWVnG3r17jb179xqSjDlz5hh79+41vv/+e8MwDOP55583hg4daut/5MgR4/333ze+++4748svvzQGDRpk1KxZ00hJSbH1+eKLLwwXFxdjxowZxuHDh40ZM2YYrq6uxq5du8p7eHbKYqxxcXHGxo0bjePHjxt79+41hg0bZri6uhpffvlleQ/PzlNPPWV4e3sbiYmJRlpamu3122+/2fpcO94TJ04YHh4exjPPPGMcOnTIeOeddww3Nzfjww8/tPVxxrktq7E649yWZKyGYdj+3d91113G4MGDjb179xoHDx607XfGeTWMshtvRZnbZcuWGa6ursb8+fPt3nP+/HlbH2ec27IaqzPO680gzBTDtm3bDEn5XtHR0YZhGEZ0dLQRHh5u63/o0CGjVatWhru7u+Hl5WX07dvX+Pbbb/Md99///rfRuHFjw83NzWjSpImxatWqchpR4cpirGPGjDHq1atnVK5c2ahdu7YRERFh7NixoxxHVbCCxinJWLRoka3PteM1DMNITEw0wsLCjMqVKxshISHGwoUL8x3b2ea2rMbqjHNb0rEW9J7g4GC7Ps42r4ZRduOtKHMbHh5+3e9peZxtbstqrM44rzfDYhj//518AAAAJsQ9MwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAQJK0fft2RUVFKTAwUBaLRWvXri3zc/74448aMmSIfHx85OHhoVatWik5OblYxyDMALeAa78pffvtt2rXrp2qVKmiVq1alfi4qampslgs2rdv303XWBJHjhyRv7+/srKyHHL+0hISEqK5c+cWqW+bNm20evXqsi0It6zs7Gy1bNlS8+bNK5fznTt3Th07dpSbm5s++eQTHTp0SLNnz1b16tWLdRzXsikPQHmKiYnR+fPnC/0tKi0tTTVq1LBtT5kyRVWrVtWRI0dUrVq1Ep8jKChIaWlphT4ZvaxNmjRJsbGx8vT0dMj5HWHy5MkaN26c+vXrV+hDXIGSioyMVGRkZKH7L1++rBdffFFLly7V+fPn1axZM82cOVNdunQp0flmzpypoKAgLVq0yNYWEhJS7OPwfwJwC/D395fVarVtHz9+XJ06dVJwcLB8fHxKfFwXFxf5+/vL1bX8fy/64YcftG7dOg0bNqzcz+1IvXv3VkZGhjZt2uToUnALGjZsmL744gutWLFC+/fv1wMPPKCePXvq6NGjJTreunXr1Lp1az3wwAPy9fVVWFiY3n777WIfhzAD3AL+fJnJYrEoOTlZ06ZNk8ViUVxcnKQ/rlsPGjRINWrUkI+Pj/r27avU1FRJUlxcnN577z199NFHslgsslgsSkxMzHeZKTExURaLRZs2bVJYWJjc3d11zz33KD09XZ988olCQ0Pl5eWlhx56SL/99putPsMwNGvWLDVo0EDu7u5q2bKlPvzww+uO6YMPPlDLli1Vt25dW9v333+vqKgo1ahRQ1WrVtUdd9yhDRs22PYfOnRIvXr1UrVq1eTn56ehQ4fql19+se3Pzc3VzJkzdfvtt8tqtapevXp65ZVXbPsPHDige+65R+7u7vLx8dHjjz+uCxcu2PbHxMSoX79++r//+z8FBATIx8dHsbGx+v3332190tPTFRUVJXd3d9WvX19Lly7NN7a4uDjVq1dPVqtVgYGBGj16tG2fi4uLevXqpeXLl1/36wOUtuPHj2v58uX697//rc6dO+u2227TuHHj1KlTJ7uVleI4ceKEFi5cqIYNG2rTpk168sknNXr0aC1ZsqRYx+EyE3CLSUtLU7du3dSzZ0+NGzdO1apV02+//aauXbuqc+fO2r59u1xdXfXyyy+rZ8+e2r9/v8aNG6fDhw8rMzPT9k2rZs2a+umnnwo8R1xcnObNmycPDw8NHDhQAwcOlNVq1bJly3ThwgX1799f//znPzVhwgRJ0osvvqjVq1fbvqlt375dQ4YMUe3atRUeHl7gObZv367WrVvbtcXGxury5cvavn27qlatqkOHDtkuo6WlpSk8PFzDhw/XnDlzdPHiRU2YMEEDBw7U1q1bJUkTJ07U22+/rddee02dOnVSWlqavv32W0nSb7/9pp49e6pdu3bavXu30tPT9dhjj2nkyJFavHixrYZt27YpICBA27Zt07FjxzRo0CC1atVKw4cPl/RH4Dl16pS2bt2qypUra/To0UpPT7e9/8MPP9Rrr72mFStW6I477tDp06f19ddf243z7rvv1qxZs4o030Bp+e9//yvDMNSoUSO79pycHNsKb2pqqurXr3/d48TGxtruycnNzVXr1q01ffp0SVJYWJgOHjyohQsX6pFHHilybYQZ4BaTd1moWrVq8vf3lyS9++67qlSpkv71r3/JYrFIkhYtWqTq1asrMTFRERERcnd3V05Oju091/Pyyy+rY8eOkqRHH31UEydO1PHjx9WgQQNJ0v33369t27ZpwoQJys7O1pw5c7R161a1b99ektSgQQN9/vnnevPNNwsNM6mpqbrrrrvs2k6ePKn77rtPzZs3tx0nz8KFC3XnnXfavmnmjTsoKEjfffedAgIC9Prrr2vevHmKjo6WJN12223q1KmTJGnp0qW6ePGilixZoqpVq0qS5s2bp6ioKM2cOVN+fn6SpBo1amjevHlycXFRkyZN1Lt3b23ZskXDhw/Xd999p08++US7du1S27ZtJUnvvPOOQkND7cbg7++vbt26yc3NTfXq1dPdd99tN846dero5MmTys3N5b4ZlJvc3Fy5uLgoOTlZLi4udvvyfmmoU6eODh8+fN3j/Pn+vYCAADVt2tRuf2hoqFatWlWs2ggzAJScnKxjx47lu5H20qVLOn78eLGP16JFC9t/+/n5ycPDwy5Y+Pn56auvvpL0x6WfS5cuqXv37nbHuHz5ssLCwgo9x8WLF1WlShW7ttGjR+upp57S5s2b1a1bN9133322WpKTk7Vt27YCb3g+fvy4zp8/r5ycHN17770Fnu/w4cNq2bKlLchIUseOHZWbm6sjR47Ywswdd9xh940+ICBABw4csB3D1dXVbkWpSZMmdn+58cADD2ju3Llq0KCBevbsqV69eikqKsruviR3d3fl5uYqJydH7u7uhX6NgNIUFhamq1evKj09XZ07dy6wj5ubm5o0aVLkY3bs2FFHjhyxa/vuu+8UHBxcrNoIMwCUm5uru+66q8D7N2rXrl3s47m5udn+22Kx2G3nteXm5trOLUnr169XnTp17Pr9+abla9WqVUvnzp2za3vsscfUo0cPrV+/Xps3b1Z8fLxmz56tUaNGKTc317aKcq2AgACdOHHiumMyDMO2anWtP7dfb6yGYeTrf62goCAdOXJECQkJ+vTTTzVixAi9+uqrSkpKsh37119/lYeHB0EGpe7ChQs6duyYbTslJUX79u1TzZo11ahRIz388MN65JFHNHv2bIWFhemXX37R1q1b1bx5c/Xq1avY53vmmWfUoUMHTZ8+XQMHDtRXX32lt956S2+99VaxjsP6JADdeeedOnr0qHx9fXX77bfbvby9vSVJlStX1tWrV0v93E2bNpXVatXJkyfznTsoKKjQ94WFhenQoUP52oOCgvTkk09q9erVevbZZ21/GXHnnXfq4MGDCgkJyXeeqlWrqmHDhnJ3d9eWLVsKrXPfvn3Kzs62tX3xxReqVKlSvnsIChMaGqorV65oz549trYjR47o/Pnzdv3c3d3117/+Vf/4xz+UmJionTt32lZ3JOmbb77RnXfeWaRzAsWxZ88ehYWF2VZFx44dq7CwML300kuS/rj8/Mgjj+jZZ59V48aN9de//lVffvnldf9fvZ42bdpozZo1Wr58uZo1a6a///3vmjt3rh5++OFiHYeVGaCCyMjIyPfhdTVr1lS9evVu+N6HH35Yr776qvr27atp06apbt26OnnypFavXq3nnntOdevWVUhIiDZt2qQjR47Ix8fHFnJulqenp8aNG6dnnnlGubm56tSpkzIzM7Vjxw5Vq1bNdv/KtXr06KHHHntMV69etV3WGTNmjCIjI9WoUSOdO3dOW7dutd2PEhsbq7ffflsPPfSQnnvuOdWqVUvHjh3TihUr9Pbbb6tKlSqaMGGCxo8fr8qVK6tjx476+eefdfDgQT366KN6+OGHNWXKFEVHRysuLk4///yzRo0apaFDh9ouMd1I48aN1bNnTw0fPlxvvfWWXF1dNWbMGLsVlsWLF+vq1atq27atPDw89P7778vd3d1u2f2zzz5TRERESb/kQKG6dOliW0EsiJubm6ZOnaqpU6eW2jn79OmjPn363NQxWJkBKojExETbb1R5r7zfpm7Ew8ND27dvV7169TRgwACFhobqb3/7my5evCgvLy9J0vDhw9W4cWO1bt1atWvX1hdffFFqtf/973/XSy+9pPj4eIWGhqpHjx76z3/+c92/iujVq5fc3Nz06aef2tquXr2q2NhYhYaGqmfPnmrcuLEWLFggSQoMDNQXX3yhq1evqkePHmrWrJmefvppeXt7226inTx5sp599lm99NJLCg0N1aBBg2x/aeTh4aFNmzbp119/VZs2bXT//ffr3nvvLfYnpS5atEhBQUEKDw/XgAED9Pjjj8vX19e2v3r16nr77bfVsWNHtWjRQlu2bNF//vMf21+L/Pjjj9qxY8ct9/k6wPVYjOtFMABwYgsWLNBHH310S32A3HPPPaeMjIxi31MAVGRcZgJgWo8//rjOnTunrKysW+aRBr6+vho3bpyjywCcCiszAADA1LhnBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmNr/Bz9R/qxYXfV8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Plotting the distribution of mean lifetimes obtained from the 500 experiments\n",
    "plt.hist(mean_lifetimes, bins=40)\n",
    "#plt.savefig('outputHistogram.pdf')\n",
    "plt.title(\"Mean muon lifetimes\")\n",
    "plt.xlabel(\"Lifetime (seconds)\")\n",
    "plt.ylabel(\"Count\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a84d4e42",
   "metadata": {},
   "source": [
    "<b>4. What do we learn from the distribution?<b/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfd5de9a",
   "metadata": {},
   "source": [
    "<b>a) How well (with what numerical precision) can you expect to estimate the true\n",
    "lifetime from any single experiment?<b/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "8015bfc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Standard deviation of the gaussain distribution\n",
    "std_of_dist = mean_lifetimes.std()\n",
    "std_of_dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "id": "23554756",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The numerical precision of estimating the true lifetime for any single experiment is \n",
      "equal to the standard deviation of multiple experiments:  6.394320414005588e-08\n"
     ]
    }
   ],
   "source": [
    "print(\"\"\"The numerical precision of estimating the true lifetime for any single experiment is \n",
    "equal to the standard deviation of multiple experiments: \"\"\", std_of_dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3a7b797",
   "metadata": {},
   "source": [
    "From any single experiment we expect to estimate the true lifetime of a muon with a numerical precision that is equal to the standard deviation of the gaussian distribution. This is because the standard deviation calculated by performing numerous experiments gives the precision of estimating the true lifetime for any single experiment. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a866c75",
   "metadata": {},
   "source": [
    "<b>b) Is the method biased? What is the bias and precision on the bias?<b/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "id": "05f700d1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0923154147048713e-06"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Mean of the gaussain distribution\n",
    "mean_of_dist = mean_lifetimes.mean()\n",
    "mean_of_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c354b07c",
   "metadata": {},
   "source": [
    "The method is biased as it uses a truncated version (0 to 10ms) of the complete exponential distribution of the decay times of muons. The exponential distribution can otherwise produce decay times between 0ms and infinity. The estimate of the mean decay time of muons is missing longer lived muons and hence the average lifetime is pulled down and becomes biased. \n",
    "\n",
    "Further indication of the bias is given by the estimated muon lifetime lying more than 1 std away from the true lifetime.\n",
    "\n",
    "The bias is given by = |true value - estimated value|, ie the absolute difference between the true and experimentally obtained values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "89767d48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The bias of the experiment is:  1.0768458529512886e-07 seconds\n"
     ]
    }
   ],
   "source": [
    "#The method is biased and the bias is given by the difference between the actual value and the obtained value\n",
    "bias = 2.2e-6-mean_of_dist\n",
    "print(\"The bias of the experiment is: \", bias, \"seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc875862",
   "metadata": {},
   "source": [
    "The precision on the bias is given by the standard error of the mean of the distribution. This is because the precision on the bias is a measure of how precise the mean of the distribution is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "id": "ee9260d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Error on the knowledge of the mean = Standard error of the mean\n",
    "SEM = std_of_dist/np.sqrt(500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "id": "8bafc9d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Precision on the bias is:  2.8596270231262186e-09\n"
     ]
    }
   ],
   "source": [
    "#The precision on the bias is given by the standard error of the mean of the distribution\n",
    "print(\"Precision on the bias is: \", SEM)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebbbe839",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}