plots.py (3122B)
1 import numpy as np 2 import matplotlib.pyplot as plt 3 import matplotlib 4 from sympy import * 5 from scipy.integrate import quad 6 7 font = {'family' : 'sans', 8 'weight' : 'bold', 9 'size' : 18} 10 11 matplotlib.rc('font', **font) 12 matplotlib.matplotlib_fname() 13 14 def main(): 15 # numpy representation 16 w_00 = 1 17 gamma = w_00/20 18 G_np = lambda w: 1/(-w**2 - 1j*gamma*w + w_00**2) 19 w_np = np.linspace(w_00-2*gamma, w_00+2*gamma, 200) 20 21 # sympy representation 22 w = Symbol('w', real=True) 23 z = Symbol('z') 24 g = Symbol('g', real=True) 25 w_0 = Symbol('w_0', real=True) 26 27 G = 1/(-w**2 - 1j*g*w + w_0**2) 28 G_n = np.abs(G.subs([(w, w_00), (w_0, w_00), (g, gamma)]))**2 29 30 #equation for half maximum solve for w 31 solutions = solve(Eq(1/2*1/(g*w_0)**2, abs(G)**2), w) 32 33 so = solve(Eq(0, 1/G.subs(w, z)), z) 34 35 a_1 = solutions[1].subs([(w_0, w_00), (g, gamma)]) 36 f_1 = abs(G.subs([(w, a_1), (g, gamma), (w_0, w_00)]))**2 37 38 a_2 = solutions[3].subs([(w_0, w_00), (g, gamma)]) 39 f_2 = abs(G.subs([(w, a_2), (g, gamma), (w_0, w_00)]))**2 40 fig, ax= plt.subplots(1, 2, figsize=[17,7]) 41 42 # Plots for |G(w)|^2 and arg(G(w)) 43 for i in range(len(ax)): 44 ax[i].spines['right'].set_visible(False) 45 ax[i].spines['top'].set_visible(False) 46 ax[i].yaxis.set_ticks_position('left') 47 ax[i].xaxis.set_ticks_position('bottom') 48 49 ax[0].plot(w_np/w_00, np.abs(G_np(w_np))**2/G_n, c='black') 50 ax[0].set_xlabel(r'$\frac{\omega}{\omega_0}$') 51 ax[0].set_ylabel(r'$\frac{|G(\omega)|^2}{|G(\omega_0)|^2}$') 52 ax[0].scatter(1, 1/(gamma*w_00)**2/G_n, c='r') 53 ax[0].scatter(a_1/w_00, f_1/G_n, c='r') 54 ax[0].scatter(a_2/w_00, f_2/G_n, c='r') 55 ax[0].plot(np.linspace(w_00-gamma/2, w_00+gamma/2, 20)/w_00, f_1*np.ones(20)/G_n) 56 57 ax[1].plot(w_np/w_00, np.angle(G_np(w_np))/np.pi, c='black') 58 ax[1].set_xlabel(r'$\frac{\omega}{\omega_0}$') 59 ax[1].set_ylabel(r'$arg(G(\omega)\frac{1}{\pi}$') 60 61 fig.tight_layout() 62 63 fig.savefig('section2.png') 64 65 66 # CHAPTER 5 67 M_rho = 0.77 68 G_rho = 0.15 69 M_pi = 0.14 70 71 def F_BW(s): 72 sigma = lambda x: np.sqrt(1 - 4*M_pi**2/x) 73 G = G_rho* s/M_rho**2* (sigma(s)/sigma(M_rho**2))**3 * np.heaviside(s- 4*M_pi**2, 0 ) 74 return M_rho**2 / (M_rho**2 - s - 1j*M_rho*G) 75 76 s = np.linspace(0.1, 1, 200) 77 delta = lambda x: np.angle(F_BW(x)) 78 79 # P.V. 80 s_0 = 4*M_pi**2 81 def integrand(s_, x): 82 return (delta(s_) - delta(x))/(s_*(s_ - x)) 83 def integral(x): 84 return quad(integrand, s_0, np.inf, args=(x))[0] 85 86 I = np.vectorize(integral) 87 88 first = s/np.pi * I(s) # first integral 89 second = delta(s) * 1/np.pi * np.log(s_0/(s-s_0)) # second integral 90 91 F = np.exp(first + second + 1j * delta(s)) 92 93 # BW - Omnes representations 94 plt.figure(figsize=[10, 7]) 95 plt.plot(s, np.abs(F), label='Omnes', c='black') 96 plt.plot(s, np.abs(F_BW(s)), label='Breit-Wigner', c='red') 97 plt.legend(loc='best', fontsize=12) 98 plt.xlabel(r'$s\ [GeV^2]$') 99 plt.ylabel(r'$|F_\pi^V(s)|$') 100 plt.savefig('omnes_bw.png') 101 102 103 if __name__ == '__main__': 104 main()