plots.py (3360B)
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 ax[i].xaxis.set_ticks(np.arange(0.9,1.1,0.05)) 49 50 ax[0].plot(w_np/w_00, np.abs(G_np(w_np))**2/G_n, c='black') 51 ax[0].set_xlabel(r'$\frac{\omega}{\omega_0}$') 52 ax[0].set_ylabel(r'$\frac{|G(\omega)|^2}{|G(\omega_0)|^2}$') 53 ax[0].scatter(1, 1/(gamma*w_00)**2/G_n, c='r') 54 ax[0].scatter(a_1/w_00, f_1/G_n, c='r') 55 ax[0].scatter(a_2/w_00, f_2/G_n, c='r') 56 ax[0].plot(np.linspace(w_00-gamma/2, w_00+gamma/2, 20)/w_00, f_1*np.ones(20)/G_n) 57 58 ax[1].plot(w_np/w_00, np.angle(G_np(w_np))/np.pi, c='black') 59 ax[1].set_xlabel(r'$\frac{\omega}{\omega_0}$') 60 ax[1].set_ylabel(r'$arg(G(\omega)\frac{1}{\pi}$') 61 62 fig.tight_layout() 63 64 fig.savefig('section2.png') 65 66 67 # CHAPTER 5 68 M_rho = 0.77 69 G_rho = 0.15 70 M_pi = 0.14 71 72 def F_BW(s): 73 sigma = lambda x: np.sqrt(1 - 4*M_pi**2/x) 74 G = G_rho* s/M_rho**2* (sigma(s)/sigma(M_rho**2))**3 * np.heaviside(s- 4*M_pi**2, 0 ) 75 return M_rho**2 / (M_rho**2 - s - 1j*M_rho*G) 76 77 s = np.linspace(0.1, 1, 200) 78 delta = lambda x: np.angle(F_BW(x)) 79 80 # P.V. 81 s_0 = 4*M_pi**2 82 def integrand(s_, x): 83 return (delta(s_) - delta(x))/(s_*(s_ - x)) 84 def integral(x): 85 return quad(integrand, s_0, np.inf, args=(x))[0] 86 87 I = np.vectorize(integral) 88 89 first = s/np.pi * I(s) # first integral 90 second = delta(s) * 1/np.pi * np.log(s_0/(s-s_0)) # second integral 91 92 F = np.exp(first + second + 1j * delta(s)) 93 94 # BW - Omnès representations 95 96 fig, ax = plt.subplots(figsize=[10, 7]) 97 ax.spines['right'].set_visible(False) 98 ax.spines['top'].set_visible(False) 99 ax.yaxis.set_ticks_position('left') 100 ax.xaxis.set_ticks_position('bottom') 101 102 ax.plot(s, np.abs(F), label='Omnès', c='black') 103 ax.plot(s, np.abs(F_BW(s)), label='Breit-Wigner', c='red') 104 ax.legend(loc='best', fontsize=12) 105 ax.set_xlabel(r'$s\ [GeV^2]$') 106 ax.set_ylabel(r'$|F_\pi^V(s)|$') 107 fig.savefig('omnes_bw.png') 108 109 110 if __name__ == '__main__': 111 main()