commit 7577c40224a027fea6622a7a99f3f2fc0f73808b
parent 69fac7f5314d23d2ad5492846f036a6ccb7f85ec
Author: miksa <milutin@popovic.xyz>
Date: Thu, 8 Apr 2021 16:12:35 +0200
t0-method done
Diffstat:
7 files changed, 159 insertions(+), 64 deletions(-)
diff --git a/sesh1/prog/Untitled Folder/Untitled.ipynb b/sesh1/prog/Untitled Folder/Untitled.ipynb
@@ -2,52 +2,39 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 23,
+ "execution_count": 52,
"id": "automatic-recommendation",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
- "from scipy.optimize import curve_fit\n",
- "import matplotlib.pyplot as plt\n",
+ "from sympy import *\n",
"global m_p; m_p = 0.13957\n",
- "from ipywidgets import interact"
+ "from sympy.utilities.lambdify import lambdify\n",
+ "\n"
]
},
{
"cell_type": "code",
- "execution_count": 31,
+ "execution_count": 118,
"id": "demographic-distributor",
"metadata": {},
- "outputs": [
- {
- "data": {
- "application/vnd.jupyter.widget-view+json": {
- "model_id": "6d9961d305ea48e0b650fda599667a73",
- "version_major": 2,
- "version_minor": 0
- },
- "text/plain": [
- "interactive(children=(FloatSlider(value=0.7, description='m_q', max=1.0), FloatSlider(value=0.2, description='…"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "data": {
- "text/plain": [
- "<function __main__.somefunc(m_q=0.7, g_q=0.2, m_w=1, g_w=0.2, e_w=0.002, a=-500, b=200, c=0.2)>"
- ]
- },
- "execution_count": 31,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
+ "outputs": [],
"source": [
- "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)\n",
+ "#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",
+ "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",
+ "metadata": {},
+ "outputs": [],
+ "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",
@@ -55,36 +42,43 @@
" part3 = c*(1 + a*s + b*s**2)**2\n",
" return part1 * part2 * part3\n",
"\n",
- "def somefunc(m_q=0.7, g_q=0.2, m_w=1, g_w=0.2, e_w=2e-3, a=-500, b=200, c=0.2): \n",
- " data = np.loadtxt('../data/SND-VFF.txt')\n",
- " s = data[:,0]\n",
- " F2 = data[:, 1]\n",
- " sigma = data[:, 2]\n",
- "\n",
- " p0 = [m_q, g_q, m_w, g_w, e_w, a, b, c] # in GeV\n",
- " popt, pcov = curve_fit(model, s, F2, p0, sigma=np.sqrt(sigma))\n",
- " popt, uncert = np.round(popt, 3), np.round(np.sqrt(np.diagonal(pcov)), 3)\n",
+ "def foo(a ,b ,c):\n",
+ " return a + b - c"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 100,
+ "id": "elect-storm",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "s_ = symbols(\"s\")\n",
+ "y = model(s_, *var)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 114,
+ "id": "capital-lounge",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "A = Matrix([y]).jacobian(var)\n",
+ "A = A.subs({var[i]:p0[i] for i in range(len(p0))})\n",
"\n",
- " plt.figure(figsize=[10, 7])\n",
- " plt.title('SND DATA FIT')\n",
- " plt.scatter(s, F2, marker='.', c='black')\n",
- " plt.plot(s, model(s, *popt), color='red')\n",
- " plt.annotate('Very Bad, A guessing game with the parameters', (0.15, 40))\n",
- " plt.annotate(r'No $\\omega$ resonance recognized by the fit', (0.15, 35))\n",
+ "sol = []\n",
"\n",
- " plt.annotate(r'$M_{\\rho} = $' + f'({popt[0]}' + r'$\\pm$' + f'{uncert[0]}) GeV', (0.7, 40))\n",
- " plt.annotate(r'$\\Gamma_{\\rho} = $' + f'({popt[1]}' + r'$\\pm$' + f'{uncert[1]}) GeV', (0.7, 36))\n",
- " plt.annotate(r'$M_{\\omega} = $' + f'({popt[2]}' + r'$\\pm$' + f'{uncert[2]}) GeV', (0.7, 34))\n",
- " plt.annotate(r'$\\Gamma_{\\omega} = $' + f'({popt[3]}' + r'$\\pm$' + f'{uncert[3]}) GeV', (0.7, 32))\n",
- " \n",
- "interact(somefunc, m_q=(0, 1, 0.1), g_q=(0, 1, 0.1), m_w=(0, 1, 0.1), g_w=(0, 1, 0.1), e_w=(2e-3, 2e-2, 5e-3),\\\n",
- " a=(-100, 100, 0.1), b=(-100, 100, 0.1), c=(0, 100, 0.1))\n"
+ "s = np.linspace(0, 1, 100)\n",
+ "for i in range(len(s)):\n",
+ " sol.append(A.subs({s_: s[i]}))\n",
+ " \n"
]
},
{
"cell_type": "code",
"execution_count": null,
- "id": "stock-julian",
+ "id": "magnetic-possibility",
"metadata": {},
"outputs": [],
"source": []
@@ -92,7 +86,23 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "elect-storm",
+ "id": "excessive-exhibition",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "brave-addiction",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "continental-snapshot",
"metadata": {},
"outputs": [],
"source": []
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
@@ -19,9 +19,9 @@ def main():
data = np.loadtxt('./data/SND-VFF.txt')
s = data[:,0]
F2 = data[:, 1]
- sigma = np.sqrt(data[:, 2])
+ sigma = data[:, 2]
- p0 = [0.7, 0.2, 1, 0.2, 2e-3, 13.4, 45.3, 8.5] # in GeV
+ 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)
@@ -31,15 +31,16 @@ def main():
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'$\Gamma_{\omega}$ bad fit', (0.15, 38))
- plt.annotate(r'$M_{\rho} = $' + f'({popt[0]}' + r'$\pm$' + f'{uncert[0]}) GeV', (0.7, 40))
- plt.annotate(r'$\Gamma_{\rho} = $' + f'({popt[1]}' + r'$\pm$' + f'{uncert[1]}) GeV', (0.7, 38))
- plt.annotate(r'$M_{\omega} = $' + f'({popt[2]}' + r'$\pm$' + f'{uncert[2]}) GeV', (0.7, 36))
- plt.annotate(r'$\Gamma_{\omega} = $' + f'({popt[3]}' + r'$\pm$' + f'{uncert[3]}) GeV', (0.7, 34))
+ 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
@@ -0,0 +1,14 @@
+#!/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
@@ -0,0 +1,29 @@
+#!/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/tp51.py b/sesh1/prog/tp51.py
@@ -0,0 +1,41 @@
+#!/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()
diff --git a/sesh1/iterative.md b/sesh1/t0-method.md
