commit 59c077e76ee49c25d952855f9bd7757f33b07aa6
parent 7577c40224a027fea6622a7a99f3f2fc0f73808b
Author: miksa <milutin@popovic.xyz>
Date: Fri, 9 Apr 2021 12:06:13 +0200
t0 method not working
Diffstat:
8 files changed, 122 insertions(+), 210 deletions(-)
diff --git a/sesh1/prog/Untitled Folder/Untitled.ipynb b/sesh1/prog/Untitled Folder/Untitled.ipynb
@@ -2,83 +2,158 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 52,
- "id": "automatic-recommendation",
+ "execution_count": 21,
+ "id": "emotional-newman",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
- "from sympy import *\n",
+ "import sympy as sp\n",
+ "import matplotlib.pyplot as plt\n",
+ "from scipy.optimize import curve_fit\n",
+ "from ipywidgets import interact\n",
"global m_p; m_p = 0.13957\n",
- "from sympy.utilities.lambdify import lambdify\n",
- "\n"
+ "\n",
+ "sig_p = lambda x: np.sqrt(1 - 4*m_p**2/x)\n",
+ "g_s = lambda s, m_q, g_q: g_q*s/m_q**2 * (sig_p(s)/sig_p(m_q**2))**3 * np.heaviside(s - 4*m_p**2, 0)\n",
+ "\n",
+ "def model(s, m_q, g_q, m_w, g_w, e_w, a, b, c):\n",
+ " part1 = (m_q)**4/((m_q**2 - s)**2 + m_q**2*g_s(s, m_q, g_q)**2)\n",
+ " part2 = 1 + (e_w * 2*s * (m_w**2 - s))/((m_w**2 - s)**2 + m_w**2*g_w**2)\n",
+ " part3 = c*(1 + a*s + b*s**2)**2\n",
+ " return part1 * part2 * part3"
]
},
{
"cell_type": "code",
- "execution_count": 118,
- "id": "demographic-distributor",
+ "execution_count": null,
+ "id": "fitting-consciousness",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "id": "western-clearing",
"metadata": {},
"outputs": [],
"source": [
- "#m_q, g_q, m_w, g_w, e_w, a, b, c = symbols(\"m_q g_q m_w g_w e_w a b c\")\n",
- "var = Matrix(symbols(\"m_q g_q m_w g_w e_w a b c\"))\n",
+ "data = np.loadtxt('../data/SND-VFF.txt')\n",
+ "x_data = data[:,0]\n",
+ "y_data = data[:, 1]\n",
+ "sigma = data[:, 2]\n",
"p0 = [0.765 ,0.115 ,0.813 ,0.041 ,-0.006 ,-1.005 ,1.014 ,1.854] # guess"
]
},
{
"cell_type": "code",
- "execution_count": 99,
- "id": "stock-julian",
+ "execution_count": 75,
+ "id": "confidential-mapping",
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[<matplotlib.lines.Line2D at 0x7f25c8f59460>]"
+ ]
+ },
+ "execution_count": 75,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAGbCAYAAAALJa6vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABIiElEQVR4nO3deVic1cH+8e8Zhm3YEhLIDpM9JjExiiZu0bovdWlrtYo1WhW1e7VvN/x1saVv375va2s3pdaldjTu+xZ345KF7KtZgZCEPRBg2Of8/oDExCQGEphnlvtzXblkhoG5H4Th5jznOcdYaxERERGRnnM5HUBEREQk3KhAiYiIiPSSCpSIiIhIL6lAiYiIiPSSCpSIiIhIL7mD+WSDBw+2Xq83mE8pIiIickSWLFlSba3NONj7glqgvF4vRUVFwXxKERERkSNijCk51Pt0Ck9ERESkl1SgRERERHpJBUpERESkl1SgRERERHpJBUpERESkl1SgRERERHpJBUpERESkl1SgRERERHpJBUpERESkl1SgRERERHpJBUpERESkl1SgRERERHpJBUpERESkl1SgRESOks/nw+v14nK58Hq9+Hw+pyOJSD9zOx1ARCSc+Xw+8vLy8Pv9AJSUlJCXlwdAbm6uk9FEpB9pBEpE5Cjk5+fvLU97+P1+5syZoxEpkQimESgRkaNQWlp60Ps7OzsBjUiJRCqNQImIHIWsrKy9b8ePmMzAc25h+C33kzT17L33+/1+8vPznYgnIv1EBUpE5CgUFBTg8XhIP//bDL329yRPO4+YpAEkjjlhv8cdaqRKRMKTTuGJiByF3NxcPq6O46WdHnYvfo6Gj+Yy8OLbiR00ar/H7TtSJSLhTyNQIiJH4ePNNbxakcQ5x2RS/UYhD91/L9TvJDZ9BJiul1iPx0NBQYHDSUWkL2kESkTkCG2qbOCbviWMHpzE3Vcdh8tlyM3NZUltHM9sjyN2wFCGp8ZSUFCgCeQiEUYjUCIiR6C4uolr/rmQGJeL+6/LISUhdu/7rr3sHABe+WApxcXFKk8iEUgFSkSkl8p2+cm9fyEdAcujN8/EOzhpv/ePy0wGYFNVoxPxRCQIVKBERHphSUktF/zhLbZX1rDmb7dx3sxjD1goMzUhlsyUeDZVqkCJRCrNgRIR6aG5i0rJf3YlbbsqqXj617RXl1ACB10oc1xmMhtVoEQilkagREQOY3dLO9+fu4yfPLOKzh3r2PHQ92ivLtn7/oMtlDkuM5nNlY1Ya4MdV0SCQCNQIiKfY0lJLd+bu5yd9S3cfu4Evn/eZdhA5wGP++xCmeMyk2ls7aBidytD0xKCFVdEgkQFSkTkIDo6A+T9+VneLo+jY3cl7sWPMGjKzWSNGklJSckBj//sQpnjMronklc2qkCJRCCdwhMR6ebz+fB6vcSmZjD6pnt4uzKBxjXvsOPB71Cy9D3y8vK46KKL8Hg8+33cwRbK3HslXmVD0PKLSPCoQImI0FWe8vLy2NmeyNDr7sYMGEHVC7+n5pW7sW3NQNdcp1deeYXCwkKys7MxxpCdnU1hYeEBaz1lpMSTkuDWUgYiEcoEc4JjTk6OLSoqCtrziYj0lNfrpTppNIMu/B4dDVVUPf1r2mu2HfA4YwyBQKBHn/NLf/+QeLeLuXkn93VcEQkCY8wSa23Owd6nESgREaBm4DEMvuSHtG5fS/m/7zhoeYLebQo8LiOZTZVNfRVRREKICpSIRL1HFpQw6Pxv49+8mIonf0mg5eDzlnq7KfDYzGSqG1vZ3dLeV1FFJESoQIlIVHtqSRn/77nVTExpp+m1u6Hz07ITGxvLoEGDPneu0+dJT4oDYHezCpRIpFGBEpGotWhrLT99ZiWnjhvEiz++lMJ7/77f5PAHH3yQ6upqAoHAEW0K7ImLAcDfduC6USIS3rQOlIhEpZKaJm55pIhR6R7+fs0JxLld5Obm9rokfZ6kuK6X2KbWjj77nCISGjQCJSJRp7mtk5seLsICD8w5kTRPbL88z54RqGaNQIlEHBUoEYk6v31lHRsrG/nL1TPwDk7qt+d57+03ADj3oi/i9Xrx+Xz99lwiElwqUCISVd5aV8EjC0q4+fTRnD4+o9+ex+fz8Ztf3Nl1wx1PSUkJeXl5KlEiEUIFSkQi3p4tWtzJ6dx43zsMTejkh+dP7NfnzM/Px99QB4ArLhHoWsk8Pz+/X59XRIJDBUpEItqeLVpKSkoYeNZNWHc8a//1Xzz1+Nx+fd7S0lICbS0AmNiE/e4XkfCnAiUiES0/Px+/309C9nSSJp9B/YInaSjb0O8jQVlZWdj2rgLl2qdA9WYlcxEJXSpQIhLRSktLIcZN+rm30r5rB/ULnvr0/n5UUFCAJyEe29GOiesqUL1dyVxEQpfWgRKRiJaVlcWuYScRO2gUFU/8fO9K4/09ErRnPan8olZcsQlkZ2dTUFDQp+tMiYhzNAIlIhHtx7/8LQNO+RpNn3xIy9alQPBGgnJzcxmeOZgb8m47opXMRSR0qUCJSEQrSzuWmNg4kja8fsR72h0NT7wbf5tWIheJNDqFJyIRa2t1E3MXb+Pak73c9bsljmRIiovRXngiEUgjUCISsf74xgbi3S6+c9Z4xzIkxsXgb1WBEok0KlAiEpFWb6/nxRU7uPG00WSkxDuWIynOjb9dp/BEIo0KlIhEpD+9uZG0xFhunj3G0RyeeLdGoEQiUI8LlDEmxhizzBjzUvft0caYhcaYTcaYx40xcf0XU0Sk5z4pb+DNdRXccKqX1IRYR7N4YmNo0iRykYjTmxGo7wHr9rn9P8Dd1tpxwC7gxr4MJiJypO59bzOeuBiuP8XrdBQ88ZpELhKJelSgjDEjgYuB+7tvG+As4KnuhzwMXN4P+UREemVbrZ8XVuzgmpOyGOBxfmDc030VnrXW6Sgi0od6OgL1J+BHQKD79iCgzlq7Z1y6DBhxsA80xuQZY4qMMUVVVVVHk1VE5LDue38zMcZw0+nOzn3awxPnpjNgaesMHP7BIhI2DlugjDFfBCqttUe0iIq1ttBam2OtzcnIyDiSTyEi0iPVja08UVTGV04YwdC0hMN/QBAkxcUAaCK5SITpyUKapwKXGmMuAhKAVODPwABjjLt7FGoksL3/YoqIHN6jC0tp6wiEzOgTdI1AATS1dTAwyflTiiLSNw47AmWt/am1dqS11gt8DXjbWpsLvANc0f2wOcDz/ZZSROQw2joC/GdBCWdMyGBsRrLTcfbyxHeNQDVrIrlIRDmadaB+DNxujNlE15yof/VNJBGR3nt19U4qG1q5/lSv01H24+k+hdekAiUSUXpVoKy171prv9j99hZr7UnW2nHW2q9aa1v7J6KIyOE99FExowcnccb40Jpr+dF77wBw8uln4PV68fl8DicSkb6glchFJOwt31bHstI65pycjctlnI6zl8/n43//+zcAGHcCJSUl5OXlqUSJRAAVKBEJe//+qJjkeDdfOWGk01H2k5+fj7+hHgAT13VVoN/vJz8/38lYItIHVKBEJKzV+dt4adVOvjRjBCkOb9vyWaWlpdj2ZgBcsQn73S8i4U0FSkTC2rPLttPWEeBrJ41yOsoBsrKysG0twKcjUHvuF5HwpgIlImHLWsvcRduYPjKNKcPTnI5zgIKCAhJiu15mTfcIlMfjoaCgwMlYItIHerKQpohISFpaWscnFQ3895ePdTrKQeXm5gKQv7wdV1wi2dnZFBQU7L1fRMKXRqBEJGzNXVRKUlwMl0wf7nSUQ8rNzSUt2cN3f/BDiouLVZ5EIoQKlIiEpd0t7by0cieXHjeC5PjQHkxPiovBr4U0RSKKCpSIhKWXV+6kub2Tq04Mvcnjn5WoAiUScVSgRCQsPb2kjHGZyUwfGXqTxz8rKd5NU1uH0zFEpA+pQIlI2CmubqKoZBdfOX4kxoTOyuOH4tEIlEjEUYESkbDzzNIyXAa+NGOE01F6xBPnxq8RKJGIogIlImElELA8vXQ7p44bzNC0hMN/QAjQCJRI5FGBEpGwsqi4lu11zVwRYvvefZ6kODf+VhUokUiiAiUiYeXpJWUkx7s5b/JQp6P0WGJcjCaRi0QYFSgRCRst7Z28trqcC6YOJTEuxuk4PZYUH0NzWyfWWqejiEgfUYESkbDx3oYqGlo7uDSEVx4/GE+cm46Apa0z4HQUEekjKlAiEjZeWrmT9KQ4Thk7yOkovbJ25TIAklIH4vV68fl8DicSkaOlAiUiYcHf1sGbayu46NihuGPC56XL5/Mx95GHu264EygpKSEvL08lSiTMhc+rkIhEtbfWVdLc3skl08Lr9F1+fj4tTbsBcMV1Lbvg9/vJz893MpaIHCUVKBEJCy+u2MGQ1HhO9KY7HaVXSktLsW3NAJjYhP3uF5HwpQIlIiFvd0s7735SxRenDcflCv2tW/aVlZVFoL0VABOXsN/9IhK+VKBEJOTNW1NBW2eAS8Ls6juAgoIC4mO6li9wdY9AeTweCgoKnIwlIkfJ7XQAEZHDeXHFDkalJzJ9ZJrTUXotNzeXqhYXf9oIrrhEsrOzKSgoIDc31+loInIUNAIlIiGttqmNDzZV88VpwzEmvE7f7XHlVy4H4P6HHqG4uFjlSSQCqECJSEh7dfVOOgM27K6+21dSXNdgvzYUFokcKlAiEtJeXLGDsRlJHDMsxekoR2zPtjN+7YcnEjFUoEQkZFXsbmHh1loumR6+p+8A4twuYmOMRqBEIogKlIiErJdX7sRa+GIYn77bwxPnVoESiSAqUCISsl5cuYPJw1IZl5nsdJSj5omLoalVp/BEIoUKlIiEpIrdLSwrreOiY4c6HaVPJMVrBEokkqhAiUhIemNtBQDnT4mcAtWgESiRiKECJSIhad7aCkYPToqI03cAKfFuncITiSAqUCIScna3tPPx5mrOmzwkrK++21dSfAyNLSpQIpFCBUpEQs476ytp77ScN2WI01H6TFK8m0aNQIlEDBUoEQk589ZWMDg5nuNGDXQ6Sp9JiXfTpIU0RSKGCpSIhJTWjk7eXV/JuZMziXFFxuk76B6BaunAWut0FBHpAypQIhJSPtpcQ1NbJ+dFyNV3eyTFu+kIWFo7Ak5HEZE+oAIlIiFl3poKkuJiOGXsIKej9KmUhK4NhXUlnkhkUIESkZARCFjeWFvBmZMyiXfHOB2nTyXFdRUoTSQXiQwqUCISMpZt20V1YyvnTY6cq+/2SIpXgRKJJCpQIhIy5q2pIDbG8IVJmU5H6XN7TuFpLSiRyKACJSIhwVrL62vKmTVmEKkJsU7H6XPz33kTgC+cdyFerxefz+dwIhE5GipQIhISNlU2Ulzjj7ir7wB8Ph+//dXPATBxiZSUlJCXl6cSJRLGVKBEJCTM6948OBLnP+Xn5+OvrwXAFZcIgN/vJz8/38lYInIU3E4HEBEBmLemnONGDWBIaoLTUfpcaWkpxHYdl4nz7H+/iIQljUCJiON21jezoqw+ova+21dWVha2rQX4dARqz/0iEp5UoETEcW/sPX0XefOfAAoKCvB4Egm0+nHFd41AeTweCgoKHE4mIkdKp/BExHHz1lQwJiOJcZnJTkfpF7m5uQDkL2rDFZdIdnY2BQUFe+8XkfCjESgRcVS9v50FW2oidvRpj9zcXMZkjeCqr19PcXGxypNImFOBEhFHvfNJJR0By/kROv9pX8nxbu2FJxIhVKBExFHz1paTmRLP9JEDnI7S75LiVKBEIoUKlIg4pqW9k3c/qeLcyUNwuYzTcfpdcoKbBm3lIhIRVKBExDEfba7G39YZkauPH0xyvJumNhUokUigAiUiQefz+fB6vVx5+2+hrZnihfOcjhQUSfEx2kxYJEKoQIlIUPl8PvLy8igp3Ubi2JNo2rSIb94aHfvCJcfH0tTa6XQMEekDKlAiElT5+fn4/X7ih08iJmkA/o0LomZfuOT4GNo6A7R2qESJhDsVKBEJqj37vyVOmIXtaKd5S9F+90ey5PiutYs1CiUS/lSgRCSo9uz/5hk/i5aSFdi25v3uj2RJewuU5kGJhDsVKBEJqoKCAlJGTiR24HD8Gz8GomdfuD0jUI0qUCJhT3vhiUhQ5ebm8k5lPG+WB2jevCiq9oVLTlCBEokUKlAiEnTl7iEcn+2iuKHW6ShBlaQRKJGIoVN4IhJU2+uaWb19d9QsnrmvvafwtBaUSNhTgRKRoHpjTTkA502O/M2DPytZk8hFIoYKlIgE1by1FYzPTGZMRrLTUYJOp/BEIocKlIgETZ2/jYVbazlvSvSNPgEkxcUAKlAikUAFSkSC5u31lXQGLOdNjr75TwCPz30MOlq567e/x+v1RsX2NSKRSgVKRIJm3poKhqYmcOyINKejBN2ePQA7W/yYuARKSkrIy4uOPQBFIpEKlIgERUt7J+9tqOLcyUNwuYzTcYJuzx6AgbZmXHEegKjZA1AkEqlAiUhQzN9YTXN7Z9TOf9qz11+gzY+JSzzgfhEJLypQIhIU89aUk5LgZtaYQU5HccSevf5sWzOueM8B94tIeFGBEpF+19EZ4M11FZw9KZPYmOh82SkoKMDj8RBo9ePqHoGKlj0ARSLRYV/JjDEJxphFxpgVxpg1xphfdd8/2hiz0BizyRjzuDEmrv/jikg4WlKyi13+9qhcfXyP3NxcCgsL8cS6MHGJZGdnU1hYGBV7AIpEop7shdcKnGWtbTTGxAIfGGNeBW4H7rbWzjXG3AvcCPyjH7OKSJiat7aCOLeL2RMynI7iqNzcXNZ4VvH66nKWFBc7HUdEjsJhR6Bsl8bum7Hd/yxwFvBU9/0PA5f3R0ARCW/WWuatLee0cYP3bmUSzVLi3VpIUyQC9GgygjEmxhizHKgE3gA2A3XW2j2vAmXAiEN8bJ4xpsgYU1RVVdUHkUUknKzb2cC22uao3PvuYJLi3bR2BGjvDDgdRUSOQo8KlLW201p7HDASOAmY1NMnsNYWWmtzrLU5GRnRPXwvEo3mrS3HGDj7GBUo+HQ/PG0oLBLeenU5jLW2DngHOBkYYIzZMx4/Etjet9FEJBLMW1NBTvZAMlLinY4SElK0obBIROjJVXgZxpgB3W8nAucC6+gqUld0P2wO8Hw/ZRSRMLWt1s/anbujdu+7g0lO6CpQdf52h5OIyNHoyQjUMOAdY8xKYDHwhrX2JeDHwO3GmE3AIOBf/RdTRMLRG2srADhX85/2mjQ0BYCVZfUOJxGRo3HYS2KstSuBGQe5fwtd86FERA5q3tpyJg5JwTs4yekoIWP04CQGJ8exuLiWa2ZqFXKRcBWdSwKLSL/b1dTGoq21Ubv33aE8+uijVK9byJPvLsXr9eLz+ZyOJCJHQAVKRPrFm+sqCFg0/2kfPp+PvLw8aj9ZjDttCGW1TeTl5alEiYQhFSgR6Rfz1lYwPC2BqSNSnY4SMvLz8/H7/bSWrQEgfuQU/H4/+fn5+z3O5/Ph9XpxuVwapRIJUSpQItLnmts6mb+xivOmDMUY43SckFFaWgpAW+VWAq1+EkZN2e9++HSUqqSkBGstJSUlGqUSCUEqUCLS597fWEVLe0Crj39GVlb3pHEboHX7WuJHdhWo9PT0vY/ZM0q1r4ONUomIs1SgRKTPzVtTQVpiLCeOTj/8g6NIQUEBsbGxALRsW0NcRjauhBQaGhr2jjDtOxq1r0PdLyLOUIESkT7V3hngzXUVnD0pk9gYvcTsKzc3l9TUrjlhrWVrAYgfOZm2tjbmfONGEkdNJW3WV0k96cvEJO9fPveOXolISNDW6CLSpxZtraW+uZ3zp+rqu4Opra0FoHXnBmxHO2knX0nqSV8mfth4jDtu7+MGnDGH5s1F7Hr7fuLa6ikoKHAqsogchP48FJE+9drqchJiXcwer83DD2bvSFJnOy0lK4gbMhYTE8PupS9R+fSv2XbPNWy/72YaFj1D/MjJjLjuf/nNPfeTm5urq/NEQoix1gbtyXJycmxRUVHQnk9EgisQsJz8u7eYMWog9379BKfjhKQ9V9n5/X5wxWBcbmxH6wGPM8bwyc56rrl/IYGA5Zph1fz8ezftN8Hc4/FQWFhIbm5uMA9BJGoYY5ZYa3MO9j6NQIlIn1lRVkfF7lbOn6qr7w4lNzeXwsJCsrOzMTaAy3Yc9HFZWVmMH5LC43mziI1x8ZfVhhYTt99jdHWeiHNUoESkz7y2phy3y3DWRBWoz5Obm0txcTGBQICHH34Yj8ez3/s9Hs/eOU9jMpJ55MaTCMTEMfji24H919XS1XkizlCBEpE+Ya3l9dXlnDx2EGmeWKfjhI39RqSMITs7+4DTcuOHpOBa/gyJo48nJefS/T5eV+eJOEMFSkT6xIaKRopr/Jw/RVff9da+I1LFxcUHndP0mxsupHXLYgaecT2xGaOB/UeqRCS4VKBEpE+8vqYcY9Dq4/3k2mtz+dXFEzDtzQw6/1tkZ3s1gVzEQVoHSkT6xGuryzk+ayCZqQlOR4lYeXOuIX1qGT98cgX/9/R8rjhhpNORRKKWRqBE5Khtq/Wzduduzp+i0af+9uUZI5iRNYDfvbqe3S3tTscRiVoqUCJy1F5fUw6g+U9B4HIZfnXpFGqaWrnnzY1OxxGJWipQInLUXl9TzqShKWQPSnI6SlSYNnIAV+WM4qGPiimubnI6jkhUUoESkaNS1dBKUckujT4F2e3nTsAdY7j7zQ1ORxGJSipQInJU3lhbgbVwgTYPDqrM1ARuOHU0L6zYwbqdu52OIxJ1VKBE5Ki8vqacrHQPk4amOB0l6tw6eyzJ8W7+MO8Tp6OIRB0VKBE5YnX+Nj7cVM2FU4dijDn8B0ifSvPEcusZY3lzXSVLSnY5HUckqqhAicgRm7e2go6A5aJjhzkdJWrdcKqXwclx/ElzoUSCSgVKRI7Yq6t2MmJAItNGpjkdJWp54tzcdPoY5m+sZsW2OqfjiEQNFSgROSL1ze18sKmai47V6TunXTsrm7TEWP76ziano4hEDRUoETkib6ytoL1Tp+9CQXK8m+tP8fLG2go+KW9wOo5IVFCBEpEj8kr36bvjRg1wOorQNRcqKS6Gv7+rUSiRYFCBEpFeq29uZ/7GKl19F0IGeOK4dlY2L67YodXJRYJABUpEeu2tdV2n7y7U6buQcuNpo3HHuLj3vc1ORxGJeCpQItJrr6zaybC0BGbo9F1IyUxN4KqcUTy9tIwddc1OxxGJaCpQItIru1vaeX9DNRdOHYbLpdN3oeaWM8ZgLRS+v8XpKCIRTQVKRHrlrXUVtHUGuHia9r4LRSMHepiW1sJD72/AnTwQr9eLz+dzOpZIxFGBEpFeeWVVOUNTE5gxaqDTUeQgfD4fb97zI2yMm5QTLqOkpIS8vDyVKJE+pgIlIj3W0NLOexuquGDqUJ2+C1H5+fk07NiM/5MPSTn+IkxcIn6/n/z8fKejiUQUFSgR6bG311fS1hHg4mm6+i5UlZaWArB74TO44pNIOe6C/e4Xkb6hAiUiPfbyyp0MSY3nhCydvgtVWVlZALSVb6SlZAUpOZeBy733fhHpGypQItIjDS3tvLuhSlffhbiCggI8Hg8A9Qufxp0ymIHHnUtBQYHDyUQii9vpACISHl5fU0FbR4BLpg93Oop8jtzcXKBrLlRp8TKo2864L97K1Vdf4nAykciiESgR6ZHnl29n5MBEjs8a4HQUOYzc3FyKi4sJBALcfctFVLbG8O6GSqdjiUQUFSgROayqhlY+2lzDpdOHa++7MPPFacMZnpbAve9pYU2RvqQCJSKH9cqqnXQGLJcdN8LpKNJLsTEubjx9DIu21rK0dJfTcUQihgqUiBzWCyt2MHFIChOHpjgdRY7A104cRVpiLIUahRLpMypQIvK5ttX6WVKyi0uP0+TxcJUU7+brs7J5fW05W6oanY4jEhFUoETkc724cgcAl+rqu7A25xQvsTEu/jl/q9NRRCKCCpSIfK4Xlu/g+KwBjEr3OB1FjkJGSjxXnDCSp5eWUdnQ4nQckbCnAiUih7ShooH15Q0afYoQN58+hvbOAA9/VOx0FJGwpwIlIof0wvIduAxcPE0FKhKMHpzEBVOG8sjHJTS2djgdRySsqUCJyEFZa3lhxQ5OHTeYjJR4p+NIH8mbPYbdLR3MXaTNhUWOhgqUiBzU8m11lNb6dfouwszIGsjM0en864OttHcGnI4jErZUoETkoJ5fvoM4t4vzpw51Oor0sVvPGMvO+hZeXLHD6SgiYUsFSkQO0NEZ4OVVOzlrYiapCbFOx5E+dubEDCYOSeG+97ZgrXU6jkhYUoESkQPM31RNVUMrl8/Q6btIZIwhb/YYPqlo4N0NVU7HEQlLKlAicoCnl5Qx0BPLWZOGOB1F+skl04czLC2B+97b7HQUkbCkAiUi+6lvbmfe2gounT6cOLdeIiJVnNvFjaeNZsGWWpZvq3M6jkjY0aujiOznpZU7aOsI8JUTRjodRfrZ107KIsFlufzHf8HlcuH1evH5fE7HEgkLKlAiAoDP58Pr9fLDvz4J9TtZ+e7LTkeSfvb8U49Tu+g57KjjiEkbSklJCXl5eSpRIj2gAiUi+Hw+8vLy2L67nfgRx7Br6avccot+kUa6/Px8ahY8DZ2dpJ54OQB+v5/8/Hxng4mEARUoESE/Px+/30/y1LOxgU6a1ryjX6RRoLS0lEBTHY2r3yLp2HNwedL23i8in08FSkS6f2EakqZ+gZbiZXQ27drnfolUWVlZAOxe/CzGHUvK8V/c734ROTQVKBEhKyuLhOxpuFMzaVz99n73S+QqKCjA4/HQUbud5g0LSDn+i3jSBlJQUOB0NJGQpwIlIhQUFJB23HkEWhpp3rgAAI/Ho1+kES43N5fCwkKys7PZvegZYhJTuP5X95Gbm+t0NJGQpwIlIlx2xVWkHDMbV9ly6GwnOzubwsJC/SKNArm5uRQXF9OyfR0nedNZ3jxImwyL9IAKlIjwyqqdtFvD03/4LwKBAMXFxSpPUeiWM8awva6ZV1btdDqKSMhTgRIRnl5SxujBSRyfNdDpKOKgL0zMZHxmMvdqk2GRw1KBEolyJTVNLNxay5dnjMAY43QccZDLZbh59hjW7dzN/I3VTscRCWkqUCJRbu7ibbgMfDVnlNNRJARcdtxwhqTGc682GRb5XCpQIlGsrSPAk0XbOGvSEIamJTgdR0JAvDuGG08bzUeba1hWusvpOCIhSwVKJIq9ta6C6sY2rpmp0Sf5VO7MbAZ6YvnL25ucjiISslSgRKLYo4tKGZ6WwBkTMp2OIiEkKd7NTaeP4e31lawqq3c6jkhIUoESiVKlNX7mb6zmqhOziHFp8rjs77qTs0lNcPOXtzc6HUUkJB22QBljRhlj3jHGrDXGrDHGfK/7/nRjzBvGmI3d/9X1zyJhZO7iUlwGrjxxpNNRJASlJMRyw6mjmbe2gnU7dzsdRyTk9GQEqgO4w1o7GZgFfMsYMxn4CfCWtXY88Fb3bREJA+2dAZ4oKuOsSZkMS0t0Oo6EqG+cOprkeDd/fUdzoUQ+67AFylq701q7tPvtBmAdMAK4DHi4+2EPA5f3U0YR6WNvrq2gurGVq0/SZsFyaGmeWOacks0rq3ayqbLB6TgiIaVXc6CMMV5gBrAQGGKt3bPefzkw5BAfk2eMKTLGFFVVVR1NVhHpI48uKmVYWgJnTtTkcfl8N542hsTYGP6qK/JE9tPjAmWMSQaeBr5vrd3vhLjtWvP/oOv+W2sLrbU51tqcjIyMoworIkdvW+2eyeOjNHlcDis9KY5rZ2XzwoodbK1ucjqOSMjoUYEyxsTSVZ581tpnuu+uMMYM637/MKCyfyKKSF/aM3n8qhO19pP0zM2njyE2xsXfNRdKZK+eXIVngH8B66y1f9znXS8Ac7rfngM83/fxRKQvafK4HImMlHiumZnFM8u2U1KjUSgR6NkI1KnA14GzjDHLu/9dBPwOONcYsxE4p/u2iISw19eUU9WgyePSe7edMZbYGMOf39S6UCIA7sM9wFr7AXCoiRJn920cEelPD31YTFa6R5PHpdcyUxOYc7KXwvlbuO3MsYwfkuJ0JBFHaSVykSixqqyeopJdzDnFq8njckRuOWMsntgY/qRRKBEVKJFo8eBHW0mKi+GrOVp5XI5MelIcN542mpdX7WTNDu2RJ9FNBUokClQ1tPLSip1cccJIUhNinY4jYezG08eQlhjLH+dtcDqKiKNUoESiwKMLS2nrDHDdKV6no0iYS0uMJW/2GN5aX8nS0l1OxxFxjAqUSIRr6wjwn4UlnDkxg7EZyU7HkQhw/SleBiXFaRRKopoKlEiEe2XVTqoaWrleo0/SR5Li3dx25lg+2FTNx5trnI4j4ggVKJEI9+BHxYwZnMTs8dpKSfrOtbOyGZqawB/mfULXbl4i0UUFSiSCLS3dxYptdVx/qheXli6QPpQQG8O3zxpHUcku3l6vnbwk+qhAiUSwhz4sJiXezZeP19IF0veuOnEUYwYn8d+vrqejM+B0HJGgUoESiVAVu1t4ZdVOrjxxFMnxh910QKTXYmNc/OiCSWyqbOSJojKn44gElQqUSIR64MOtBKxlzslep6NIBDt/yhBysgfyxzc20NTa4XQckaBRgRKJQPXN7fgWlHLRscPIGuRxOo5EMGMMP7v4GKobWyl8f4vTcUSCRgVKJAL5FpbQ2NrBrWeMdTqKRIHjswZy8bHDKHx/C5W7W5yOIxIUKlAiEaalvZMHPihm9oQMpo5IczqORIkfXTCRjkCAu9/U4poSHVSgRCLMU0vKqG5s5dYzxjgdRaJI9qAkrp2VzeOLt7GhosHpOCL9TgVKJIJ0dAYofH8L00cN4OQxg5yOI1Hmu2eNJyneze9eXe90FJF+pwIlEkFeXV1Oaa2f284YizFaOFOCa2BSHN88cxxvr6/kw03VTscR6VcqUCIRwlrLP97dzJiMJM6bPMTpOBKlbjjVy6j0RH75whratbimRDAVKJEIMX9jNWt37ubW2WO1bYs4JiE2htkpVWysbGTwrC/j9Xrx+XxOxxLpcypQIhHiH+9uZmhqApfNGO50FIliPp+PP/3wGzRvWULaadewraqevLw8lSiJOCpQIhFgWekuPt5Sw42njSbeHeN0HIli+fn5+P1+at8qxLjjGHjGHPx+P/n5+U5HE+lTKlAiEeDPb21kgCeWq2dmOR1FolxpaSkAHbXb2V30AsnTziVu2IS994tEChUokTC3pGQX735SxS2zx2rTYHFcVtanJb7+o7l0NNSQfs6tZGVlO5hKpO+pQImEubvf2MCgpDiuO1m/oMR5BQUFeDxd+y/atmbq3n2Q+OET+NIdv3c4mUjfUoESCWMLt9TwwaZqbj1jLEkafZIQkJubS2FhIdnZ2RhjGNy0lSxPB+/XDaC+ud3peCJ9RgVKJIzd/eYGMlLiuXaWRp8kdOTm5lJcXEwgEKC4uJi/33gmu/xt/N/rnzgdTaTPqECJhKmPNlWzYEst3zxzLIlxuvJOQtfUEWlcf8po/rOwhCUltU7HEekTKlAiYchayx/f2MCQ1HiuPklX3knou+O8CQxPS+QnT6+irUMrlEv4U4ESCUPvb6ymqGQX3/7COBJiNfokoS8p3s2vL5/CxspG7ntvs9NxRI6aCpRImNkz+jRiQCJXnjjK6TgiPXbWpCFcPG0Yf3l7E5urGp2OI3JUVKBEwszb6ytZsa2Ob581TquOS9j5xSWTSYh18bNnVmGtdTqOyBFTgRIJIx2dAX736nqy0j1cccJIp+OI9FpmSgI/u+gYFm6t5YmibU7HETliKlAiYeTJJWVsrGzkJxdOIjZGP74Snq7MGcVJo9MpeHkdVQ2tTscROSJ6BRYJE02tHfxh3gZOyB7IhVOHOh1H5Ii5XIbffulYWtoD3PmcTuVJeFKBEgkT972/herGVvIvPgZjjNNxRI7KuMxk7jhvAq+vqeC55dudjiPSaypQIiHO5/PhPWY6f3ptNZQuYd38V5yOJNInbjp9DDnZA/n582vYWd/sdByRXlGBEglhPp+PvLw8GkefiXG52P7K38nLy8Pn8zkdTeSoxbgM//fV6XR0Wn701EqdypOwogIlEsLy8/NpT8ok6dhzaFjyEh31Ffj9fvLz852OJtInvIOT+NnFxzB/YzX/WVjqdByRHlOBEglhpaWlDPzCNwi0NFH/8eP73S8SKa6dmcXp4wfz25fXUVLT5HQckR5RgRIJYVknnU/i6OOp/2gugZZPV27OytL+dxI5jDH8/oppuGMMdzyxgs6ATuVJ6FOBEglRbR0BMs6/jc66chqWvbz3fo/HQ0FBgYPJRPresLREfnXpFIpKdlH4/han44gclgqUSIi6/4MtVLXGcMNxyWSPHIExhuzsbAoLC8nNzXU6nkif+9KMEVx07FD+MO8TlpbucjqOyOcywbzqIScnxxYVFQXt+UTC1fa6Zs75w3vMnjCY+76e43QckaCpb27n4nvmA/Dyd08nLTHW4UQSzYwxS6y1B30R1giUSAi668U1APz8kikOJxEJrrTEWP5y9QzK61v4ydNa2kBClwqUSIh5Z30lr6+p4Dtnj2PEgESn44gE3YysgfzX+RN5dXU5Pi1tICFKBUokhLS0d/KLF9YwNiOJm04b43QcEcfcfPoYzpiQwV0vrWXdzt1OxxE5gAqUSAj5x7ubKa318+vLphLn1o+nRC+Xy/CHK6czIDGWbz+6FH9bh9ORRPajV2iREFFS08Q/3tvMJdOHc8q4wU7HEXHc4OR4/nTVcWypbuLOZ1drPpSEFBUokRBgreXO51YTF+PizouPcTqOSMg4Zdxgvn/2BJ5Ztp1/f1zidByRvVSgRELAE0XbmL+xmh9fMJEhqQlOxxEJKd85axznHDOEX7+0loVbapyOIwKoQIk4bmd9M795aR2zxqSTOzPb6TgiIcflMvzxqulkDfLwrUeXsrO+2elIIipQIk6y1vLTZ1bREbD8/ivTcbmM05FEQlJqQiyFX8+hpT3Arf9ZSkt7p9ORJMqpQIk46KklZbz7SRU/vmAiWYM8TscRCWnjMpP5w5XTWbGtjp8/r0nl4iwVKBGHlNe3cNdLaznJm851J3udjiMSFs6fMpTvnjWOJ4rK+I8W2RQHqUCJOMBay8+eXUV7Z4DfXzFNp+5EeuH750zgrEmZ/PKFNby/ocrpOBKlVKBEHPDssu28vb6S/zp/Et7BSU7HEQkrLpfhnqtnMGFICt/0LWV9uVYql+BTgRIJsm21fn7x/Bpysgdy/Slep+OIhKXkeDcPXJ9DUnwM33hwMRW7W5yOJFFGBUokiDo6A3xv7jIA7r7qOGJ06k7kiA1LS+Rfc06krrmdGx9erO1eJKhUoESC6M9vbWRpaR0FXz6WUem66k7kaE0dkcbfrjmetTt2893HltEZ0JV5EhwqUCJBsmBLDX99ZxNXnDCSS6cPdzqOSMT4wqRMfnXpFN5cV8ldL67R8gYSFG6nA4hEg11Nbfzg8eV4ByXxq0unOB1HJOJ8/WQvpbV+/jl/K4OS4/nu2eOdjiQRTgVKpJ9Za/nx0yupbmzlmdtOJSleP3Yi/eGnFx5DbVM7f3xjA6kJbq4/dbTTkSSC6ZVcpJ/5FpYyb20F+Rcdw7Ej05yOIxKxXC7D/3zlWHa3tPPLF9eS5onlSzNGOh1LIpTmQIn0o5Vlddz10lpmT8jgxtP017BIf3PHuPjL1TM4ZewgfvjkSt5YW+F0JIlQKlAi/aS6sZVbH1lCRnI8f7rqOK02LhIkCbExFF6Xw9ThqXzr0aV8tLna6UgSgVSgRPpBR2eA7zy6jOqmNu699gTSk+KcjiQSVZLj3Tx0w0lkp3u46eEiFhfXOh1JIowKlEg/+J/X1vPxlhp++6VjNe9JxCEDk+L4z00zGZqawJwHFrFoq0qU9B0VKJE+9sKKHfxz/lauOzmbK07QBFYRJw1JTWBu3iyGpiVw/YOLWLClxulIEiFUoET60Lqdu/nxUyvJyR7InRdPdjqOiACZ3SVqWFoCNzy4mI83q0TJ0VOBEukjNY2t3PLIElIS3Pw993ji3PrxEgkVmSkJzM07mZEDE7nhoUWaWC5HTa/wIn2gpb2TvEeWUL67hXu/fgKZqQlORxKRz8hIiefRm2eRle7hGw8t5u31WuJAjpwKlMhRCgQsdzy5giUlu/jTVcdxfNZApyOJyCFkpMTz2M2zGJ+Zws3/XsKzy8qcjiRh6rAFyhjzgDGm0hizep/70o0xbxhjNnb/V78xJCr5fD7GXP59Xl65E5Y/y66VbzsdSUQOY1ByPI/ePJOZo9P5weMreOCDrU5HkjDUkxGoh4ALPnPfT4C3rLXjgbe6b4tEFZ/Px3f//DhMPo+GZa9S8vq/yMvLw+fzOR1NRA4jJSGWB64/kQumDOWul9biveTbuFwuvF6vfoalRw5boKy17wOfXTzjMuDh7rcfBi7v21gioS//L/8h+cybaN5SRO0b/wDA7/eTn5/vcDIR6YmE2BhONetpXvMWTLmQged9k5LSbfpDSHrEWGsP/yBjvMBL1tqp3bfrrLUDut82wK49tw/ysXlAHkBWVtYJJSUlfRJcxEkrttVxyR/foKOunPJHf4xta977PmMMgUDAwXQi0lNer5eSkhIGzL6OtJOvxL95MdUv/J6sYZkUFxc7HU8cZoxZYq3NOdj7jnoSue1qYIdsYdbaQmttjrU2JyMj42ifTsRxn5Q3MOfBRbja/VQ+9av9yhNAVlaWQ8lEpLdKS0sBqHv/39S89lcSRx/P0Nzfs31X82E+UqLdkRaoCmPMMIDu/1b2XSSR0FVc3cS1/1pIXIyLO06IIz6w/4usx+OhoKDAoXQi0lv7/sHTuOI1Kp/8Be60TEZc/ydWbKtzLpiEvCMtUC8Ac7rfngM83zdxRELXjrpmcu9fSEdnAN9NM/nODVdTWFhIdnY2xhiys7MpLCwkNzfX6agi0kMFBQV4PJ69t1uKl1P35J0MTE3iqsKPeXXVTgfTSSg77BwoY8xjwJnAYKAC+AXwHPAEkAWUAFdaaw+7S2NOTo4tKio6usQiDqhubOXKez+mqqGVR2+epQ2CRSKIz+cjPz+f0tJSsrKyKCgo4PzLriDv30UsLa3je2eP53tnj8flMk5HlSD7vDlQPZpE3ldUoCQc7Wpq45r7F7K1upFHbpzJid50pyOJSBC0tHeS/+xqnl5axlmTMrn7yuNI88Q6HUuCqF8nkYtEsurGVq7+5wI2VzVS+PUclSeRKJIQG8P/fXUav758KvM3VnHp3z5g3c7dTseSEKECJXIIFbtbuOq+jymuaeKBOScye4KuIhWJNsYYvj4rm7l5s2hu6+RLf/+Q55dvdzqWhAAVKJGD2FHXzFX3fUx5fQsP33ASp40f7HQkEXHQCdnpvPTd05g2YgDfm7ucO59bRUt7p9OxxEEqUCKfsa3Wz5X3fUxNYxv/vnEmM8cMcjqSiISAzJQEfDfP5ObTR/OfBaVc/rcP2VjR4HQscYgKlMg+NlU2cuV9H9PQ0oHv5pmckK19skXkU7ExLvIvnsyDN5xIVUMrl/z1A+YuKiWYF2RJaFCBEulWVFzLV/7xEe2dAR67eRbTRg5wOpKIhKgvTMzk1e+dTk52Oj95ZhXffmwZ9c3tTseSIFKBEgFeW11O7v0LGZQUx7PfPJXJw1OdjiQiIS4zNYF/f+MkfnTBRF5bXc75d7/P+xuqnI4lQaICJVHv4Y+Kuc23hMnDU3nqtlMYle45/AeJiAAul+GbZ47j2W+eQnKCm+seWMTPnl1FY2uH09Gkn6lASdQKBCy/e3U9v3hhDWdPGsKjN80iPSnO6VgiEoamjRzAS985jbzZY3hsUSkX/vl9FmypcTqW9CMVKIlK/rYOvv3YUu59bzO5M7O499rjSYyLcTqWiISxhNgYfnbRMTx5y8m4jOHqfy7gF8+vpqFFc6MikQqURJ1ttX6+/PePeG11OT+7aBK/uXwq7hj9KIhI38jxpvPq905nzsle/r2ghHP/+D6vryl3Opb0Mf3WkKiyYEsNl/3tQ7bXNfPA9SeSN3ssxmiDUBHpW544N7+8dArP3HYKAzyx3PLIEm55pIjy+hano0kfUYGSqPGfBSVce/9CBnhief5bp3LmxEynI4lIhJuRNZAXv3MaP75gEu9+UsU5f3yPBz/cSkdnwOlocpRUoCTiNbd18qOnVnDnc6s5ffxgnvvWqYzJSHY6lohEidgYF7edOZZ5P5jNjKwB/OrFtVx8zwd8tKna6WhyFFSgJKJtqmzgsr99wJNLyvjOWeO4f86JpCbEOh1LRKJQ9qAk/v2Nk7jv6yfQ1NbBNfcv5Ju+JZTt8jsdTY6ACpRErGeWlnHh3e+xoWQnFY//nL/knc/cxx51OpaIRDFjDOdPGcqbt5/B7edO4O31lZzzx/f44xsbaNLaUWHFBHP/npycHFtUVBS055Po1NzWyS9eWM0TRWW0bV9L5XO/o7OxFgCPx0NhYSG5ubkOpxQRge11zfz2lXW8vHIng5Pj+f4547nqxFHE6srgkGCMWWKtzTno+1SgJJKs2FbHD55YztbqJuzqVyl5+R9g95+smZ2dTXFxsTMBRUQOYmnpLv77lXUsLt7FmIwkfnzBJM6bPERXCTtMBUoiXntngL++vYm/vrOJzJR4/veK6cyemHnQHdKNMQQCugJGREKLtZY311Xyu1fXsbmqieOzBnDHeRM5ZewgFSmHfF6Bcgc7jEhf21TZyO1PLGdlWT1fmjGCX146hbTEWLKysigpKTng8VlZWQ6kFBH5fMYYzp08hC9MzOCJojLueWsjufcvZObodO44byInjU53OqLsQydZJWx1Biz3z9/CxffMZ1utn7/nHs/dVx1HWmLXVXYFBQV4PPtvDOzxeCgoKHAirohIj7hjXFwzM4t3/+tMfnnJZLZUN3HlfR9z7f0LWVJSu99jfT4fXq8Xl8uF1+vF5/M5lDr66BSehKU1O+r56TOrWFlWz1mTMvndl48lMzXhgMf5fD7y8/MpLS0lKyuLgoICTSAXkbDS3NaJb2EJ/3h3MzVNbcwcnc5tZ46lbPEb3HJLHn7/p8sg6EKZvqU5UBIxmts6+dNbG7h//lYGemL5xSVT+OK0YZofICIRr6m1g7mLt3H//C3srG+BXWVUzX8U//oP9rtYRhfK9B0VKIkI8zdWcedzqymp8XNlzkh+dtExDPDEOR1LRCSo2joCPL98O9+/72ViB42iva6c3Qufpmn1W9iONl0o04dUoCSs/eXBx/jDm1tg1HHQUMmN0zz8v1uucjqWiIijvN7RVMYNJW3WV4kfPpHOpl00LH2ZATVrKF6/0ul4EUFX4UlYam7r5Jt/fY63dyRgh0yi/r2H2b34OX4XH8uY5A6d4xeRqFZQ8Bvy8vIof2QB8aOOJW3WVxhw+rXEGMv35y7julO8zBg1QFMc+olGoCTkWGt5edVO/vuV9Wyva6Zp7bvsevdBOhtq9j5G5/hFRA68UOYHP/9v6jKm8VRRGQ2tHRw7Io05p3j54rRhJMTGOB037OgUnoSNjzfX8LvX1rNiWx3HDEvl3T/cRsu2NQc8Tuf4RUQOrbG1g2eXbeffHxWzsbKRgZ5YvnZSFl87cRTZg5Kcjhc2VKAk5K0v383/vLqedz6pYlhaArefO4EvHz+SsWNGH3QxTI1AiYgcnrWWjzfX8PDHxbyxtoKAhVlj0rkyZxQXTh1GYpxGpT6PCpSErJKaJu55axPPLCsjJd7Nt74wjjmnePcONft8PvLytM6JiMjRKq9v4emlZTxRtI2SGj8p8W4uOW44V+aMYvrINM2VOggVKAk5JTVN/PXtTTyzbDtul2HOKV6+eebYgy5LoMUwRUT6jrWWRVtrebxoG6+s2klLe4CJQ1K44oSRXDJ9OEPTPl2UONpff1WgxHF7fgh37G5j+DnfwDVmFrHuGHJnZnPrGWMOuoq4iIj0r4aWdl5auZPHF29j+bY6jIGZo9O5dPoImtZ/wPe/eXNUnwFQgRJH+Xw+bvtZAfHTL8Yz6XRsZwctq9/kV1efxm3XX+N0PBERAbZWN/HC8h08v2I7W6qaINCJf8sSmta+R/Omhdj2FiC65qCqQElQfHao9ze/KWDMyedzzV0PwrDJBFr9NCx/lYai5+lsrI2qH0IRkXBhrWXNjt3MvvYHeI6ZjTs1g0BbC81bivBv+JiWLUV0tjQ6HTMoVKCk3+032TvGTdLE00g76XJih4yjs3EXu5e8QMOyV7CtTXs/RksRiIiELq/XS0lJKfEjJ5N0zGwSJ5yMOzkdOjs445hhXDB1KOccM4SMlPi9HxNpc6ZUoKTPHOqHw+v1UlbbSPJxF5Iy/XxikgbSXlNG7Jb5sHUBJVs3H/C5NAIlIhK6DrwK2pA6ZjqX3HYnpYF0Smv9GAM52QM5f8pQmjYu5GffuSmi5kypQEmfONSSAj+9+1/8+ZUVeCacDMbQvGkxDUtfpKV4BcbAI488oqUIRETC0KH+aLbWsr68gdfXlPPa6nLWlzcA0F67g+YtRTRvXkzLttXQ2R7WfyyrQEmf6BrO7VrU0sQnkXTMbFJmXERc5mhobaJ++Ws0LnuFjvqKvR+z5wcn0oZ1RUTkU9tq/Uw550oSxuaQkDUNV2w8gbYWWkqW07JlCeveeZoRAxL3Pj5cfieoQEmvHOob2xXjJj5rGsnHnoNnwskYdxxtFZtpWPoy9/zwOr59q0aZRESi1Z4/so07jvisaXjG5pA49kTcaUMAGD04iVPHDSKwYx333PltmnZV7f3YUP19oQIlPXaw03TJw0bzlTv+l3dLmsEzkM7mBprWvkfTqjdoq9isUSYRETnkNI/f3HM/qRNn8eGmahZsqcHf1okNdNJWvomW4uU0lyyndft6skcOD7lTfSpQ0mN7/oJwJabimXAKSVPOJGHUVAgEGJ/ayZIn/8quNe9DZzsQun81iIhI8B3uD+m2jgBpo6cRnz2dBO904odPwrhisB3ttO7cwB3XXcpJowdxQvZAkuPdDh5JFxUo6ZF6fzsjT+5a7DIhezrGFUN7TRmNq9/Cv+Zt2ndXa5RJRESOyn7zaeMSSRg1lfhRU0kbezwxmWPoDFhcBqYMT+Ok0emcNDqdE73ppCcduNVXf1OBkgPsKULbyqsYNfMiJp6XyxZ/HO2dlva6cvzr5tO0/n3aK7cCWnJARET6xudtEn/5FVexrLSORVtrWLi1luXb6mjt6FovcHxmMgM6alj8yuPsWPUhw1Pc/f5HvApUFDvYiNHudsOdf38cd/bxXSNN7lg6G6qZPTqZaWnt3PX9GzUZXERE+k1Pz2a0dnSyqqyehVtree7D1XxS04YrPgmA8sd+Skz15n79/aQCFeEO9Y24b8uPzfDiGTeTpIknEztkHADtu3bSvGkh/k8+7JrAl52lyeAiIhKS9qyMHjt4FPHDJ9G07j1se2u/niFRgYpghxoK/cs/Crnr3sfYnTQSz9gTcQ8YCkDrjvX4Ny6kedNC2qtL9/tc2lpFRERClcvl4mCdpT9/d31egXJ+irscoDcjQPn5+d3lyRCbOZrE0ceTOPp47lqVDLNvI7mthZbSldQveBL/poUEmuoO+bxZWVn9c0AiIiJHKSsra+/k88/e7wQVqBDz2RGlkpIS8vLyAA4oUdtq/dSkjmfQxV8hcfQMYpIGAtBWuZXdRc+T6t/OtiVvQ2fHfh83aNAgmpubDxi1Kigo6M9DExEROWIFBQUHPePi1O8uncILMfte3rmv7Oxs5i9dy4ItNSzYUsuCLTVsr2sGoLOprmsxsq1LaSleRmfTLrKzsw/5zVZYWAigeU4iIhJWgj1HV3Ogwsjec7zGRezgLOKHTyJ+xDEkZE3duxx+elIcs8akM2vMIKrXfsz/++7N+P1Nez/HvlfNaUK4iIjIkdEcqCDpTVn57GPvvOu3jJ15DlkX3UZz8jDih03EFe8BukaYYmq38quvn8WsMYMYn5mMy2W6PtHJXoYkBA75vLm5uSpMIiIifUwjUH3k8xYG+2yBefgRH9+5swA7YBTxIyYRN3wScYO7JsEZLG2VW2kpW0vr9vW07lhPXNturcMkIiISZDqFt4/+OqV1yLlLo8fy4vtLWLW9nlXb61mzo56VJdUQEwtAZ/PurqK0fR0DO3exZv6rPPfU4zrtJiIi4rCoKFA9KUa9GSXqLZfLBbGJxGVkE5vhJW7oOOKGjCUuIxvTXZZSE9xMHZHGa48W0lq+ibaKzXTs2rH3c2gdJhERkdAR8QWqp8Xo865w680qpg8/4uMX//c3qjvjGTxmGhNmfoEN5Q2QlL73MZ3Nu2kr30RicxV/L/gpU4enMSo9EWNMn+UQERGR/hPxBaqnhaS3q5g+9IiPX/7f36huczN4zGSOPe18KlpiqPBbTEzX/Hvb2U7nrh2MSnOzcfF7+HdspK2qmM7dVYcc3erPkTARERHpGxF/FV5paWmP7j/kKqYTprJgSw2bqxrZXNnE5qpGVm4tp7YtFXNhPoO7H7e8tBJbt52m7Rtoqy6hvaqE9trtEOjAnZ3NnwsKuk4jNlTvXYfpYIVoz32a5yQiIhKeomIEqjNg2VnfzANPvMif/+Uj4EnHPWAYsQOH4R4wbO9yAQCJsTGMyUhizcdvU7ftE9prymivKaNj1w5sR+shM2j+koiISGSJ+FN4e06JtbqTiMsYjXvgUBIHj2L6qefQkTCAsl3NtHXuU246O2ivLye2pZ5Zx47j/FNmMDYjmbGZyQxLTcDlMoc83Xcomr8kIiISWSL+FN7eU2JPL4cJZwIQ77LEpqYxbpCH86YMJXuQh+x0D1mDPAxLSyRmz0KUh3Co033aR05EREQiYgRqj81VjdQ3t5Od7iE9KQ5jPr8kfZ7Pm+gNmr8kIiIS6SJ+BGqPsRnJffa5DjfRW4VJREQkekXUCJSIiIhIX/m8EShXsMOIiIiIhDsVKBEREZFeUoESERER6SUVKBEREZFeUoESERER6SUVKBEREZFeUoESERER6SUVKBEREZFeUoESERER6aWjKlDGmAuMMZ8YYzYZY37SV6FEREREQtkRFyhjTAzwN+BCYDJwtTFmcl8FExEREQlVRzMCdRKwyVq7xVrbBswFLuubWCIiIiKh62gK1Ahg2z63y7rv248xJs8YU2SMKaqqqjqKpxMREREJDf0+idxaW2itzbHW5mRkZPT304mIiIj0O/dRfOx2YNQ+t0d233dIS5YsqTbGlBzFc4a6wUC10yEcFO3HD/oa6Pij+/hBXwMdf2Qdf/ah3mGstUf0GY0xbmADcDZdxWkxcI21ds0RfcIIYIwpstbmOJ3DKdF+/KCvgY4/uo8f9DXQ8UfP8R/xCJS1tsMY823gdSAGeCCay5OIiIhEj6M5hYe19hXglT7KIiIiIhIWtBJ53yp0OoDDov34QV8DHb9E+9dAxx8ljngOlIiIiEi00giUiIiISC+pQImIiIj0kgrUETjcJsrGmNuNMWuNMSuNMW8ZYw65jkQ46sHx32qMWWWMWW6M+SDS9kjs6SbaxpivGGOsMSbiLuntwffA9caYqu7vgeXGmJucyNlfevI9YIy5svt1YI0x5tFgZ+xPPfj/f/c+/+83GGPqHIjZr3rwNcgyxrxjjFnW/bvgIidy9pceHH929++/lcaYd40xI53I2a+stfrXi390LdmwGRgDxAErgMmfecwXAE/327cBjzudO8jHn7rP25cCrzmdO5jH3/24FOB9YAGQ43RuB74Hrgf+6nRWB49/PLAMGNh9O9Pp3ME8/s88/jt0LXPjePYgfw8UArd1vz0ZKHY6d5CP/0lgTvfbZwGPOJ27r/9pBKr3DruJsrX2HWutv/vmArpWaY8UPTn+3fvcTAIi6UqFnm6i/Wvgf4CWYIYLkmjfSLwnx38z8Ddr7S4Aa21lkDP2p97+/78aeCwoyYKnJ18DC6R2v50G7Ahivv7Wk+OfDLzd/fY7B3l/2FOB6r0ebaK8jxuBV/s1UXD1dBPpbxljNgO/B74bpGzBcNjjN8YcD4yy1r4czGBB1NOfga90D98/ZYwZdZD3h6ueHP8EYIIx5kNjzAJjzAVBS9f/evwa2D19YTSf/iKNFD35GvwSuNYYU0bXeonfCU60oOjJ8a8Avtz99peAFGPMoCBkCxoVqH5kjLkWyAH+1+kswWat/Zu1dizwY+BOp/MEizHGBfwRuMPpLA57EfBaa6cBbwAPO5wn2Nx0ncY7k64RmH8aYwY4GcghXwOestZ2Oh3EAVcDD1lrRwIXAY90vz5Eix8CZxhjlgFn0LXlW0R9H0TT/8y+0qNNlI0x5wD5wKXW2tYgZQuG3m4iPRe4vD8DBdnhjj8FmAq8a4wpBmYBL0TYRPLDfg9Ya2v2+b6/HzghSNmCoSc/A2XAC9badmvtVrr2DR0fpHz9rTevAV8j8k7fQc++BjcCTwBYaz8GEujaaDcS9OQ1YIe19svW2hl0/S7EWlsXtIRBoALVe4uB8caY0caYOLpeIF7Y9wHGmBnAfXSVp0ia+wA9O/59f1FcDGwMYr7+9rnHb62tt9YOttZ6rbVeuubAXWqtLXImbr/oyffAsH1uXgqsC2K+/nbY4weeo2v0CWPMYLpO6W0JYsb+1JPjxxgzCRgIfBzkfMHQk69BKXA2gDHmGLoKVFVQU/afnrwGDN5nxO2nwANBztjvVKB6yVrbAezZRHkd8IS1do0x5i5jzKXdD/tfIBl4svsy3gNeXMJVD4//292Xbi8HbgfmOJO27/Xw+CNaD78G3+3+HlhB1xy4651J2/d6ePyvAzXGmLV0TaD9L2ttjTOJ+1Yvfga+Bsy13ZdhRZIefg3uAG7u/hl4DLg+Ur4WPTz+M4FPjDEbgCFAgSNh+5G2chERERHpJY1AiYiIiPSSCpSIiIhIL6lAiYiIiPSSCpSIiIhIL6lAiYiIiPSSCpSIiIhIL6lAiYiIiPTS/wfEi1KAsLsuYQAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 720x504 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
"source": [
- "sig_p = lambda x: sqrt(1 - 4*m_p**2/x)\n",
- "g_s = lambda s, m_q, g_q: g_q*s/m_q**2 * (sig_p(s)/sig_p(m_q**2))**3 * Heaviside(s, 4*m_p**2)\n",
"\n",
- "def model(s, m_q, g_q, m_w, g_w, e_w, a, b, c):\n",
- " part1 = (m_q)**4/((m_q**2 - s)**2 + m_q**2*g_s(s, m_q, g_q)**2)\n",
- " part2 = 1 + (e_w * 2*s * (m_w**2 - s))/((m_w**2 - s)**2 + m_w**2*g_w**2)\n",
- " part3 = c*(1 + a*s + b*s**2)**2\n",
- " return part1 * part2 * part3\n",
+ "data = np.loadtxt('../data/SND-VFF.txt')\n",
+ "cov_stat = np.loadtxt('../data/BABAR-StatCov.txt')\n",
+ "x_data = data[:,0]\n",
+ "y_data = data[:, 1]\n",
"\n",
- "def foo(a ,b ,c):\n",
- " return a + b - c"
+ "p0 = [0.9, 0.2, 0.81, 0.04, 0.02, -1, 0.84, 1.55] # in GeV\n",
+ "\n",
+ "p, pcov = curve_fit(model, x_data, y_data, p0)#, sigma=cov_stat)\n",
+ "x_model = np.linspace(x_data[0], x_data[-1], 200)\n",
+ "\n",
+ "plt.figure(figsize=[10, 7])\n",
+ "plt.scatter(x_data, y_data, c='black')\n",
+ "plt.plot(x_model, model(x_model, *p))"
]
},
{
"cell_type": "code",
- "execution_count": 100,
- "id": "elect-storm",
+ "execution_count": 64,
+ "id": "numerous-concord",
"metadata": {},
- "outputs": [],
- "source": [
- "s_ = symbols(\"s\")\n",
- "y = model(s_, *var)"
- ]
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[<matplotlib.lines.Line2D at 0x7f25c96d1520>]"
+ ]
+ },
+ "execution_count": 64,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOZUlEQVR4nO3cf6zd9V3H8edrVNCFjV8Fxij1onTRsiVuOYEt/kKBUkigiyMGzLLOoE0WMXGosYZEGOwPUDfMMvxRBxmSOMAlumvm0jB+hGQB5ABzUpS1AxwFNgpFDCEDu73943wnl7tb7rk9p+f09vN8JDc93+/303ve+eS2z57vubepKiRJ7XrLtAeQJE2XIZCkxhkCSWqcIZCkxhkCSWrcimkPsC9WrlxZMzMz0x5DkpaNlStXsnXr1q1VtX7+tWUZgpmZGfr9/rTHkKRlJcnKhc57a0iSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGjeWECRZn+SxJDuSbF7g+mFJbu2u359kZt711UleTvIH45hHkjS8kUOQ5BDgeuBcYC1wcZK185ZdArxYVacA1wHXzrv+aeAro84iSVq6cbwiOA3YUVWPV9VrwC3AhnlrNgA3dY+/CJyZJABJPgg8AWwbwyySpCUaRwhOBJ6ac7yzO7fgmqraA7wEHJPkcOCPgE8s9iRJNiXpJ+nv2rVrDGNLkmD6bxZfCVxXVS8vtrCqtlRVr6p6xx577P6fTJIasWIMn+Np4KQ5x6u6cwut2ZlkBXAE8AJwOnBhkj8FjgR+kOR7VfXZMcwlSRrCOELwALAmyckM/sK/CPiNeWtmgY3AvcCFwJ1VVcAv/nBBkiuBl42AJE3WyCGoqj1JLgW2AocAN1bVtiRXAf2qmgVuAG5OsgPYzSAWkqQDQAb/MF9eer1e9fv9aY8hSctKkgerqjf//LTfLJYkTZkhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGjSUESdYneSzJjiSbF7h+WJJbu+v3J5npzp+d5MEk/979+qvjmEeSNLyRQ5DkEOB64FxgLXBxkrXzll0CvFhVpwDXAdd2558Hzq+q9wAbgZtHnUeStDTjeEVwGrCjqh6vqteAW4AN89ZsAG7qHn8RODNJqurhqnqmO78N+Ikkh41hJknSkMYRghOBp+Yc7+zOLbimqvYALwHHzFvzIeChqnp1DDNJkoa0YtoDACQ5lcHtonVvsmYTsAlg9erVE5pMkg5+43hF8DRw0pzjVd25BdckWQEcAbzQHa8C/hH4SFV9a29PUlVbqqpXVb1jjz12DGNLkmA8IXgAWJPk5CSHAhcBs/PWzDJ4MxjgQuDOqqokRwJfBjZX1dfGMIskaYlGDkF3z/9SYCvwH8BtVbUtyVVJLuiW3QAck2QHcBnww28xvRQ4BfiTJF/vPo4bdSZJ0vBSVdOeYcl6vV71+/1pjyFJy0qSB6uqN/+8P1ksSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0zBJLUOEMgSY0bSwiSrE/yWJIdSTYvcP2wJLd21+9PMjPn2h935x9Lcs445pEkDW/kECQ5BLgeOBdYC1ycZO28ZZcAL1bVKcB1wLXd710LXAScCqwH/rL7fJKkCVkxhs9xGrCjqh4HSHILsAF4dM6aDcCV3eMvAp9Nku78LVX1KvBEkh3d57t3DHP9iE/88zYefeZ/9senlqT9bu07384V55869s87jltDJwJPzTne2Z1bcE1V7QFeAo4Z8vcCkGRTkn6S/q5du8YwtiQJxvOKYCKqaguwBaDX69W+fI79UVJJWu7G8YrgaeCkOcerunMLrkmyAjgCeGHI3ytJ2o/GEYIHgDVJTk5yKIM3f2fnrZkFNnaPLwTurKrqzl/UfVfRycAa4F/HMJMkaUgj3xqqqj1JLgW2AocAN1bVtiRXAf2qmgVuAG7u3gzezSAWdOtuY/DG8h7gd6rq+6POJEkaXgb/MF9eer1e9fv9aY8hSctKkgerqjf/vD9ZLEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1DhDIEmNMwSS1LiRQpDk6CS3J9ne/XrUXtZt7NZsT7KxO/fWJF9O8p9JtiW5ZpRZJEn7ZtRXBJuBO6pqDXBHd/wGSY4GrgBOB04DrpgTjD+vqp8B3gv8fJJzR5xHkrREo4ZgA3BT9/gm4IMLrDkHuL2qdlfVi8DtwPqqeqWq7gKoqteAh4BVI84jSVqiUUNwfFU92z3+DnD8AmtOBJ6ac7yzO/f/khwJnM/gVYUkaYJWLLYgyVeBdyxw6fK5B1VVSWqpAyRZAXwB+ExVPf4m6zYBmwBWr1691KeRJO3FoiGoqrP2di3Jd5OcUFXPJjkBeG6BZU8DZ8w5XgXcPed4C7C9qv5ikTm2dGvp9XpLDo4kaWGj3hqaBTZ2jzcCX1pgzVZgXZKjujeJ13XnSPJJ4Ajg90acQ5K0j0YNwTXA2Um2A2d1xyTpJfkcQFXtBq4GHug+rqqq3UlWMbi9tBZ4KMnXk/zWiPNIkpYoVcvvLkuv16t+vz/tMSRpWUnyYFX15p/3J4slqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXGGQJIaZwgkqXEjhSDJ0UluT7K9+/Wovazb2K3ZnmTjAtdnkzwyyiySpH0z6iuCzcAdVbUGuKM7foMkRwNXAKcDpwFXzA1Gkl8DXh5xDknSPho1BBuAm7rHNwEfXGDNOcDtVbW7ql4EbgfWAyQ5HLgM+OSIc0iS9tGoITi+qp7tHn8HOH6BNScCT8053tmdA7ga+BTwymJPlGRTkn6S/q5du0YYWZI014rFFiT5KvCOBS5dPvegqipJDfvESX4O+Omq+niSmcXWV9UWYAtAr9cb+nkkSW9u0RBU1Vl7u5bku0lOqKpnk5wAPLfAsqeBM+YcrwLuBj4A9JI82c1xXJK7q+oMJEkTM+qtoVngh98FtBH40gJrtgLrkhzVvUm8DthaVX9VVe+sqhngF4BvGgFJmrxRQ3ANcHaS7cBZ3TFJekk+B1BVuxm8F/BA93FVd06SdABI1fK73d7r9arf7097DElaVpI8WFW9+ef9yWJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGpaqmPcOSJdkF/Ne05xiTlcDz0x7iAOFevM69eCP343X7uhfPA1TV+vkXlmUIDiZJ+lXVm/YcBwL34nXuxRu5H6/bH3vhrSFJapwhkKTGGYLp2zLtAQ4g7sXr3Is3cj9eN/a98D0CSWqcrwgkqXGGQJIaZwgmJMn6JI8l2ZFk8wLXL0vyaJJvJLkjyU9OY85JWGwv5qz7UJJKctB+2+Awe5Hk17uvjW1J/n7SM07KEH9GVie5K8nD3Z+T86Yx5yQkuTHJc0ke2cv1JPlMt1ffSPK+kZ6wqvzYzx/AIcC3gJ8CDgX+DVg7b82vAG/tHn8MuHXac09rL7p1bwPuAe4DetOee4pfF2uAh4GjuuPjpj33FPdiC/Cx7vFa4Mlpz70f9+OXgPcBj+zl+nnAV4AA7wfuH+X5fEUwGacBO6rq8ap6DbgF2DB3QVXdVVWvdIf3AasmPOOkLLoXnauBa4HvTXK4CRtmL34buL6qXgSoqucmPOOkDLMXBby9e3wE8MwE55uoqroH2P0mSzYAf1cD9wFHJjlhX5/PEEzGicBTc453duf25hIGtT8YLboX3cvck6rqy5McbAqG+bp4F/CuJF9Lcl+SH/nvAQ4Sw+zFlcCHk+wE/gX43cmMdkBa6t8pb2rFyONorJJ8GOgBvzztWaYhyVuATwMfnfIoB4oVDG4PncHgVeI9Sd5TVf89zaGm5GLg81X1qSQfAG5O8u6q+sG0B1vufEUwGU8DJ805XtWde4MkZwGXAxdU1asTmm3SFtuLtwHvBu5O8iSD+5+zB+kbxsN8XewEZqvqf6vqCeCbDMJwsBlmLy4BbgOoqnuBH2fwH7C1aKi/U4ZlCCbjAWBNkpOTHApcBMzOXZDkvcDfMIjAwXofGBbZi6p6qapWVtVMVc0weL/kgqrqT2fc/WrRrwvgnxi8GiDJSga3ih6f4IyTMsxefBs4EyDJzzIIwa6JTnngmAU+0n330PuBl6rq2X39ZN4amoCq2pPkUmArg++OuLGqtiW5CuhX1SzwZ8DhwD8kAfh2VV0wtaH3kyH3oglD7sVWYF2SR4HvA39YVS9Mb+r9Y8i9+H3gb5N8nMEbxx+t7ltoDjZJvsDgHwAru/dErgB+DKCq/prBeyTnATuAV4DfHOn5DtJ9lCQNyVtDktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktS4/wOewbL27S45zAAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": []
},
{
"cell_type": "code",
- "execution_count": 114,
- "id": "capital-lounge",
+ "execution_count": 81,
+ "id": "chief-kentucky",
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[0.032, 0. , 0. , ..., 0. , 0. , 0. ],\n",
+ " [0. , 0.032, 0. , ..., 0. , 0. , 0. ],\n",
+ " [0. , 0. , 0.032, ..., 0. , 0. , 0. ],\n",
+ " ...,\n",
+ " [0. , 0. , 0. , ..., 0.013, 0. , 0. ],\n",
+ " [0. , 0. , 0. , ..., 0. , 0.013, 0. ],\n",
+ " [0. , 0. , 0. , ..., 0. , 0. , 0.013]])"
+ ]
+ },
+ "execution_count": 81,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
- "A = Matrix([y]).jacobian(var)\n",
- "A = A.subs({var[i]:p0[i] for i in range(len(p0))})\n",
- "\n",
- "sol = []\n",
- "\n",
- "s = np.linspace(0, 1, 100)\n",
- "for i in range(len(s)):\n",
- " sol.append(A.subs({s_: s[i]}))\n",
- " \n"
+ "np.eye(len(x_data)) * np.array([0.032 if i<=2 else 0.013 for i in range(len(x_data))])"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "magnetic-possibility",
+ "id": "administrative-india",
"metadata": {},
"outputs": [],
"source": []
@@ -86,7 +161,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "excessive-exhibition",
+ "id": "neither-store",
"metadata": {},
"outputs": [],
"source": []
@@ -94,7 +169,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "brave-addiction",
+ "id": "elder-celebration",
"metadata": {},
"outputs": [],
"source": []
@@ -102,7 +177,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "continental-snapshot",
+ "id": "serious-forum",
"metadata": {},
"outputs": [],
"source": []
diff --git a/sesh1/prog/data/CMD2-VFF.txt b/sesh1/prog/data/CMD2-VFF.txt
@@ -1,5 +1,5 @@
-# 43+10+29 = 82 data points
-# normalization uncertainty, totally correlated: 0.6% for the first 43 points, 0.7% for the next 10 points, 0.8% for the last 29 points
+# 43+10+29 = 82 data points
+# normalization uncertainty, totally correlated: 0.6% for the first 43 points, 0.7% for the next 10 points, 0.8% for the last 29 points
# s/GeV |F(s)|^2 sigma_stat(|F(s)|^2)
0.37271025 8.10930218 1.14188265
0.38502025 9.87674407 0.732981756
@@ -82,4 +82,4 @@
0.8836 5.00944749 0.115959433
0.9025 4.59104999 0.0991893516
0.917764 4.12106848 0.0862549217
-0.9409 3.81892364 0.0881290071
-\ No newline at end of file
+0.9409 3.81892364 0.0881290071
diff --git a/sesh1/prog/fit_single.png b/sesh1/prog/fit_single.png
Binary files differ.
diff --git a/sesh1/prog/fit_single.py b/sesh1/prog/fit_single.py
@@ -1,46 +0,0 @@
-#!/usr/bin/env python3.9
-
-import numpy as np
-from scipy.optimize import curve_fit
-import matplotlib.pyplot as plt
-
-global m_p; m_p = 0.13957
-
-sig_p = lambda x: np.sqrt(1 - 4*m_p**2/x)
-g_s = lambda s, m_q, g_q: g_q*s/m_q**2 * (sig_p(s)/sig_p(m_q**2))**3 * np.heaviside(s, 4*m_p**2)
-
-def model(s, m_q, g_q, m_w, g_w, e_w, a, b, c):
- part1 = (m_q)**4/((m_q**2 - s)**2 + m_q**2*g_s(s, m_q, g_q)**2)
- part2 = 1 + (e_w * 2*s * (m_w**2 - s))/((m_w**2 - s)**2 + m_w**2*g_w**2)
- part3 = c*(1 + a*s + b*s**2)**2
- return part1 * part2 * part3
-
-def main():
- data = np.loadtxt('./data/SND-VFF.txt')
- s = data[:,0]
- F2 = data[:, 1]
- sigma = data[:, 2]
-
- p0 = [0.7, 0.2, 1, 0.2, 2e-3, 14.4, 25, 12.2] # in GeV
- popt, pcov = curve_fit(model, s, F2, p0, sigma=sigma)
- popt, uncert = np.round(popt, 3), np.round(np.sqrt(np.diagonal(pcov)), 3)
-
- s_model = np.linspace(s[0], s[-1], 500)
-
- plt.figure(figsize=[10, 7])
- plt.title('SND DATA FIT')
- plt.scatter(s, F2, marker='.', c='black')
- plt.plot(s_model, model(s_model, *popt), color='red')
-
- plt.annotate(r'$M_{\rho} = $' + f'({popt[0]}' + r'$\pm$' + f'{uncert[0]}) GeV', (0.2, 40))
- plt.annotate(r'$\Gamma_{\rho} = $' + f'({popt[1]}' + r'$\pm$' + f'{uncert[1]}) GeV', (0.2, 38))
- plt.annotate(r'$M_{\omega} = $' + f'({popt[2]}' + r'$\pm$' + f'{uncert[2]}) GeV', (0.2, 36))
- plt.annotate(r'$\Gamma_{\omega} = $' + f'({popt[3]}' + r'$\pm$' + f'{uncert[3]}) GeV', (0.2, 34))
-
- plt.savefig('fit_single.png')
-
- print(*popt, sep='\n')
-
-
-if __name__ == "__main__":
- main()
diff --git a/sesh1/prog/t0-method/t0-imp.py b/sesh1/prog/t0-method/t0-imp.py
@@ -1,14 +0,0 @@
-#!/usr/bin/env python3.9
-
-import numpy as np
-import sympy as sp
-from t0-method import *
-
-def main():
- p0 = [0.765 ,0.115 ,0.813 ,0.041 ,-0.006 ,-1.005 ,1.014 ,1.854] # guess
- data = np.loadtext('../data/
-
-
-
-if __name__ == "__main__":
- main()
diff --git a/sesh1/prog/t0-method/t0-method.py b/sesh1/prog/t0-method/t0-method.py
@@ -1,29 +0,0 @@
-#!/usr/bin/env python3.9
-
-import numpy as np
-import sympy as sp
-
-def design_matrix(model, x_data, p0, var_str):
- var = sp.Matrix(sp.symbols(var_str))
- x = symbols("x")
-
- A = sp.Matrix(model(x, *var)).jacobian(var)
- A = A.subs({var[i]: p0[i] for i in range(len(p0))})
-
- return np.array([np.array([A.subs({x: x_data[i]]})) for i in range(len(s))])
-
-def cov_syst(model, x_data, p0, cov_relsyst)
- return cov_relsyst @ model(x_data, *p0) @ model(x_data, *p0).T
-
-
-def t0_fit(model, var_str, x_data, y_data, p0, cov_stat, cov_relsyst, iterations=100):
-
- A = design_matrix(model, x_data, p0, var_str)
-
- for i in range(iterations):
- P = np.linalg.inv(cov_syst(model, x_data, p0, cov_relsyst) + cov_stat)
- p0 = (A.T @ P @ A).inv() @ A.T @ P @ y
-
- dp = (A.T@P@A).inv()**(1/2)
-
- return p0, dp
diff --git a/sesh1/prog/tp43.py b/sesh1/prog/tp43.py
@@ -1,32 +0,0 @@
-#!/usr/bin/env python3.9
-
-from sympy import *
-
-def main():
- p_, n, nprime = symbols('p_, n, nprime')
- vr = Matrix([p_, n, nprime]) # array of variables
-
- p = Matrix([1, 0.8, 1.2]) # initial guess
- y = Matrix([0.8, 1.2, 1, 1]) # observables
- P = eye(4)*25 # inverse convariace matrix
- f = Matrix(([n*p_, nprime*p_, n, nprime])) # parameters
-
-
- for _ in range(10):
-
- f0 = f.subs({p_:p[0], n:p[1], nprime:p[2]})
- I0 = y - f0
-
- # design matrix
- A = f.jacobian(vr).subs({p_:p[0], n:p[1], nprime:p[2]})
-
- dp=(A.T @ P @ A).inv() @ A.T @ P @ I0
-
- p = p + dp
-
- cov = (A.T @ P @ A).inv()
- for i in range(3):
- print(f'{vr[i]}:\t {p[i]}\t\\pm {sqrt(cov[i, i])}')
-
-if __name__ == "__main__":
- main()
diff --git a/sesh1/prog/tp51.py b/sesh1/prog/tp51.py
@@ -1,41 +0,0 @@
-#!/usr/bin/env python3.9
-
-import numpy as np
-
-def main():
-
- # one experiment
- sig = 0.2; s = 0.2
- p0 = 0.9 # guess
-
- y0 = np.array([0.8, 1.2])
- A = np.array([1, 1])
- cov_sta = np.array([[sig**2, 0], [0, sig**2]])
-
- for i in range(10):
- cov_sys = np.array([[s**2*p0**2, s**2*p0**2], [s**2*p0**2, s**2*p0**2]])
- P = np.linalg.inv(cov_sys+cov_sta)
- p0 = (A.T@P@A)**(-1) * A.T@P@y0
-
- dp = (A.T@P@A)**(-1/2)
-
- print('One Experiment')
- print(f'p = {p0} \pm {round(dp, 2)}\n')
-
- # two experiments
-
- p0 = 0.9 # guess
-
- for i in range(10):
- cov_sys = np.array([[s**2*p0**2, 0], [0, s**2*p0**2]])
- P = np.linalg.inv(cov_sys+cov_sta)
- p0 = (A.T@P@A)**(-1) * A.T@P@y0
-
- dp = (A.T@P@A)**(-1/2)
-
- print('Two Experiments')
- print(f'p = {p0} \pm {round(dp, 2)}\n')
-
-
-if __name__ == "__main__":
- main()
