tprak

Theoretical Physics Practical Training
git clone git://popovic.xyz/tprak.git
Log | Files | Refs

fit_single.py (1571B)


      1 #!/usr/bin/env python3.9
      2 
      3 import numpy as np
      4 from scipy.optimize import curve_fit
      5 import matplotlib.pyplot as plt
      6 
      7 global m_p; m_p = 0.13957
      8 
      9 sig_p = lambda x: np.sqrt(1 - 4*m_p**2/x)
     10 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)
     11 
     12 def model(s, m_q, g_q, m_w, g_w, e_w, a, b, c):
     13     part1 = (m_q)**4/((m_q**2 - s)**2 + m_q**2*g_s(s, m_q, g_q)**2)
     14     part2 = 1 + (e_w * 2*s * (m_w**2 - s))/((m_w**2 - s)**2 + m_w**2*g_w**2)
     15     part3 = c*(1 + a*s + b*s**2)**2
     16     return part1 * part2 * part3
     17 
     18 def main():
     19     data = np.loadtxt('./data/SND-VFF.txt')
     20     s = data[:,0]
     21     F2 = data[:, 1]
     22     sigma = data[:, 2]
     23 
     24     p0 = [0.7, 0.2, 1, 0.2, 2e-3, 14.4, 25, 12.2]   # in GeV
     25     popt, pcov = curve_fit(model, s, F2, p0, sigma=sigma)
     26     popt, uncert = np.round(popt, 3), np.round(np.sqrt(np.diagonal(pcov)), 3)
     27 
     28     s_model = np.linspace(s[0], s[-1], 500)
     29 
     30     plt.figure(figsize=[10, 7])
     31     plt.title('SND DATA FIT')
     32     plt.scatter(s, F2, marker='.', c='black')
     33     plt.plot(s_model, model(s_model, *popt), color='red')
     34 
     35     plt.annotate(r'$M_{\rho} = $' + f'({popt[0]}' + r'$\pm$' + f'{uncert[0]}) GeV', (0.2, 40))
     36     plt.annotate(r'$\Gamma_{\rho} = $' + f'({popt[1]}' + r'$\pm$' + f'{uncert[1]}) GeV', (0.2, 38))
     37     plt.annotate(r'$M_{\omega} = $' + f'({popt[2]}' + r'$\pm$' + f'{uncert[2]}) GeV', (0.2, 36))
     38     plt.annotate(r'$\Gamma_{\omega} = $' + f'({popt[3]}' + r'$\pm$' + f'{uncert[3]}) GeV', (0.2, 34))
     39 
     40     plt.savefig('fit_single.png')
     41 
     42     print(*popt, sep='\n')
     43 
     44 
     45 if __name__ == "__main__":
     46     main()