tprak

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

multiple_fit.py (3667B)

      1 
      2 import numpy as np
      3 import sympy as sp
      4 from scipy.linalg import block_diag
      5 from matplotlib import pyplot as plt
      6 from scipy.special import gammaincc
      7 
      8 from method import *
      9 from model import *
     10 
     11 import matplotlib.pylab as pylab
     12 params = {'legend.fontsize': 'x-large',
     13           'figure.figsize': (15, 5),
     14          'axes.labelsize': 'x-large',
     15          'axes.titlesize':'x-large',
     16          'xtick.labelsize':'x-large',
     17          'ytick.labelsize':'x-large'}
     18 pylab.rcParams.update(params)
     19 
     20 global p0;
     21 p0 = [0.9, 0.2, 0.81, 0.04, 0.02, -1, 0.84, 1.55]   # in GeV
     22 
     23 def main():
     24     # Preparation
     25     SND_data = np.loadtxt('../data/SND-VFF.txt')
     26     SND_x = SND_data[:, 0]
     27     SND_y = SND_data[:, 1]
     28     SND_statcov = (np.eye(len(SND_x)) * SND_data[:, 2])**2
     29     SND_relsystcov = (np.eye(len(SND_x)) * \
     30                       np.array([0.032 if i<=2 else 0.013 for i in range(len(SND_x))]))**2
     31 
     32     CMD2_data = np.loadtxt('../data/SND-VFF.txt')
     33     CMD2_x = CMD2_data[:, 0]
     34     CMD2_y = CMD2_data[:, 1]
     35     CMD2_statcov = (np.eye(len(CMD2_x)) * CMD2_data[:, 2])**2
     36     CMD2_relsystcov = []
     37     for i in range(len(CMD2_x)):
     38         if i<=43:
     39             CMD2_relsystcov.append(0.006)
     40         elif i>10 and i<=53:
     41             CMD2_relsystcov.append(0.007)
     42         else:
     43             CMD2_relsystcov.append(0.008)
     44     CMD2_relsystcov = (np.eye(len(CMD2_x)) * np.array(CMD2_relsystcov))**2
     45 
     46     KLOE_data = np.loadtxt('../data/KLOE-VFF.txt')
     47     KLOE_x = KLOE_data[:, 0]
     48     KLOE_y = KLOE_data[:, 1]
     49     KLOE_statcov = np.loadtxt('../data/KLOE-StatCov.txt')
     50     KLOE_relsystcov = np.loadtxt('../data/KLOE-RelSystCov.txt')
     51 
     52     BABAR_data = np.loadtxt('../data/BABAR-VFF.txt')
     53     BABAR_x = BABAR_data[:, 0]
     54     BABAR_y = BABAR_data[:, 1]
     55     BABAR_statcov = np.loadtxt('../data/BABAR-StatCov.txt')
     56     BABAR_relsystcov = np.loadtxt('../data/BABAR-RelSystCov.txt')
     57 
     58     x_data = np.block([SND_x, CMD2_x, KLOE_x, BABAR_x])
     59     y_data = np.block([SND_y, CMD2_y, KLOE_y, BABAR_y])
     60     cov_stat = block_diag(SND_statcov, CMD2_statcov, KLOE_statcov, BABAR_statcov)
     61     cov_relsyst = block_diag(SND_relsystcov, CMD2_relsystcov, KLOE_relsystcov, BABAR_relsystcov)
     62     var_str = "m_q g_q m_w g_w e_w a b c"
     63 
     64     # Fit
     65     p, dp, chi_sq = t0_fit(model, var_str, x_data, y_data, p0, cov_stat,\
     66                     cov_relsyst)
     67 
     68 
     69     dof = len(x_data) - len(p)
     70     p_value = 1 - gammaincc(chi_sq[-1]/2, dof/2)
     71     chi_min = chi_sq[-1]/dof
     72 
     73     with open(f'./results/multi.txt', 'w+') as f:
     74         for i in range(len(p)):
     75             f.write(f'{p[i]}\t\\pm {dp[i]}\n')
     76 
     77         f.write(f'chi_sq_min_dof = {chi_min}\n')
     78         f.write(f'p_value = {p_value}\n')
     79 
     80 
     81     # Plot
     82     x_model = np.linspace(x_data[0], x_data[-1], 500)
     83     x = sp.symbols('x')
     84     la_mod = model(x, *p)
     85     la_mod = sp.lambdify(x, la_mod, modules=['numpy', {'Heaviside': np.heaviside}])
     86     y_model = la_mod(x_model)
     87 
     88     plt.figure(figsize=[10, 7])
     89 
     90     plt.errorbar(SND_x, SND_y, yerr=np.sqrt(np.diag(SND_statcov)), fmt='.', ms=8, label='SND', c='blue', alpha=0.4)
     91     plt.errorbar(CMD2_x, CMD2_y, yerr=np.sqrt(np.diag(CMD2_statcov)), fmt='.', ms=8, label='CMD2', c='green', alpha=0.4)
     92     plt.errorbar(KLOE_x, KLOE_y, yerr=np.sqrt(np.diag(KLOE_statcov)), fmt='.', ms=8, label='KLOE', c='orange', alpha=0.4)
     93     plt.errorbar(BABAR_x, BABAR_y, yerr=np.sqrt(np.diag(BABAR_statcov)), fmt='.', ms=8, label='BABAR', c='red', alpha=0.4)
     94 
     95     plt.plot(x_model, y_model, lw=2, c='black', label='Fit')
     96 
     97     plt.title('t0-Multifit')
     98     plt.legend(loc='best')
     99     plt.ylabel(r'$|F_{\pi}^V(s)|^2|$')
    100     plt.xlabel(r'$\sqrt{s}$[GeV]')
    101 
    102     plt.savefig('./plots/multi.png')
    103     plt.close()
    104 
    105 
    106 if __name__ == "__main__":
    107     main()