{ "cells": [ { "cell_type": "markdown", "id": "94e3caf0", "metadata": {}, "source": [ "# Att anpassa en funktion till ett diagram" ] }, { "cell_type": "markdown", "id": "0c1da02a", "metadata": {}, "source": [ "I denna övning ska vi pröva på att anpassa en funktion till ett histogram. Funktionsanpassning är ett mycket användbart verktyg när man vill simplifiera statistik för beräkningar, eller jämföra resultat med en teoretisk modell. Inom partikelfysiken är detta användbart när vi vill avgöra position, höjd och bredd för en pik i ett histogram.\n", "\n", "Innan du tacklar detta exempel rekommederas att du har bekantat dig med dokumentet **Intro 3: Funktionsanpassning i Python**.\n", "\n", "Vi utgår från CMS' partikelfysikdata gällande sådana kollisioner där två myoner uppstår. Datafilen *Dimuon_DoubleMu.csv* är hämtad från CERN:s opendata-portal och kan hittas i samma mapp som denna notebook." ] }, { "cell_type": "markdown", "id": "8fef6da1", "metadata": {}, "source": [ "## Hämta och rita data" ] }, { "cell_type": "markdown", "id": "7bf75996", "metadata": {}, "source": [ "Vi börjar med att ta in de relevanta funktionspaketen och läsa in datafilen." ] }, { "cell_type": "code", "execution_count": 1, "id": "8a0a14e5", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
RunEventtype1E1px1py1pz1pt1eta1phi1...type2E2px2py2pz2pt2eta2phi2Q2M
016561774601703G9.6987-9.51040.36621.86339.51750.19453.1031...G9.76337.3277-1.15246.34737.41780.7756-0.1560117.4922
116561775100943G6.2039-4.26660.4565-4.47934.2910-0.91213.0350...G9.66907.2740-2.8211-5.71047.8019-0.6786-0.3700111.5534
216561775587682G19.2892-4.2121-0.651618.81214.26222.1905-2.9881...G9.82444.3439-0.47358.79854.36971.4497-0.108619.1636
316561775660978G7.0427-6.3268-0.26853.08026.33250.4690-3.0992...G5.58574.47480.8489-3.23194.5546-0.66050.1875112.4774
416561775947690G7.27510.1030-5.5331-4.72125.5340-0.7736-1.5522...G7.3181-0.39886.94082.28256.95230.32271.6282114.3159
\n", "

5 rows × 21 columns

\n", "
" ], "text/plain": [ " Run Event type1 E1 px1 py1 pz1 pt1 eta1 \\\n", "0 165617 74601703 G 9.6987 -9.5104 0.3662 1.8633 9.5175 0.1945 \n", "1 165617 75100943 G 6.2039 -4.2666 0.4565 -4.4793 4.2910 -0.9121 \n", "2 165617 75587682 G 19.2892 -4.2121 -0.6516 18.8121 4.2622 2.1905 \n", "3 165617 75660978 G 7.0427 -6.3268 -0.2685 3.0802 6.3325 0.4690 \n", "4 165617 75947690 G 7.2751 0.1030 -5.5331 -4.7212 5.5340 -0.7736 \n", "\n", " phi1 ... type2 E2 px2 py2 pz2 pt2 eta2 phi2 \\\n", "0 3.1031 ... G 9.7633 7.3277 -1.1524 6.3473 7.4178 0.7756 -0.1560 \n", "1 3.0350 ... G 9.6690 7.2740 -2.8211 -5.7104 7.8019 -0.6786 -0.3700 \n", "2 -2.9881 ... G 9.8244 4.3439 -0.4735 8.7985 4.3697 1.4497 -0.1086 \n", "3 -3.0992 ... G 5.5857 4.4748 0.8489 -3.2319 4.5546 -0.6605 0.1875 \n", "4 -1.5522 ... G 7.3181 -0.3988 6.9408 2.2825 6.9523 0.3227 1.6282 \n", "\n", " Q2 M \n", "0 1 17.4922 \n", "1 1 11.5534 \n", "2 1 9.1636 \n", "3 1 12.4774 \n", "4 1 14.3159 \n", "\n", "[5 rows x 21 columns]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "data = pd.read_csv('Dimuon_DoubleMu.csv')\n", "data.head()" ] }, { "cell_type": "markdown", "id": "88caeb4a", "metadata": {}, "source": [ "Sedan ritar vi in statistiken från kolumnen **M** (Myonernas sammanlagda massa, mätt i GeV) i ett histogram." ] }, { "cell_type": "code", "execution_count": 2, "id": "51642fb7", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAARJ0lEQVR4nO3df6zddX3H8edrLSKKTJAL6dpmt27NNiCbQsO6uZhl6Ki4rCwZSZc4moWkCcFNly1Licl0fzTBZXOTZJB04CiO2DXoQjOCk1SNWULAi6BQatcqDq509Dqn4pKh4Ht/nE/n4fbcS+85995zzr3PR3Ly/Z739/s59/PJh57X+f44h1QVkiT9xLA7IEkaDQaCJAkwECRJjYEgSQIMBElSs3bYHejXhRdeWJOTk8PuhiSNlUcfffRbVTXRa9vYBsLk5CRTU1PD7oYkjZUk/zHXNk8ZSZIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkoBVHgiTu+8fdhckaWSs6kCQJP2YgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSc2rBkKSjyU5meTJrtoFSR5Mcqwtz+/adnOS40mOJrm6q35FkifatluTpNXPTvJPrf5wkslFHqMk6QycyRHCXcC2WbXdwKGq2gwcas9JcgmwA7i0tbktyZrW5nZgF7C5PU695g3Af1fVzwJ/A3y438FIkvr3qoFQVV8Avj2rvB3Y19b3Add21fdX1YtV9TRwHLgyyTrgvKp6qKoKuHtWm1OvdS9w1amjB0nS8un3GsLFVXUCoC0vavX1wLNd+0232vq2Prv+ijZV9RLwXeBNvf5okl1JppJMzczM9Nl1SVIvi31Rudcn+5qnPl+b04tVe6tqS1VtmZiY6LOLkqRe+g2E59tpINryZKtPAxu79tsAPNfqG3rUX9EmyVrgJzn9FJUkaYn1GwgHgZ1tfSdwX1d9R7tzaBOdi8ePtNNKLyTZ2q4PXD+rzanX+l3gs+06gyRpGa19tR2SfAL4deDCJNPAB4FbgANJbgCeAa4DqKrDSQ4ATwEvATdV1cvtpW6kc8fSOcAD7QFwJ/DxJMfpHBnsWJSRSZIW5FUDoap+b45NV82x/x5gT4/6FHBZj/r/0gJFkjQ8flNZkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqBgqEJH+c5HCSJ5N8Islrk1yQ5MEkx9ry/K79b05yPMnRJFd31a9I8kTbdmuSDNIvSdLC9R0ISdYDfwRsqarLgDXADmA3cKiqNgOH2nOSXNK2XwpsA25Lsqa93O3ALmBze2zrt1+SpP4MespoLXBOkrXA64DngO3AvrZ9H3BtW98O7K+qF6vqaeA4cGWSdcB5VfVQVRVwd1cbSdIy6TsQquqbwF8BzwAngO9W1WeAi6vqRNvnBHBRa7IeeLbrJaZbbX1bn10/TZJdSaaSTM3MzPTbdUlSD4OcMjqfzqf+TcBPAa9P8p75mvSo1Tz104tVe6tqS1VtmZiYWGiXJUnzGOSU0TuAp6tqpqp+CHwK+FXg+XYaiLY82fafBjZ2td9A5xTTdFufXZckLaNBAuEZYGuS17W7gq4CjgAHgZ1tn53AfW39ILAjydlJNtG5ePxIO630QpKt7XWu72ojSVoma/ttWFUPJ7kX+BLwEvAYsBc4FziQ5AY6oXFd2/9wkgPAU23/m6rq5fZyNwJ3AecAD7SHJGkZ9R0IAFX1QeCDs8ov0jla6LX/HmBPj/oUcNkgfZEkDcZvKkuSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQbCHCZ33z/sLkjSsjIQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkZlUGgreUStLpVmUgSJJOZyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSgAEDIckbk9yb5KtJjiT5lSQXJHkwybG2PL9r/5uTHE9yNMnVXfUrkjzRtt2aJIP0S5K0cIMeIXwU+HRV/TzwS8ARYDdwqKo2A4fac5JcAuwALgW2AbclWdNe53ZgF7C5PbYN2C9J0gL1HQhJzgPeDtwJUFU/qKrvANuBfW23fcC1bX07sL+qXqyqp4HjwJVJ1gHnVdVDVVXA3V1tJEnLZJAjhDcDM8A/JHksyR1JXg9cXFUnANryorb/euDZrvbTrba+rc+unybJriRTSaZmZmYG6LokabZBAmEtcDlwe1W9Ffgf2umhOfS6LlDz1E8vVu2tqi1VtWViYmKh/ZUkzWOQQJgGpqvq4fb8XjoB8Xw7DURbnuzaf2NX+w3Ac62+oUddkrSM+g6EqvpP4NkkP9dKVwFPAQeBna22E7ivrR8EdiQ5O8kmOhePH2mnlV5IsrXdXXR9VxtJ0jJZO2D7PwTuSfIa4OvAH9AJmQNJbgCeAa4DqKrDSQ7QCY2XgJuq6uX2OjcCdwHnAA+0hyRpGQ0UCFX1OLClx6ar5th/D7CnR30KuGyQvkiSBuM3lSVJgIEgSWpWbSBM7r5/2F2QpJGyagNBkvRKBoIkCTAQJEmNgSBJAgwESVJjIMzDO5EkrSYGgiQJMBAkSY2BIEkCDARJUmMgSJIAA+FVeaeRpNXCQJAkAQbCGfEoQdJqYCBIkgADQZLUGAiSJMBAkCQ1BsIZ8sKypJXOQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIELEIgJFmT5LEk/9KeX5DkwSTH2vL8rn1vTnI8ydEkV3fVr0jyRNt2a5IM2i9J0sIsxhHC+4AjXc93A4eqajNwqD0nySXADuBSYBtwW5I1rc3twC5gc3tsW4R+LbrJ3ff7BTVJK9ZAgZBkA/Bu4I6u8nZgX1vfB1zbVd9fVS9W1dPAceDKJOuA86rqoaoq4O6uNpKkZTLoEcLfAn8G/KirdnFVnQBoy4tafT3wbNd+0622vq3Prp8mya4kU0mmZmZmBuy6JKlb34GQ5LeAk1X16Jk26VGreeqnF6v2VtWWqtoyMTFxhn9WknQm1g7Q9m3Abye5BngtcF6SfwSeT7Kuqk6000En2/7TwMau9huA51p9Q4+6JGkZ9X2EUFU3V9WGqpqkc7H4s1X1HuAgsLPtthO4r60fBHYkOTvJJjoXjx9pp5VeSLK13V10fVcbSdIyGeQIYS63AAeS3AA8A1wHUFWHkxwAngJeAm6qqpdbmxuBu4BzgAfaQ5K0jBYlEKrq88Dn2/p/AVfNsd8eYE+P+hRw2WL0RZLUH7+pLEkCDARJUmMg9MFvK0taiQwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCB0De/iyBppTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCB0NOZfgvZbytLWkkMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBsLAvNNI0kphIEiSAANBktQYCJIkYIBASLIxyeeSHElyOMn7Wv2CJA8mOdaW53e1uTnJ8SRHk1zdVb8iyRNt261JMtiwJEkLNcgRwkvAn1TVLwBbgZuSXALsBg5V1WbgUHtO27YDuBTYBtyWZE17rduBXcDm9tg2QL8kSX3oOxCq6kRVfamtvwAcAdYD24F9bbd9wLVtfTuwv6perKqngePAlUnWAedV1UNVVcDdXW3GgncaSVoJFuUaQpJJ4K3Aw8DFVXUCOqEBXNR2Ww8829VsutXWt/XZ9V5/Z1eSqSRTMzMzi9F1SVIzcCAkORf4JPD+qvrefLv2qNU89dOLVXuraktVbZmYmFh4ZyVJcxooEJKcRScM7qmqT7Xy8+00EG15stWngY1dzTcAz7X6hh51SdIyGuQuowB3Akeq6iNdmw4CO9v6TuC+rvqOJGcn2UTn4vEj7bTSC0m2tte8vqvN2PA6gqRxN8gRwtuA3wd+I8nj7XENcAvwziTHgHe251TVYeAA8BTwaeCmqnq5vdaNwB10LjR/DXhggH4tiG/kktSxtt+GVfVv9D7/D3DVHG32AHt61KeAy/rtiyRpcH5TWZIEGAiLytNPksaZgSBJAgwESVJjICwyTxtJGlcGgiQJMBAkSY2BIEkCDARJUmMgzLIYF4W9sCxpHBkIkiTAQFgyHiVIGjd9/7jdSuKbtyR5hLCkDBpJ48RAkCQBBsKS8yhB0rgwELQsDEZp9BkIy2C1vxmu9vFL48JAkCQBBsKy8VOypFFnICyj1RgKq3HM0rgyEJbZan6DXM1jl8aBgSBJAgyEoZjcfb+fliWNHANBkgQYCEO10o8SVvr4pJXGQBgyTx9JGhUGwohYLcGwGsYojSsDYcSslDfM+caxUsYorTT+D3JGUPcb5jduefcQe9If3/Cl8TQygZBkG/BRYA1wR1XdMuQujYTZb66jHhCGgTS+RiIQkqwB/g54JzANfDHJwap6arg9Gz3zveEOIyz6DYBT7UY94KTVZCQCAbgSOF5VXwdIsh/YDhgICzCOn8579bk7JCZ3329oSMtkVAJhPfBs1/Np4Jdn75RkF7CrPf1+kqN9/r0LgW/12XbUrLix5MOvLM5+PiZW3LwMuxOLZKWMZZBx/PRcG0YlENKjVqcVqvYCewf+Y8lUVW0Z9HVGgWMZTY5lNK2UsSzVOEblttNpYGPX8w3Ac0PqiyStSqMSCF8ENifZlOQ1wA7g4JD7JEmrykicMqqql5K8F/hXOredfqyqDi/hnxz4tNMIcSyjybGMppUyliUZR6pOO1UvSVqFRuWUkSRpyAwESRKwygIhybYkR5McT7J72P1ZqCTfSPJEkseTTLXaBUkeTHKsLc8fdj97SfKxJCeTPNlVm7PvSW5u83Q0ydXD6XVvc4zlQ0m+2ebm8STXdG0b5bFsTPK5JEeSHE7yvlYfu7mZZyxjNzdJXpvkkSRfbmP5i1Zf2nmpqlXxoHOx+mvAm4HXAF8GLhl2vxY4hm8AF86q/SWwu63vBj487H7O0fe3A5cDT75a34FL2vycDWxq87Zm2GN4lbF8CPjTHvuO+ljWAZe39TcA/976PHZzM89Yxm5u6Hw369y2fhbwMLB1qedlNR0h/P/PY1TVD4BTP48x7rYD+9r6PuDa4XVlblX1BeDbs8pz9X07sL+qXqyqp4HjdOZvJMwxlrmM+lhOVNWX2voLwBE6vxwwdnMzz1jmMspjqar6fnt6VnsUSzwvqykQev08xnz/sYyiAj6T5NH2Mx4AF1fVCej8gwAuGlrvFm6uvo/rXL03yVfaKaVTh/JjM5Ykk8Bb6XwaHeu5mTUWGMO5SbImyePASeDBqlryeVlNgXBGP48x4t5WVZcD7wJuSvL2YXdoiYzjXN0O/AzwFuAE8NetPhZjSXIu8Eng/VX1vfl27VEbqfH0GMtYzk1VvVxVb6Hzyw1XJrlsnt0XZSyrKRDG/ucxquq5tjwJ/DOdQ8Lnk6wDaMuTw+vhgs3V97Gbq6p6vv0D/hHw9/z4cH3kx5LkLDpvoPdU1adaeSznptdYxnluAKrqO8DngW0s8byspkAY65/HSPL6JG84tQ78JvAknTHsbLvtBO4bTg/7MlffDwI7kpydZBOwGXhkCP07Y6f+kTa/Q2duYMTHkiTAncCRqvpI16axm5u5xjKOc5NkIskb2/o5wDuAr7LU8zLsq+nLfOX+Gjp3HnwN+MCw+7PAvr+Zzl0EXwYOn+o/8CbgEHCsLS8Ydl/n6P8n6Byu/5DOp5kb5us78IE2T0eBdw27/2cwlo8DTwBfaf84143JWH6NzqmFrwCPt8c14zg384xl7OYG+EXgsdbnJ4E/b/UlnRd/ukKSBKyuU0aSpHkYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUvN/l8SETVGQ6HsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.hist(data['M'], bins=300)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "710365d4", "metadata": {}, "source": [ "Vi märker att datan har en liten pik mellan 80 och 100 GeV, så vi vill rita ett nytt histogram över detta intervall. Obs: Denna gång sparar vi `plt.hist()`-funktionens värden. Funktionen återger tre värden, och vi vill använda de två första; staplarnas höjder (antal händelser) och gränser.\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "8a25cf09", "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcY0lEQVR4nO3dfZQddZ3n8feHoDw7CjRsCMQObNAFViL2MjOyMiAPAiLIDEIYRkHRyBEWccZZgriKzHAGRtHd1TEYBAMcDCAIgjA8Hh5klqcOJCQhYAJEDGSS5mF4EA+a8N0/6teVonNvd3Wn69bt25/XOffcur+qX9U31ZX7vb/6Vf1KEYGZmRnARnUHYGZm7cNJwczMck4KZmaWc1IwM7Ock4KZmeWcFMzMLFdZUpC0k6S7JC2RtFjSl1P51pJul7Q0vb+nUOdMScskPSnpY1XFZmZmjamq+xQkTQQmRsQjkrYC5gGfBE4EXoqI8yTNBN4TEWdI2g2YC+wN7ADcAewaEWsrCdDMzNazcVUrjoiVwMo0/ZqkJcAk4Ehgv7TYpcDdwBmp/MqIeBN4RtIysgRxf7NtbLvtttHd3V3Rv8DMrDPNmzfvhYjoajSvsqRQJKkb+CDwILB9ShhExEpJ26XFJgEPFKqtSGVNdXd309vbO/oBm5l1MEm/aTav8o5mSVsC1wKnR8Srgy3aoGy9c1uSZkjqldTb19c3WmGamRkVJwVJ7yBLCFdExM9T8arU39Df77A6la8AdipU3xF4fuA6I2J2RPRERE9XV8PWj5mZjVCVVx8JuBhYEhHfLcy6ATghTZ8A/KJQPl3SJpKmAFOBh6qKz8zM1ldln8I+wKeBhZLmp7KvAecBV0s6CXgW+BRARCyWdDXwOLAGOMVXHpmZtVaVVx/dR+N+AoADmtQ5Fzi3qpjMzGxwvqPZzMxyTgpmZpZzUjAzs5yTgpmZ5VpyR7OZDU/3zJvy6eXnfbzGSGy8cUvBzMxyTgpmZpbz6SOzMcSnlaxqbimYmVnOScHMzHJOCmZmlnNSMDOznJOCmZnlnBTMzCznpGBmZjknBTMzyzkpmJlZzknBzMxylSUFSZdIWi1pUaHsKknz02t5/7ObJXVL+n1h3oVVxWVmZs1VOfbRHOAHwGX9BRFxbP+0pAuAVwrLPxUR0yqMx8zMhlBZUoiIeyV1N5onScAxwEer2r6ZmQ1fXX0KHwFWRcTSQtkUSY9KukfSR2qKy8xsXKtr6OzjgLmFzyuByRHxoqQPAddL2j0iXh1YUdIMYAbA5MmTWxKsWVU8FLa1m5a3FCRtDPwlcFV/WUS8GREvpul5wFPAro3qR8TsiOiJiJ6urq5WhGxmNm7UcfroQOCJiFjRXyCpS9KENL0zMBV4uobYzMzGtSovSZ0L3A+8T9IKSSelWdN5+6kjgH2BxyQtAK4BTo6Il6qKzczMGqvy6qPjmpSf2KDsWuDaqmIxM7Ny/Ixmsw7gDmsbLR7mwszMcm4pmLW5YivArGpuKZiZWc5JwczMck4KZmaWc1IwM7OcO5rNWsCXjNpY4aRg1iZ8lZG1A58+MjOznJOCmZnlnBTMzCznpGBmZjknBTMzyzkpmJlZzknBzMxyTgpmZpZzUjAzs1yVz2i+RNJqSYsKZWdLek7S/PQ6rDDvTEnLJD0p6WNVxWVmZs1VOczFHOAHwGUDyr8XEd8pFkjaDZgO7A7sANwhadeIWFthfGaVqnrYCg+LYVWorKUQEfcCL5Vc/Ejgyoh4MyKeAZYBe1cVm5mZNVbHgHinSvoM0Av8XUS8DEwCHigssyKVmXUc/8K3dtbqjuZZwC7ANGAlcEEqV4Nlo9EKJM2Q1Cupt6+vr5IgzYaje+ZN+ctsrGtpUoiIVRGxNiLeAi5i3SmiFcBOhUV3BJ5vso7ZEdETET1dXV3VBmxmNs60NClImlj4eBTQf2XSDcB0SZtImgJMBR5qZWxmZlZhn4KkucB+wLaSVgDfBPaTNI3s1NBy4IsAEbFY0tXA48Aa4BRfeWRm1nqVJYWIOK5B8cWDLH8ucG5V8ZiZ2dB8R7OZmeWcFMzMLOekYGZmOScFMzPLOSmYmVmujmEuzDqW72q2sc4tBTMzyzkpmJlZzqePzDpM8RTW8vM+XmMkNha5pWBmZrlBk4KkjYqP0zQzs842aFJIQ1wvkDS5RfGYmVmNyvQpTAQWS3oI+F1/YUQcUVlUZmZWizJJ4VuVR2FmZm1hyKQQEfdIei8wNSLukLQ5MKH60MzMrNWGvPpI0heAa4AfpaJJwPUVxmRmZjUpc0nqKcA+wKsAEbEU2K7KoMzMrB5lksKbEfGH/g+SNiZ7nKaZmXWYMh3N90j6GrCZpIOALwE3DlVJ0iXA4cDqiNgjlX0b+ATwB+Ap4LMR8R+SuoElwJOp+gMRcfJw/zFmo813B9t4U6alMBPoAxYCXwRuBr5eot4c4JABZbcDe0TEB4BfA2cW5j0VEdPSywnB2k73zJvyl1mnKnP10VvARcBFkrYGdoyIIU8fRcS9qQVQLLut8PEB4OjhhWvWHpwYrFMNmRQk3Q0ckZadD/RJuici/nYDt/054KrC5ymSHiXr0P56RPxqA9dvZgU+FWZllOlT+JOIeFXS54GfRMQ3JT22IRuVdBawBrgiFa0EJkfEi5I+BFwvafeIeLVB3RnADIDJkz36htlg3KKx4SrTp7CxpInAMcAvN3SDkk4g64A+vv80VES8GREvpul5ZJ3QuzaqHxGzI6InInq6uro2NBwzMysokxTOAW4FlkXEw5J2BpaOZGOSDgHOAI6IiDcK5V2SJqTpnYGpwNMj2YaZmY1cmY7mnwE/K3x+GviroepJmgvsB2wraQXwTbKrjTYBbpcE6y493Rc4R9IaYC1wckS8NOx/jZmZbZCmSUHS9xnkJrWIOG2wFUfEcQ2KL26y7LXAtYOtz8zMqjdYS6G3ZVGYmVlbaJoUIuLS4mdJW0TE75otb2ZmY1+ZUVL/XNLjZMNQIGlPST+sPDIzM2u5Mlcf/W/gY0D/JaMLyDqGzcysw5RJCkTEbwcUra0gFjMzq1mZO5p/K+nDQEh6J3Aa6VSSmZl1ljIthZPJHrQzCVgBTEufzcysw5S5ee0F4PgWxGJmZjWr7OY1MzMbewY7fdQLzAM2BfYiG+9oKdnpI3c0m5l1oCFvXpN0IrB/RPwxfb4QuK1ZPTMzG7vKdDTvAGxV+LxlKjMzsw5T5pLU84BHJd2VPv8FcHZlEZmZWW3KXH30E0n/CvxpKpoZEf9ebVhmZlaHUnc0AxOAPuBlYFdJHubCzKwDDdlSkHQ+cCywGHgrFQdwb4VxmZlZDcr0KXwSeF9EvFlxLGZmVrMyp4+eBt5RdSBmZla/Mi2FN4D5ku4E8tbCUHc0S7oEOBxYHRF7pLKtgauAbmA5cExEvJzmnQmcRHZj3GkRcetw/zFmZXXPvCmfXn7ex2uMpB7j/d9vzZVpKdwA/APw/8jucO5/DWUOcMiAspnAnRExFbgzfUbSbsB0YPdU54eSJpTYhpmZjaIyl6ReOtQyTerdK6l7QPGRwH5p+lLgbuCMVH5l6rd4RtIyYG/g/pFs28zKc6vBispcfTQV+CdgN7JxkACIiJ1HsL3tI2Jlqr9S0napfBLwQGG5FanMzMxaqOnpI0n3pcmfALOANcD+wGXA5aMchxqUNRyhVdIMSb2Sevv6+kY5DDOz8W2wPoXD0vtmEXEnoIj4TUScDXx0hNtbJWkiQHpfncpXADsVltsReL7RCiJidkT0RERPV1fXCMMwM7NGBksKP03vb0oSsFTSqZKOArYbpN5gbgBOSNMnAL8olE+XtImkKcBU4KERbsNsWLpn3pS/zMa7pkkhIg5Pk18hGxn1NOBDwKdZ98XelKS5ZB3F75O0QtJJZIPrHSRpKXBQ+kxELAauBh4HbgFOiQg/s8HMrMXKXH30YJp8Dfhs2RVHxHFNZh3QZPlzgXPLrt/MzEZfmauPdgW+SnbDWb58RIy0X8HMzNpUmTuafwZcCPwYP4bTOpz7FWy8K5MU1kTErMojMTOz2jVNCmmcIoAbJX0JuI63j330UsWxmZlZiw3WUphHdgNZ/41lf1+YF8BI7mg2M7M21jQpRMSUVgZiZmb1K/s4TjMzGwecFMzMLOekYGZmucGuPtprsIoR8cjoh2NmZnUa7OqjCwaZF4x8pFQzM2tTg119tH8rAzGz+vkpbFbmjmYk7cH6T167rKqgzMysHmUGxPsm2XOVdwNuBg4F7iN7ApuZmXWQMlcfHU023PW/R8RngT2BTSqNyszMalEmKfw+It4C1kh6F9kjND3EhZlZByrTp9Ar6d3ARWTjIb2OH5VpZtaRyjx57Utp8kJJtwDviojHqg3LzMzqMOTpI0l39k9HxPKIeKxYNlyS3idpfuH1qqTTJZ0t6blC+WEj3YaZmY3MYHc0bwpsDmwr6T2sG0L7XcAOI91gRDwJTEvbmAA8R/ashs8C34uI74x03WZmtmEGO330ReB0sgQwj3VJ4VXgX0Zp+wcAT0XEbyQNubCZmVWr6emjiPg/6ZkKX42InSNiSnrtGRE/GKXtTwfmFj6fKukxSZek1omZmbXQkH0KEfF9SR+W9NeSPtP/2tANS3oncATws1Q0C9iF7NTSSpqMvSRphqReSb19fX0bGoaZmRWUuaP5crIv6/nA2lQcbPgdzYcCj0TEKoD+97TNi4BfNqoUEbOB2QA9PT2xgTGYmVlBmfsUeoDdImK0v4CPo3DqSNLEiFiZPh4FLBrl7ZmZ2RDKJIVFwH8iO6UzKiRtDhxE1pnd758lTSNrhSwfMM/MzFqgTFLYFnhc0kPAm/2FEXHESDcaEW8A2wwo+/RI12dmZqOjTFI4u+ogzMysPZQZ5uKe4mdJ+wB/DdzTuIZZeyo+QMaGZ+C+8wN4OlfZh+xMI0sExwDPANdWGJOZmdVksGEudiW7uew44EXgKkB+TKeZWecarKXwBPAr4BMRsQxA0ldaEpWZmdVisKTwV2QthbvSkNlXsm78I7O25YfPjw73wYxPg419dF1EHAu8H7gb+AqwvaRZkg5uUXxmZtZCZcY++l1EXBERhwM7kg13MbPqwMzMrPXKPKM5FxEvRcSPIuKjVQVkZmb1GVZSMDOzzuakYGZmOScFMzPLOSmYmVnOScHMzHKlxj4yG6t8A1Zr+cbBsc8tBTMzy7mlYGaVcKthbHJLwczMcrW0FCQtB14D1gJrIqJH0tZkw3N3kz2j+ZiIeLmO+Gzscd+B2eios6Wwf0RMi4ie9HkmcGdETAXuxOMrmZm1XDv1KRwJ7JemLyUbmfWMuoIxs+bcX9C56mopBHCbpHmSZqSy7SNiJUB6366m2MzMxq26Wgr7RMTzkrYDbpf0RNmKKYnMAJg8eXJV8ZmZjUu1tBQi4vn0vhq4DtgbWCVpIkB6X92k7uyI6ImInq6urlaFbGY2LrS8pSBpC2CjiHgtTR8MnAPcAJwAnJfef9Hq2Kw9NbuyyOeyxz73TbSfOk4fbQ9cJ6l/+z+NiFskPQxcLekk4FngUzXEZmOIL0M1G30tTwoR8TSwZ4PyF4EDWh2PmZmt4zuazcws56RgZmY5JwUzM8s5KZiZWa6dhrkwy/nKorHDf6vO4paCmZnlnBTMzCznpGBmZjn3Kdio8rAFZmObWwpmZpZzS8Fawi0Is7HBScHMWsqXsLY3nz4yM7OcWwpmVjm3DsYOJwVrG/7iMKufk4K1nDudzdqXIqLuGEasp6cnent76w7DCvxr30bDcH8s+IfG8EiaFxE9jea5o9nMzHItTwqSdpJ0l6QlkhZL+nIqP1vSc5Lmp9dhrY7NzGy8q6NPYQ3wdxHxiKStgHmSbk/zvhcR36khJjMzo4akEBErgZVp+jVJS4BJrY7DzMzWV+vVR5K6gQ8CDwL7AKdK+gzQS9aaeLnG8MYtd9qZjV+1dTRL2hK4Fjg9Il4FZgG7ANPIWhIXNKk3Q1KvpN6+vr5WhWtmNi7UkhQkvYMsIVwRET8HiIhVEbE2It4CLgL2blQ3ImZHRE9E9HR1dbUuaDOzcaDlp48kCbgYWBIR3y2UT0z9DQBHAYtaHdtY4FM7ZlalOvoU9gE+DSyUND+VfQ04TtI0IIDlwBdriM3MbFyr4+qj+wA1mHVzq2PpVFW1JtxKsTr4uGstj300Tg0cjqLZf7Zmw1b4P6rVzUOqVMNJwTaY/3PaaPMxVR8nhTHMv9bNbLQ5KXQI/7Ky8WC4x7l/OA2fR0k1M7OcWwpJs18gG/LrYrR+pbgVYGat4qTQIlUknQ3Z7kiXM7PO5qQwDFV/sfv8p1l1/P+rHCeFIdT1C9q/3M1ar0ziGOz/ZickGyeFmvnL38zaiZOCmY07PpXUnJOCmY1rHsrl7ZwUzMxGSZmLUdo92TgpjIJ2/yObjSdV9NO1uu+vzu8UJwUzszbW6nucnBRGma8mMrOBxtL3gpOCmVlN2vHUs5OCmVkbaJfWRNuNkirpEElPSlomaWbd8ZiZjSdt1VKQNAH4F+AgYAXwsKQbIuLxKrbXLpnZzKxdtFtLYW9gWUQ8HRF/AK4Ejqw5JjOzcaPdksIk4LeFzytSmZmZtUBbnT4C1KAs3raANAOYkT6+LunJDdjetsALG1C/Ko5reBzX8Diu4WnLuHT+BsX13mYz2i0prAB2KnzeEXi+uEBEzAZmj8bGJPVGRM9orGs0Oa7hcVzD47iGZ7zF1W6njx4GpkqaIumdwHTghppjMjMbN9qqpRARaySdCtwKTAAuiYjFNYdlZjZutFVSAIiIm4GbW7S5UTkNVQHHNTyOa3gc1/CMq7gUEUMvZWZm40K79SmYmVmNOjIpSPqKpMWSFkmaK2lTSVtLul3S0vT+niZ1Kxtmo0lc35b0hKTHJF0n6d1N6i6XtFDSfEm9LYjrbEnPpe3Nl3RYk7qt3l9XFWJaLml+k7pV7q8vp5gWSzo9lbXD8dUornY4vhrF1Q7HV6O4ajm+JF0iabWkRYWypseUpDPTPnlS0searLPUMbmeiOioF9nNbs8Am6XPVwMnAv8MzExlM4HzG9SdADwF7Ay8E1gA7FZxXAcDG6ey8xvFleYtB7Zt4f46G/jqEHVbvr8GLHMB8I0W7689gEXA5mR9cncAU9vg+GoWV93HV7O46j6+GsZV1/EF7AvsBSwqlDU8poDd0r7YBJiS9tGEBusc8phs9OrIlgLZH3kzSRuT/dGfJxsu49I0/1Lgkw3qVT3MxnpxRcRtEbEmzX+A7N6MVmu0v8po+f7qnyFJwDHA3FHcXhn/BXggIt5If7d7gKOo//hqGFcbHF/N9lcZLd9f/TNbfXxFxL3ASwOKmx1TRwJXRsSbEfEMsIxsXw1U5phcT8clhYh4DvgO8CywEnglIm4Dto+IlWmZlcB2DapXNszGIHEVfQ7412arAG6TNE/ZXd2jYoi4Tk2nHS5p0vSsc399BFgVEUubrYIK9hfZr8t9JW0jaXPgMLIbLms9vgaJq6jlx9cQcdV2fA0RF9R3fBU1O6bK7pcyx+R6Oi4ppIPrSLJm1Q7AFpL+pmz1BmWjcnnWUHFJOgtYA1zRZBX7RMRewKHAKZL2rTiuWcAuwDSyL+ULGlVvUNaS/QUcx+C/4irZXxGxhOw0zO3ALWTN+DWDVlqnsv01VFx1HV+DxFXr8VXi71jL8VVSZfsFOjApAAcCz0REX0T8Efg58GFglaSJAOl9dYO6Qw6zUUFcSDoBOBw4PtIJwIEi4vn0vhq4jsbNxVGLKyJWRcTaiHgLuKjJ9uraXxsDfwlc1axyhfuLiLg4IvaKiH3JmvxLqf/4ahZX3cdXw7ja4PgabH/VenwVNDumyu6XMsfkejoxKTwL/JmkzdN5wQOAJWTDZZyQljkB+EWDulUOs9EwLkmHAGcAR0TEG40qStpC0lb902Sdh4saLTuKcU0sLHNUk+21fH+leQcCT0TEikYVK95fSNouvU8m+/KYS/3HV8O42uD4ahZX3cdXs78j1Hx8FTQ7pm4ApkvaRNIUso77h4ZRf3Ab0mPeri/gW8ATZH+oy8l66bcB7iT7NXAnsHVadgfg5kLdw4Bfk/Xon9WCuJaRnR+cn14XDoyL7OqLBem1uEVxXQ4sBB5LB9fEdthfqXwOcPKAZVu5v34FPJ7Wf0Aqa4fjq1Fc7XB8NYqrHY6v9eKq6/giS0grgT+StQROanZMpeXPSvvkSeDQQvmPgZ7BjsmhXr6j2czMcp14+sjMzEbIScHMzHJOCmZmlnNSMDOznJOCmZnlnBSsVpJeb8E2zpF04AjrTlOTETzbkaRrJO2cpreUNEvSU5IeTUMyfGGI+ncPHHVT0umSfiipS9ItVcZv9XNSsI4maUJEfCMi7hjhKqaRXSvf9iTtTjZa5tOp6MfAy2Sjf34QOATYeojVzCW7SaxoOjA3IvqAlZL2GcWwrc04KVhbkLRf+pV6jbLx/69Q5lBJVw9Y7sY0PUtSr7Lx8L9VWGa5pG9Iug/4lKQ5ko5O874h6WFl4+jPTndL9/9CPl/SQ5J+Lekj6S7ac4BjlY2bf+yAmE+UdL2kGyU9I+lUSX+bfpU/IGnrtNwX0jYXSLpW2QBsSPpUimOBpHtT2e4phvnKBoubmsqvT7/0F6v5AGzHk+5albQL2dALX49sKAkiGzLk/EL8f5/ieqyw/64BDpe0SVqmm+zGrfvS/OvTdqxTjeYdgn75NdwX8Hp63w94hWwcl42A+4H/TjZ89rPAFmm5WcDfpOn+u4YnAHcDH0iflwP/s7CNOcDRxTpp+nLgE2n6buCCNH0YcEeaPhH4QZPYTyS7Y3groCvFf3Ka9z3g9DS9TaHOPwL/I00vBCal6Xen9++TjVEE2TMENhvwb92M7A7vbRrEcw/wX9P0EcB1g+z3g8me8au0v38J7Jvm3QQcmaZnAt8u1JsELKz7uPGrupdbCtZOHoqIFZH9sp0PdEc21v0twCeUDVT2cdaN4XKMpEeAR4HdyR4+0q/ZYGb7S3pQ0kLgo6lev5+n93lAd8mY74qI1yI7tfIKcGMqX1hYxx6SfpW2eXxhm/8GzEnn+SeksvuBr0k6A3hvRPw+lZ8maQHZMxF2IhvvZqCJQF+jICWdlVof/QOnHZxejwKPAO8vrLN4Cmk6bx8tdDVZy8E61MZ1B2BW8GZhei3rjs+rgFPIRrJ8OCJeSwOBfRX4bxHxsqQ5wKaF+r8buHJJmwI/JBsb5reSzh5Qp3/7xW0PJ+a3Cp/fKqxjDvDJiFgg6USyVhERcbKkPyVLdPMlTYuIn0p6MJXdKunzaV0HAn8eEW9IuntA3P1+Xyh/HNhT0kYR8VZEnAucW+jYF/BPEfGjBuu5HviupL3IWiqPFOZtmrZjHcotBRsL7iZ7VOEXWNcCeBfZF/8rkrYnG9d+KP1fmC9I2hI4ukSd18hOD22Ircg6aN9B4Xy8pF0i4sGI+AbwArBTunLo6Yj4v2QDxX0A+BPg5ZQQ3g/8WZPtLAH+M0BELAN6gX+UNCFtb1PWjcV/K/C5tB+QNElp1NCIeJ1sn1/C+s8U2JVqRgS1NuGkYG0vItaSnfM+NL0TEQvITn0sJvvy+rcS6/kPsrH7F5L9Gn64xObvAnZr1NE8DP8LeJDsgS5PFMq/rezh74uAe8lG3TwWWKTsgfHvBy4jO322saTHgH8gO4XUyE2kVkjyebKRMpdJmkf2HOIzACJ7it1PgfvTaa1reHvymwvsSfYIzKL903asQ3mUVLMOIWkzsiS2T0qkVWzjXrJO6JerWL/Vz0nBrIMou/FsSUQ8W8G6u8gSzvWjvW5rH04KZmaWc5+CmZnlnBTMzCznpGBmZjknBTMzyzkpmJlZzknBzMxy/x9gViJaG7A6cgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "undre_grans = 80 # Pro tip: Dessa värden kan skrivas in direkt i koden nedan, men det är praktiskt att\n", "ovre_grans = 100 # samla parametrar på ett ställe ifall man vill justera något i efterhand.\n", "antal_staplar = 100 # Då behöver man inte söka genom koden för att se var man borde ändra\n", "\n", "plt.figure()\n", "handelser, granser, _ = plt.hist(data['M'], bins=antal_staplar, range=(undre_grans, ovre_grans))\n", "plt.xlabel('Invariant massa (GeV)')\n", "plt.ylabel('Antal händelser')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "c267cc73", "metadata": {}, "source": [ "## Skapa en modell för funktionen" ] }, { "cell_type": "markdown", "id": "17738f3b", "metadata": {}, "source": [ "Vi vill anpassa en funktion till fördelningen ovan. Inom partikelfysik används ofta [Breit-Wigner -fördelningen](https://en.wikipedia.org/wiki/Relativistic_Breit%E2%80%93Wigner_distribution) - en fördelning som påminner om normalfördelningen. I detta exempel använder vi dock Gauss' normalfördelning, eftersom den är bekant. \n", "\n", "Normalfördelningens allmänna funktion med standardavvikelsen $\\sigma$ och medelvärdet $\\mu$ ser ut så här:\n", "\n", "$$f(x) = \\frac{1}{\\sigma \\sqrt{2\\pi}}{e^{\\frac{-(x-\\mu)^2}{2\\sigma^2}}}$$\n", "\n", "Vi *kan* anpassa den allmänna normalfördelningen till statistiken, men vi får ett betydligt bättre resultat om vi beaktar bakgrundsstrålningen också. Vi kan anta att bakgrundsstrålningen avtar linjärt i intervallet, så som skissen nedan antyder. Vi vill alltså lägga till en linjär funktion $bx + c$ till fördelningen.\n", "\n", "\n", "\n", "Funktionen som ska anpassas blir då\n", "\n", "$$f(x) = ae^{\\frac{-(x-\\mu)^2}{2\\sigma^2}} + bx + c$$\n", "\n", "där $a$ är en koefficient som ersätter $\\frac{1}{\\sigma \\sqrt{2\\pi}}$, $\\mu$ är väntevärdet, $\\sigma^2$ variansen och $bx$ och $c$ är bakgrundsstrålningen.\n", "\n", "Vi definierar funktionen i python. Denna funktion har 5 obekanta variabler, så utan hjälpmedel skulle vi inte ha någon annan möjlighet än att gissa oss fram till en anpassning. Python klarar optimeringen fint." ] }, { "cell_type": "code", "execution_count": 4, "id": "93cbec04", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "def gauss(x, *p):\n", " a, mu, sigma, b, c = p\n", " return a*np.exp(-(x-mu)**2/(2.*sigma**2)) + b*x + c" ] }, { "cell_type": "markdown", "id": "df0e3ee8", "metadata": {}, "source": [ "## Anpassa funktionen" ] }, { "cell_type": "markdown", "id": "e6441b79", "metadata": {}, "source": [ "Nu ska vi anpassa funktionen med hjälp av **scipy**-paketets **optimize.curve_fit**-kommando. Vi börjar med att bestämma staplarnas mittvärden, så att anpassningen träffar rätt." ] }, { "cell_type": "code", "execution_count": 6, "id": "51208e99", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Väntevärde = 90.8583262531103\n", "Pikbredd = 4.493052846553603\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxa0lEQVR4nO3deZgU5dX38e/p7mGRHRkQERx2BATUCS6IAUEUd03iEp+oiQkxaox5skg0JsbEaEyMz5tNhbjFBXcFZVFAFjUIDPsu26ggwoAIItt093n/qJqhGbpnema6uno5n+vqq7vv7qr6TVHMmbuWu0RVMcYYYwACfgcwxhiTOawoGGOMqWRFwRhjTCUrCsYYYypZUTDGGFPJioIxxphKnhUFEekoIjNEZJWIrBCRn7jtrUVkqoisdZ9bxUzzKxFZJyJrRORcr7IZY4yJT7y6TkFE2gPtVXWhiDQDFgCXAtcDn6vq/SIyGmilqreLSG9gHDAQOBaYBvRQ1YgnAY0xxhwh5NWMVXULsMV9/aWIrAI6AJcAQ9yvPQXMBG53259X1QPARhFZh1Mg5iRaRps2bbSoqMijn8AYY3LTggULtqtqYbzPPCsKsUSkCDgJmAu0cwsGqrpFRNq6X+sAfBAz2Sa3LaGioiJKSkpSH9gYY3KYiHyU6DPPDzSLSFPgFeA2Vd1d3VfjtB2xb0tERolIiYiUlJWVpSqmMcYYPC4KIlKAUxCeVdVX3eat7vGGiuMO29z2TUDHmMmPAz6tOk9VHaOqxapaXFgYt/djjDGmjrw8+0iAx4BVqvrXmI8mANe5r68Dxse0XyUiDUWkM9AdmOdVPmOMMUfy8pjCIOA7wDIRWey23QHcD7woIjcAHwPfAlDVFSLyIrASCAM325lHxhiTXl6effQe8Y8TAAxLMM29wL1eZTLGGFM9u6LZGGNMJSsKxhhjKllRMMYYU8mKgjHGmEppuaLZGFM7RaMnVr4uvf8CH5OYfGM9BWOMMZWsKBhjjKlku4+MySK2W8l4zXoKxhhjKllRMMYYU8mKgjHGmEpWFIwxxlSyomCMMaaSFQVjjDGVrCgYY4ypZEXBGGNMJSsKxhhjKllRMMYYU8mzoiAij4vINhFZHtP2gogsdh+lFfduFpEiEdkX89kjXuUyxhiTmJdjHz0J/AP4T0WDql5Z8VpEHgR2xXx/vaoO8DCPMcaYGnhWFFR1togUxftMRAS4Ajjbq+UbY4ypPb+OKQwGtqrq2pi2ziKySERmichgn3IZY0xe82vo7KuBcTHvtwCdVHWHiJwCvC4ifVR1d9UJRWQUMAqgU6dOaQlrjFdsKGyTadLeUxCREHA58EJFm6oeUNUd7usFwHqgR7zpVXWMqharanFhYWE6IhtjTN7wY/fRcGC1qm6qaBCRQhEJuq+7AN2BDT5kM8aYvObZ7iMRGQcMAdqIyCbgt6r6GHAVh+86AjgLuEdEwkAEuFFVP/cqmzGZojl7+EbwXbrLZnjiX7D/C+h3JU1pzx6O8jueyUNenn10dYL26+O0vQK84lUWYzLSpgVMangHx8l2dmpTiPSGBk1g6l38t2Fjnomcw/8LX+53SpNn7B7NxqSbKswbA2/dCbTgsgO/Y5F2p/T77oHmTxcx6+HR3BSa4PQgIhdCsKDaWdoBa5MqNsyFMek26wGY/EvoNowLDvyRRdr98M+PPYkfl9/KXeXXc05wAbz+I4hG/Mlq8o71FIxJp42zYeZ90O9KuPQRdi2dnPCrT0dG0Ix9/HLZCzy3aAd3hG8AJH1ZTV6yomBMuuzZBq98H9p0hwv+CoGaO+r/ilxCU9nHTaEJLNGuvBAZmoagJp/Z7iNj0iEahVdHwf5d8K0noWHTpCd9IHwlc6O9GB0aR2uOuJ7TmJSyomBMOsx7FDbMgJEPQLs+tZxY+HX592jKPn4Ves6TeMZUsKJgjNf273IOLncZCidfW6dZrNXjGBu5gG+FZjNQVqU4oDGHWFEwxmv//Tvs+5wLVw2j6FeT6jybv4UvY5O24Q8Fj1NAOIUBjTnEDjQb46U922DOP3kzchrLtUu1X4291iCe/TTkt+XX8ViDB7kiOBO4JGUxjalgPQVjvDT7zxA+wF/C30rJ7KZHT2ZxtCujgm9CxHoLJvWsKBjjlc83QskTcPK1lGr7FM1UeDh8MccHtsGq8SmapzGHWFEwxivv/x8EgvD121M627ejp7A+2h7ee8gZMsOYFLKiYIwX9u2EpS/Cid+C5qnqJTiUAI9ELoLPlsH6d1I6b2OsKBjjhUXPQPleOPWHnsx+fGQQNDvW6S0Yk0JWFIxJtWgE5o2FTmfAMSd6soiDFMDpN0Ppu7BpgSfLMPnJioIxqbb2bfjiIzh1lLfLOeU6KGgCCx73djkmr1hRMCbV5j7q7NrpdaG3y2nYDPpeBstfown7vF2WyRtWFIxJpbIPnTGOir9X441xUuKka6H8K84PzvV+WSYveFYURORxEdkmIstj2u4Wkc0isth9nB/z2a9EZJ2IrBGRc73KZYynSh6HQAGccn16ltdxILTpwZXBmelZnsl5Xg5z8STwD+A/VdofUtW/xDaISG/gKqAPcCwwTUR6qKrdbspkj0g5LHsJep4HTQtrHLaivirm/4NgMXcWPEdX2cx67eDpMk3u86ynoKqzgc+T/PolwPOqekBVNwLrgIFeZTPGE+tnwN7t0O+qtC721chgyjXojodkTP34MSDeLSJyLVAC/ExVdwIdgA9ivrPJbTMmeyx9Hhq3gu4jqv1aqnsQO2jB9OjJXB58lz+Hr0zpvE3+SfeB5oeBrsAAYAvwoNse78azca/fF5FRIlIiIiVlZWWehDSmNopGT6Tv6JfYv2wC9LkcQg3SnuGFyBAKZTdnBxalfdkmt6S1KKjqVlWNqGoUGMuhXUSbgI4xXz0O+DTBPMaoarGqFhcWFnob2JgkjQzOo5GUQ//07jqqMDvajzJtzsXBOb4s3+SOtBYFEYkdBOYyoOLMpAnAVSLSUEQ6A92BeenMZkx9XBp4n43RdnDc13xZfoQgUyIDnZ7Cwa98yWByg5enpI4D5gA9RWSTiNwAPCAiy0RkKTAU+CmAqq4AXgRWAlOAm+3MI5MtjmEHpwdW8nrkTJB4e0LTY2L0NI6SA/DhW75lMNnPswPNqnp1nObHqvn+vcC9XuUxxiuXBt8nIMpr0TOdv3J8Mi/ai23akrYrXoW+l/uYxGQzu6LZmHq6IPgBi6Ld+Fjb+ZojSoBJkYGwdioc+NLXLCZ7WVEwpj52lnJioJTJEX+OJVT1ZuQ0CO+HNVP8jmKylBUFY+pj1RsATI5mxrWWC7SHMxjfilf9jmKylBUFY+pj5QSWR4v4xOddRxWUAPS5FNZNg/27/I5jspAVBWPqavensGkekyOZ0Uuo1OdyiByE1ZP8TmKykB/DXBiTG1a9CcCU6KHjCV4PgpeU44qhWXtY/SYMiHcSoDGJWU/BmLpaOR4Ke2XeyKQi0HMkrH8Hyu3mO6Z2rCgYUxd7yuDj/8IJF/udJL6eF0D5Xtgwy+8kJstYUTCmLla/CRqF3plXFIpGT6THY3v4UhvDGjuuYGrHioIxdbH6TWhVBO36+p0kroMUMCvaDz6cAtGo33FMFqm2KIhIIPZ2msYY4MAe2Djb2UXj41hHNZkaOQX2bIXNC/yOYrJItUXBHeJ6iYh0SlMeYzLfhhnOKZ89z/M7SbVmRAeABGFNBpwRZbJGMruP2gMrRGS6iEyoeHgdzJiMtWYyNGwBnU73O0m1dtMUigbZ9QqmVpK5TuF3nqcwJltEI87Q1N2HQ7DA7zQ163kBTLkddqyHo7v6ncZkgRp7Cqo6CygFCtzX84GFHucyJjNtXgB7t0OPkX4nSU6v851nOwvJJKnGoiAiPwBeBh51mzoAr3uYyZjMtWays5+++3C/kySnZSdo29tuvGOSlswxhZuBQcBuAFVdC7T1MpQxGevDKXD8GdC4ld9Jktd9BHw8B/bv9juJyQLJFIUDqnqw4o2IhAD1LpIxGWpnKWxbCT0y+6yjI3QfAdGwc9aUMTVI5kDzLBG5A2gsIucANwFv1DSRiDwOXAhsU9W+btufgYuAg8B64Luq+oWIFAGrgDXu5B+o6o21/WGMSbXYAe5KL/3EedEzS44nVOh4qnO21Nq3ofclfqcxGS6ZnsJooAxYBvwQmAT8OonpngSq/kk1Feirqv2AD4FfxXy2XlUHuA8rCCbjzJ74DOuj7Sn682q/o9ROMATdznZu06nWyTfVq7Gn4F7ANhYYKyKtgeNUa96yVHW22wOIbXs75u0HwDdrF9cYfzRmP6cGVvGfyAggQ4bIro3u58KK12DLEjh2gN9pTAZL5uyjmSLS3C0Ii4EnROSvKVj294DJMe87i8giEZklIoNTMH9jUuaMwAoaSti5SjgbdRtOVIW//PMf2VfQTFolc0yhharuFpHvA0+o6m9FZGl9FioidwJh4Fm3aQvQSVV3iMgpwOsi0kdVjzhdQkRGAaMAOnWy0TdMegwNLGaPNqIk2tPvKLUSWwBeb9CFs4OL+EfkMh8TmUyXzDGFkIi0B64A3qzvAkXkOpwD0NdU7IZS1QOqusN9vQDnIHSPeNOr6hhVLVbV4sLCwvrGMSYJypDgEt6P9uUgWXAVcwIzIgMYIOtphZ2aahJLpijcA7wFrFPV+SLSBVhbl4WJyHnA7cDFqro3pr1QRILu6y5Ad2BDXZZhTKp1l80cJ9uzd9eRa0Z0AAFRvh6oV0ff5LhkDjS/BLwU834D8I2aphORccAQoI2IbAJ+i3O2UUNgqjhDDlecenoWcI+IhIEIcKOqfl7rn8YYDwwNLAJgZqS/z0nqZ5l2pkyb8/XgEr+jmAyWsCiIyN+p5iI1Vb21uhmrarw7hj+W4LuvAK9UNz9j/DI0sIRV0U58xtF+R6kXJcDsaD+GBJY4N94J2D22zJGq6ymUpC2FMRmqGXspDqxhbOQCv6OkxKxIf74RfA+2LIYOJ/sdx2SghEVBVZ+KfS8iTVT1K+8jGZM5BgWWUyARZkQG+B0lJd6LnkhUhcC66VYUTFzJXKdwuoisxBmGAhHpLyL/8jyZMRlgSGAxu/UoFmp3v6OkxOc0Z6l2hnVT/Y5iMlQyOxX/DzgXqDhldAnOgWFjcps6p6LOjp5IhKDfaVJmVrQ/bJoP+3b6HcVkoKSONKnqJ1WaIh5kMSazbF3OMbLT+SWaQ2ZF+oNGYcNMv6OYDJRMUfhERM4AVEQaiMjPcXclGZPT1k0Dsv9U1KqWaFdo1KLy5zMmVjJF4UacG+10ADYBA9z3xuS2tdNYGT2eMrLohjpJiBCELkNh3XQbNdUcIZl7NG9X1WtUtZ2qtlXV/6kYksKYnLV/N3zyATNzbNdRpW7D4cstzk2DjInh2cVrxmS1DTMhGs65XUeVug1zntdOhXZ9/M1iMkp1PYUSYAHQCDgZZ7yjtTi7j+xAs8lt66ZBw+Y5cyrqEZofC237wPrpficxGabGi9dE5HpgqKqWu+8fAd5ONJ0xWU/VKQpdvk54UTKjy2epbsNg7iNwYA80bOp3GpMhkjnQfCzQLOZ9U7fNmNy0bRXs3gzdzvE7ibe6DYPIQSh9z+8kJoMkUxTuBxaJyJMi8iSwEPijp6mM8VPF1b7dhvubw2udToeCo2wXkjlMMkNnPyEik4FT3abRqvqZt7GM8dG6aVB4ArTogHMH2hwVaghFg+16BXOYZMfODQJlwE6gh4jYMBcmNx34Ej6aA91zfNdRhW7D4fMNzsMYkugpiMifgCuBFUDUbVZgtoe5jPHHxtkQLc+jouCemrpuOgzs4m8WkxGSObXiUqCnqh7wOIsx/ls7FRo0hY6n+Z0kPVp3gVZFsP4dGPgDv9OYDJDM7qMNkMV3KzcmWZWnog6BUAO/06SHCHQd5vSQwgf9TmMyQDI9hb3AYhGZDlT2Fmq6ollEHgcuBLapal+3rTXwAlAElAJXqOpO97NfATfgXBh3q6q+VdsfxphkFY2eWPm69H73rmpla2DXJzD4Zz6lSp/Dfv7rh0PJY/DJXOg82MdUJhMk01OYAPwe+C/OFc4Vj5o8CZxXpW00MF1VuwPT3feISG/gKqCPO82/RCR3BrA32aHiVNR8OZ5QofNgCITsLCQDJHdK6lM1fSfBdLNFpKhK8yXAEPf1U8BM4Ha3/Xn3uMVGEVkHDATm1GXZxtTJ2qnuqajH+Z0krYp+O5txBT1o8e5r9D7nd37HMT5L5nac3UXkZRFZKSIbKh51XF47Vd0C4D63dds7ALE38tnkthmTHgf2wMdzoHuOX7CWwMxof3oHPoLdW/yOYnyWsCiISMW1708ADwNhYCjwH+DpFOeQOG1xR2gVkVEiUiIiJWVlZSmOYfLWxtnOkA/dR/idxBeVd5ezq5vzXnU9hfPd58aqOh0QVf1IVe8Gzq7j8raKSHsA93mb274J6BjzveOAT+PNQFXHqGqxqhYXFhbWMYYxVazLs1NRq1itHflMWzm70Exeq64oPOc+HxARAdaKyC0ichmHdvvU1gTgOvf1dcD4mParRKShiHQGugPz6rgMY2qlaPSbbJo/gbf29cqfU1GPIM69mzfMgEjY7zDGRwmLgqpe6L78Kc7IqLcCpwDf4dAv9oREZBzOgeKeIrJJRG7AGVzvHBFZC5zjvkdVVwAvAiuBKcDNqmr3bDBp0UM2cZxs553oSX5H8dXMaH/Yvws2l/gdxfgombOP5rovvwS+m+yMVfXqBB8NS/D9e4F7k52/MalydmARQO7eZS1J70f7ggSdU1M75eduNJPc2Uc9RGSMiLwtIu9UPNIRzph0GBpczIro8Wyltd9RfLWbJtBxoB1XyHPJXNH8EvAI8G/sNpwmxzRnD6fIhzwcuRg4/ErfvNRtGLzzB9hTBk3tRI58lMwVzWFVfVhV56nqgoqH58mMSYOvB5YSkigzIgP8jpIZKu42Z6em5q3qrlNo7Y5V9IaI3CQi7Sva3HZjst7Q4GJ2aDMWaze/o2SGY/pBk0LbhZTHqtt9tADnArKKC8t+EfOZAjb4uslqAaIMCSxmZnQA0aTvN5XjAgGnt7BmknNqajCZPcwmlyT8F1fVzukMYky6DZB1tJY9tuuoqh4jYMlzsGk+HH+632lMmtmfRyZvDQ0uJqLCrGg/v6Nklq5nO6OmrrXR6/ORFQWTt84OLGKB9mA3Tf2OklkatYBOp8OHb/udxPjAioLJS+3ZQZ/AR0yPnOx3lMzUfQRsWwFffFLzd01Oqe7so5Ore6QzpDGpNjzonFU9NXqKz0kyVI9znee11lvIN9WdWvBgNZ8pdR8p1RjfnRNYwPpoezbosX5HyUxtekDL452i8LUb/E5j0qi6s4+GpjOIMWmzfxenBVbyeGSk30kyzmH3bj7rXFj4NJTvg4LGPqYy6ZTUMQUR6SsiV4jItRUPr4MZ45l102kgEabZ8YTqdT8Xwvug9L2av2tyRjID4v0W+Lv7GAo8AFzscS5jvLNmMju0GQu1h99JMlvRmVBwFHw4xe8kJo2S6Sl8E2e4689U9btAf6Chp6mM8UqkHNa+xTuRk+wq5poUNIIuQ2DNFNC4d8c1OSiZ/xX7VDUKhEWkOc4tNG2IC5OdPp4D+3cxzc46Sk7P82H3Jvhsqd9JTJokUxRKRKQlMBZnPKSF2K0yTbZaPQmCDZkdPdHvJNmhx3mAOOvN5IUai4Kq3qSqX6jqIzi30LzO3Y1kTHZRhTUTocsQ9tHI7zTZoWkhdDzVWW8mLyRzoLlyYHVVLVXVpbFttSUiPUVkccxjt4jcJiJ3i8jmmPbz67oMY+LasgS++BhOuMjvJNml1/nw2TJn3ZmcV90VzY3c+ya0EZFWMfdSKALqfMWPqq5R1QGqOgA4BdgLvOZ+/FDFZ6pq/VWTWqsmOPcg7nWB30myS093fa2Z7G8OkxbV9RR+iHMMoZf7XPEYD/wzRcsfBqxX1Y9SND9j4lOFleOh82A4yu4RVSttujlXOK+2XUj5IGFRUNX/595T4eeq2kVVO7uP/qr6jxQt/ypgXMz7W0RkqYg8LiKtUrQMY2DbKtixDk6wS2zqpOf58NH7sO8Lv5MYjyVzoPnvInKGiHw7lVc0i0gDnIvgXnKbHga6AgOALSQYe0lERolIiYiUlJWV1TeGyRcrxwNixxPqqtcFEA3bbTrzQDIHmp8G/gKcCXzNfRSnYNkjgYWquhVAVbeqasS9JmIsMDDeRKo6RlWLVbW4sLAwBTFMXlg1AY4/A5q29TtJdupQDE3awuo3/U5iPJbMDViLgd6qKb+k8Wpidh2JSHtV3eK+vQxYnuLlmXy1fS1sWwkjH/A7SfYKBJyzkJa+ZAPk5bhkLl5bDhyTyoWKyFE41zy8GtP8gIgsE5GlOGMs/TSVyzR5bOV459l2HdVP70uh/CtYN83vJMZDyfQU2gArRWQecKCiUVXrfMROVfcCR1dp+05d52dMtVaOh+MGQnO7d0K9FA2Gxq1hxetWYHNYMkXhbq9DGOOZ7WudcXvO/aPfSbJfMAQnXAjLX4Xy/c6AeSbn1FgUVHVW7HsRGQR8G5gVfwpjMsiylwCBPpcfdgMZUzsV625w4FiebrAH1k+3iwBzVLI32RkgIg+ISCnwB2CVp6mMSQVVpyh0Pguat/c7TU6YE+3NTm3q7EIyOam6YS56iMhvRGQV8A/gE0BUdWgKL14zxjufLoTPN8CJ3/I7Sc4IE+LtSLEz5EX4QM0TmKxTXU9hNc4wFBep6pmq+ncgkp5YxqTA0pcg2BB621XMqTQpeioc/BLWv+N3FOOB6o4pfANnGIoZIjIFeB6QtKQyph6KRk8kQJS5DZ+jsPcIaNTC70hZKdExmP9G+zjrdMXr0HNkekMZz1U39tFrqnolzoB4M3GuG2gnIg+LyIg05TOmTs4IrKBQdtmuIw+UE3JOSV09EQ7u9TuOSbFkxj76SlWfVdULgeOAxcBor4MZUx+XBN5ntzaG7uf6HSU39bvS2YX0oQ2nnWtqdedyVf1cVR9V1bO9CmRMfTXiAOcF5zMlMtDOpffK8WdC8w6w9EW/k5gUq1VRMCYbnB+YSzPZx8uRs/yOkrsCATjxm86QF19t9zuNSSErCibnXBmayYboMczTXn5HyW39rnSG017xWs3fNVnDioLJLdvXcWpgNS9GhmAny3msXR9o1xeWvuB3EpNCVhRMbln0NGEN8IrtOkqPflfApvmwY73fSUyKWFEwuSNSDoufY0b0JMpo6Xea/ND3m4C4Y0yZXJDMKKnGZIe1b8NX23g+cmgUdhsEz2MtOkDnwbBkHJz1S4ruOHSKaun9NmBeNrKegskdC/8DTY9hZnSA30nyy0nXws5SKJ3tdxKTAtZTMLnhi0+cm8oPupXItKDfafLLCRdB41aw4Eng0BXksb006zVkD+spmNwwf6zzXHyDvznyUUEj6H81rHqT1uz2O42pJ1+KgoiUuvdjXiwiJW5baxGZKiJr3edWfmQzWejgV3zx3r+ZGD6FovuX+p0mP518HUTL+UbQdiFlOz97CkNVdYCqFrvvRwPTVbU7MB0bX8kka+kLtJSveCJ8nt9J8lfbXtDpdK4OvgOo32lMPWTS7qNLgKfc108Bl/oXxWQNVfjgEZZFiyjRnn6nyRtFoydWPiqdcj1dAp9xWsBuzJjN/CoKCrwtIgtEZJTb1k5VtwC4z219ymayyfp3YPsaHg+PxK5g9lnvS9ilR/Ht4HS/k5h68Ovso0Gq+qmItAWmisjqZCd0i8gogE6dOnmVz2SLuY9Ak7ZM3H+a30lMQWNejQzmmuA02rKTbdhhwWzkS09BVT91n7cBrwEDga0i0h7Afd6WYNoxqlqsqsWFhYXpimwy0bZVzgVrX7uBgxT4ncYAT0TOI0SUa0Nv+x3F1FHaewoi0gQIqOqX7usRwD3ABOA64H73eXy6s5nMlOiq5NLiV6BBUxg4CqbMSXMqE8/H2o63o8VcE5zOP8OXsI/q72dh1zJkHj96Cu2A90RkCTAPmKiqU3CKwTkishY4x31vTFxdZTPRZa/yr71nU3SPFYRM8u/wSFrJHr4RfNfvKKYO0t5TUNUNQP847TuAYenOY7LTzaHx7KcB/w6f73cUU0WJ9mRxtAvfC07m2cgwNKNOcjQ1sX8tk3WKZAuXBN7n6chwPqe533HMEYTHwufTJfAZZwcW+R3G1JIVBZN1bgmNp5wQY8MX+h3FJDA5OpDNejSjQjZKbbaxomCyShf5lEsD7/FsZDjbaeF3HJNAmBCPhc/n1MBqThW7mC2bWFEwWeWO0LPsoyEPhy/2O4qpwbORYWzTlvxvwUvOlecmK9jQ2SYjxTsNdVBgGcODi7iv/GrrJWSQRKcMH6AB/wxfwu8KnoKNs6HL19OczNSF9RRMVggS4a7QM3wcLeSJiA18ly2ejwxli7aGGX+03kKWsKJgssKVwZn0CnzCfeFv29XLWaSit8AnHzjjVJmMZ0XBZLzmfMX/hl5ibrQXk6MD/Y5jaunFyBBofhzMuNd6C1nAioJJqbhDKtfTr0PP0JI93FN+LTYSavY5SAEMGQ2bF8Cyl/2OY2pgRcFktCGBxVwRmsUjkYtYoUV+xzF1NeDb0H4ATL0LDuzxO42phhUFkxZ16UE0Yy9/LPg3H0Y78Lfw5R6mM54LBGHkA/DlFnj3Qb/TmGrYKakmY90RepZ27OTG8p/aweVc0OlU6HcVB979G+dM78jH2s7vRCYO6ymYjDQiMJ+rQzMYE7mQpdrV7zgmVYbfTTkh7go97XcSk4D1FEzm2b6WBwseYXG0Cw+Fv+l3GpMCsbsNfxC8nDsLnuOiyH95I3qGj6lMPFYUTMYoGj2Ro9jP6w3u4mgJ8aODttsoFz0eGcnI4Dx+X/AEcw+cYLftzDC2+8ikXeKDzsoDBY/SVT7llvJb2cLRvuQz3ooQ5GflP6Ih5TxQMAawaxcyifUUjGeSOdMo9ju/DL3AhcG53Fd+NXOifbyMZny2UdtzX/hq7il4iquj7zAuMqxet+a023qmjvUUTEa4KTiem0ITeCY8jEcjdp+EfPB05Bzei/Th16Fn6C6b/I5jXGkvCiLSUURmiMgqEVkhIj9x2+8Wkc0isth92H0W88S1wbf4ZcELvBYZxF3h72JXLecHJcDPyn/EVzRmbMGDtMAuassEfvQUwsDPVPUE4DTgZhHp7X72kKoOcB+TfMhm0kq5ITiJewqe4u3IKfyi/Id2P988s5XW3HjwNo6V7fy94O8EifgdKe+l/X+gqm5R1YXu6y+BVUCHdOcw/goQ5beh/3BXwTNMjnyNH5f/mLAd4spLC7UHvw5/j7OCy7g99LzfcfKer/8LRaQIOAmYCwwCbhGRa4ESnN7ETh/j5S2vD9o1YR//V/BPzgkuZEz4Au4LX209hDz3YmQofaSUUaGJbNWWgB0s9otv/xNFpCnwCnCbqu4GHga6AgOALUDcAVJEZJSIlIhISVlZWbrimhQ5RdYwucFozg4s4q7y6/lj+BorCAaAe8LX8mbkVO4qeBbmjfU7Tt7ypacgIgU4BeFZVX0VQFW3xnw+Fngz3rSqOgYYA1BcXGwnOGeJAsLcGnqVm4Lj2axtuOLgb1igPf2OZTJIhCC3ld9MA8KMmPRzCBbAKdf7HSvvpL0oiIgAjwGrVPWvMe3tVXWL+/YyYHm6s2WD7DsfWxkZmMftoecpCmzlpfBZ/C58LXs4yu9gJgOFCXFL+a182PMZeOMnsG8nDLoNxM5ISxc/egqDgO8Ay0Rksdt2B3C1iAzAubyxFPihD9lMyiiDA8u4NfQqXwt8yJrocVx38HZmRfv7HcxkuIMUwJXPwPibYNrdsGMdXPAQhBr4HS0vpL0oqOp7xD8R3U5BTRGvehPJzLcJ+7g4+F++G5xCj8BmtmlLRpd/n5ciXydCMGVZTI4raATfeAyO7gaz/sSckoXcWn4L8++/xu9kOc/OAcxTVYegSPRLPtFQFYcViLvPhHXTYMVrLGg4hUZSzoro8fzvwRt5M3q6DWpn6kYEht7BbVO/5P6CsbzV8Jewshn0vhhIbhgVU3v5XRR2bYbwfgg1hGBDp3sabOi8t32YcQlROkoZveUjTgqs5fTASiL3lRIUZZu2ZFJkKG9ETmeB9sCuTDZ1dfgv/DNZdrAzDxX8i9Yvfgf6Xw0j/uBbtlyX30Vh0i9gTYK/NgIFbrFoUOW5ong0qOazKgUmWBCnreq0VedZcHhbsIHTFiPlu4ki5bB/F8fLZ7TmSwrlCwplF+1lB8fLVjrKNrrKFprJPgAOapBF2p1/RC7jvUhfFmgPonZ6qfHAeu3A5Qd/x7pzV8Dsv8DqifwweAFPRs7lAHasIZVENXvP6iwuLtaSkpK6z6D0fdj1CYQPQOSg+3wAwgerPMd+Huf5iLaYaTWauh8YKNcg5YQIE6ScQ68jGiBMkCgBIgSIEiCK0LdDC5y/2BXUfUTDrN/6BQWEaSBhGnGQlqGwkzeOsAbYrG34SNtRqsewUo9nZfR41mhH+w9p0q6rbOaO0HMMCy5ik7bh0fCFvBw5i300Ag7/Ayn7ztZLDxFZoKrF8T7L755C0SDvlxFxf9lGDtZQbMoPtUXDh4pLJOwWHuf1mOkrKSBCAWFC7nOQKEGJUkCYAFGCKEGiCFH6Nm2Lc0KXgASc3WKBICs/K6OcEOXREPtowPVnngANmkKjFvx0/EZ20pQybUmZtmQHze0gsckY67UDN5T/gjMiy/lF6EV+X/AkPwu9xLORYbwYGeJ3vKyX30UhRqKDVvX56yJVf6Ucnq1fraYtvSb+cn+86PCf9/pzDn3vtdfsAJ7JfP+N9uWyg305WT7k+6FJ3Bh8g5tDE2Dss3DiFdDzPL8jZiUrCmniRdGpz3Lr+j1jMs1C7cFN5T04hh1cFJzDnZFlMOV2mHI70xu0Z3a0H/OjPWFXf2jewU4iqYEVhVrw+he77f80pu4+42jGRi7kzhsfhu3rYN00Ppo4jquCM/hu6C146G98pq1YFe3E0MFfh7YnONdBtCqCJoVWLFxWFGrg11/Q9pe7MfXQphu06cb3Xu9IAWF6ycecFFjLgMB6esknMPcR51hdhYIm0KID720tYCutKdMW3DjyVGjSBhq3hsatoHFLaNSCE+59n300JN4p17nwx5wVBZ/ZL39jvFVOiGXahWWRLvzHvYdP6W/Phc83wM6NsLMUPt8Iuzdz1LaVnCYracMumBp3TE5WNXLOyPuKRnxFI/ZqI/bS0CkUzzwOBY2h4CgINXKfGzptoYZOW9XnqqerhxolPn09Db0ZKwrGmLxTdOdbla9L7z80zNrllX+kKaV3nwVflcG+L2D/Tvd5F/e9Pp9mspcm7Hceso/GHOQoOQB7P4fyvVC+z7kwtnw/hPcd3iupj6BbLG58F1p3Ts08q7CiYIzJa/F760LR3e9WvovdLfToK8cknFfpqAS7j6IRp0hUnIZe7haKiraKx2Gnqu+Pf/p6+AA0alHXH7dGVhSMMSZFqj0ZpUETaNAk408osaKQApn+j2xMPvHiOF26j/35+TvFioIxxmSwdF/jZEUhxexsImNMVdn0e8GKgjHG+CQTdz1bUTDGmAyQKb2JjBv8XkTOE5E1IrJOREb7nccYY/JJRvUURCQI/BM4B9gEzBeRCaq60ovlZUplNsaYTJFpPYWBwDpV3aCqB4HngUt8zmSMMXkj04pCB+CTmPeb3DZjjDFpkFG7j4h/p/fD7hcqIqOAUe7bPSKyph7LawNsr8f0XrFctWO5asdy1U5G5pI/1SvX8Yk+yLSisAnoGPP+OODT2C+o6hhgTCoWJiIlie5T6ifLVTuWq3YsV+3kW65M2300H+guIp1FpAFwFTDB50zGGJM3MqqnoKphEbkFeAsIAo+r6gqfYxljTN7IqKIAoKqTgElpWlxKdkN5wHLVjuWqHctVO3mVS1S15m8ZY4zJC5l2TMEYY4yPcrIoiMhPRWSFiCwXkXEi0khEWovIVBFZ6z63SjCtZ8NsJMj1ZxFZLSJLReQ1EWmZYNpSEVkmIotFpCQNue4Wkc3u8haLyPkJpk33+nohJlOpiCxOMK2X6+snbqYVInKb25YJ21e8XJmwfcXLlQnbV7xcvmxfIvK4iGwTkeUxbQm3KRH5lbtO1ojIuQnmmdQ2eQRVzakHzsVuG4HG7vsXgeuBB4DRbtto4E9xpg0C64EuQANgCdDb41wjgJDb9qd4udzPSoE2aVxfdwM/r2HatK+vKt95EPhNmtdXX2A5cBTOMblpQPcM2L4S5fJ7+0qUy+/tK24uv7Yv4CzgZGB5TFvcbQro7a6LhkBndx0F48yzxm0y3iMnewo4/8iNRSSE84/+Kc5wGU+5nz8FXBpnOq+H2Tgil6q+raph9/MPcK7NSLd46ysZaV9fFR+IiABXAONSuLxknAB8oKp73X+3WcBl+L99xc2VAdtXovWVjLSvr4oP0719qeps4PMqzYm2qUuA51X1gKpuBNbhrKuqktkmj5BzRUFVNwN/AT4GtgC7VPVtoJ2qbnG/swVoG2dyz4bZqCZXrO8BkxPNAnhbRBaIc1V3StSQ6xZ3t8PjCbqefq6vwcBWVV2baBZ4sL5w/ro8S0SOFpGjgPNxLrj0dfuqJlestG9fNeTybfuqIRf4t33FSrRNJbtektkmj5BzRcHduC7B6VYdCzQRkf9JdvI4bSk5PaumXCJyJxAGnk0wi0GqejIwErhZRM7yONfDQFdgAM4v5QfjTR6nLS3rC7ia6v+K82R9qeoqnN0wU4EpON34cLUTHeLZ+qopl1/bVzW5fN2+kvh39GX7SpJn6wVysCgAw4GNqlqmquXAq8AZwFYRaQ/gPm+LM22Nw2x4kAsRuQ64ELhG3R2AVanqp+7zNuA14ncXU5ZLVbeqakRVo8DYBMvza32FgMuBFxJN7OH6QlUfU9WTVfUsnC7/WvzfvhLl8nv7ipsrA7av6taXr9tXjETbVLLrJZlt8gi5WBQ+Bk4TkaPc/YLDgFU4w2Vc537nOmB8nGm9HGYjbi4ROQ+4HbhYVffGm1BEmohIs4rXOAcPl8f7bgpztY/5zmUJlpf29eV+NhxYraqb4k3o8fpCRNq6z51wfnmMw//tK26uDNi+EuXye/tK9O8IPm9fMRJtUxOAq0SkoYh0xjlwP68W01evPkfMM/UB/A5YjfMP9TTOUfqjgek4fw1MB1q73z0WmBQz7fnAhzhH9O9MQ651OPsHF7uPR6rmwjn7Yon7WJGmXE8Dy4Cl7sbVPhPWl9v+JHBjle+mc329C6x05z/MbcuE7SterkzYvuLlyoTt64hcfm1fOAVpC1CO0xO4IdE25X7/TnedrAFGxrT/Gyiubpus6WFXNBtjjKmUi7uPjDHG1JEVBWOMMZWsKBhjjKlkRcEYY0wlKwrGGGMqWVEwvhKRPWlYxj0iMryO0w6QBCN4ZiIReVlEurivm4rIwyKyXkQWuUMy/KCG6WdWHXVTRG4TkX+JSKGITPEyv/GfFQWT00QkqKq/UdVpdZzFAJxz5TOeiPTBGS1zg9v0b2AnzuifJwHnAa1rmM04nIvEYl0FjFPVMmCLiAxKYWyTYawomIwgIkPcv1JfFmf8/2fFMVJEXqzyvTfc1w+LSIk44+H/LuY7pSLyGxF5D/iWiDwpIt90P/uNiMwXZxz9Me7V0hV/If9JROaJyIciMti9ivYe4Epxxs2/skrm60XkdRF5Q0Q2isgtIvK/7l/lH4hIa/d7P3CXuUREXhFnADZE5FtujiUiMttt6+NmWCzOYHHd3fbX3b/0V0jiAdiuwb1qVUS64gy98Gt1hpJAnSFD/hST/xdurqUx6+9l4EIRaeh+pwjnwq333M9fd5djclUqrxC0hz1q+wD2uM9DgF0447gEgDnAmTjDZ38MNHG/9zDwP+7riquGg8BMoJ/7vhT4ZcwyngS+GTuN+/pp4CL39UzgQff1+cA09/X1wD8SZL8e54rhZkChm/9G97OHgNvc10fHTPMH4Mfu62VAB/d1S/f57zhjFIFzD4HGVX7WxjhXeB8dJ88s4ET39cXAa9Ws9xE49/gVd32/CZzlfjYRuMR9PRr4c8x0HYBlfm839vDuYT0Fk0nmqeomdf6yXQwUqTPW/RTgInEGKruAQ2O4XCEiC4FFQB+cm49USDSY2VARmSsiy4Cz3ekqvOo+LwCKksw8Q1W/VGfXyi7gDbd9Wcw8+orIu+4yr4lZ5vvAk+5+/qDbNge4Q0RuB45X1X1u+60isgTnnggdcca7qao9UBYvpIjc6fY+KgZOG+E+FgELgV4x84zdhXQVh48Wug2n52ByVMjvAMbEOBDzOsKh7fMF4GackSznq+qX7kBgPwe+pqo7ReRJoFHM9F9VnbmINAL+hTM2zCcicneVaSqWH7vs2mSOxryPxszjSeBSVV0iItfj9IpQ1RtF5FScQrdYRAao6nMiMtdte0tEvu/OazhwuqruFZGZVXJX2BfTvhLoLyIBVY2q6r3AvTEH9gW4T1UfjTOf14G/isjJOD2VhTGfNXKXY3KU9RRMNpiJc6vCH3CoB9Ac5xf/LhFphzOufU0qfmFuF5GmwDeTmOZLnN1D9dEM5wBtATH740Wkq6rOVdXfANuBju6ZQxtU9W84A8X1A1oAO92C0As4LcFyVgHdAFR1HVAC/EFEgu7yGnFoLP63gO+56wER6SDuqKGqugdnnT/OkfcU6IE3I4KaDGFFwWQ8VY3g7PMe6T6jqktwdn2swPnl9X4S8/kCZ+z+ZTh/Dc9PYvEzgN7xDjTXwl3AXJwbuqyOaf+zODd/Xw7Mxhl180pguTg3jO8F/Adn91lIRJYCv8fZhRTPRNxeiOv7OCNlrhORBTj3Ib4dQJ272D0HzHF3a73M4cVvHNAf5xaYsYa6yzE5ykZJNSZHiEhjnCI2yC2kXixjNs5B6J1ezN/4z4qCMTlEnAvPVqnqxx7MuxCn4Lye6nmbzGFFwRhjTCU7pmCMMaaSFQVjjDGVrCgYY4ypZEXBGGNMJSsKxhhjKllRMMYYU+n/A2d5rp+pI/paAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from scipy.optimize import curve_fit\n", "\n", "# Vi får staplarnas mitt genom att beräkna medelvärdet av deras gränser\n", "medel = (granser[:-1] + granser[1:])/2\n", "\n", "# För att göra en anpassning behövs någon sorts startvärden.\n", "# Vi gissar rimliga värden på a, mu, sigma, b och c i den ordningen.\n", "p0 = [100, 90, 1, 1, 1]\n", "\n", "# Vi beräknar de optimala koefficienterna \n", "koefficienter, _ = curve_fit(gauss, medel, handelser, p0=p0)\n", "\n", "# Vi beräknar en anpassning med hjälp av koefficienterna\n", "anpassning = gauss(medel, *koefficienter)\n", "\n", "# Vi ritar ut histogramet tillsammans med den anpassade funktionen.\n", "plt.hist(data['M'], bins=antal_staplar, range=(undre_grans, ovre_grans))\n", "plt.plot(medel, anpassning, label='anpassning')\n", "plt.xlabel('Invariant massa (GeV)')\n", "plt.ylabel('Antal händelser')\n", "\n", "# Vi skriver ut väntevärdet och pikens bredd (FWHM, https://sv.wikipedia.org/wiki/Halvv%C3%A4rdesbredd)\n", "print('Väntevärde = ', koefficienter[1])\n", "print('Pikbredd = ', 2*np.sqrt(2*np.log(2))*koefficienter[2])\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4627b51c", "metadata": {}, "source": [ "Med dessa resultat som grund kan vi uppskatta Z-bosonens massa till 90,86 GeV. Du kan testa själv: Hur ser anpassningen ut om du tar bort den linjära korrigeringen av bakgrundsstrålningen, eller om du använder [Breit-Wigner -fördelningen](https://en.wikipedia.org/wiki/Relativistic_Breit%E2%80%93Wigner_distribution) istället? Ger Breit-Wigners en bättre anpassning än Gauss' fördelning? Du kan också testa hur antalet staplar kan påverka resultatet." ] } ], "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.8" } }, "nbformat": 4, "nbformat_minor": 5 }