-
Notifications
You must be signed in to change notification settings - Fork 0
/
parton_emission_rates.py
112 lines (81 loc) · 3.06 KB
/
parton_emission_rates.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
import numpy as np
import scipy #.optimize
import scipy.interpolate
import scipy.integrate
hbarc=0.1973
######################################################
############## Parton energy loss rates ##############
######################################################
class energy_loss_rates:
def __init__(self, alpha_s, N_f, mD_factor=1, K_factor_fct_elastic=lambda T: 1, K_factor_fct_inel=lambda T: 1):
self.alpha_s = alpha_s
self.g_s = np.sqrt(4*np.pi*alpha_s)
self.N_f = N_f
self.N_c = 3
self.C_A = self.N_c
self.sqrt_Nc_Nf_factor=np.sqrt(self.N_c/3.+self.N_f/6)
self.mD_factor=mD_factor
self.K_factor_fct_elastic=K_factor_fct_elastic
self.K_factor_fct_inel=K_factor_fct_inel
def m_D(self, T):
factor=np.sqrt(self.N_c/3.+self.N_f/6)
return self.mD_factor*self.g_s*T*factor
# For gluons
def qhat_eff(self, omega,T):
# Solution of qhat_eff = q_soft^22(mu_T) with mu_T^2=C0 Sqrt(2 omega qhat_eff)
# with C_0=2 e^{2-\gamma +\frac{\pi }{4}}
# and q_soft^22(mu_T)=norm Log(mu_T^2/m_D^2)
# with norm=alpha_s*C_A*T*m_D(T)**2
mD=self.m_D(T)
mD2=mD*mD
mD4=mD2*mD2
norm=self.alpha_s*self.C_A*T*mD2
C_0=18.1983
return -0.5*norm*np.real(scipy.special.lambertw(-mD4/(norm*C_0**2*omega),k=-1))
#print(qhat_eff(1e-2*1000,.3))
# should be about 0.052 with N_f=0 and alpha_s=0.1
def dGamma_domega_inel(self, p, omega,T):
if (p == 0)or(omega>=p):
res=0.0
else:
one_over_p=1./p
one_over_pi=0.318309886183790671537767526745
#z=omega/p
z=one_over_p*omega
qhat_eff_val=self.qhat_eff(omega,T)
#qhat_eff_val=0.04 #qhat_eff(1e-2*p,T)
#print(qhat_eff_val)
#res=self.alpha_s*self.N_c/(np.pi*p)*np.power(z*(1-z),-3./2.)*np.sqrt(qhat_eff_val/p)
res=self.alpha_s*self.N_c*one_over_p*one_over_pi*np.sqrt(qhat_eff_val*one_over_p/(z*z*z*(1-z)*(1-z)*(1-z)))
return res
# From https://arxiv.org/pdf/1006.2379.pdf
def dGamma_domega_elastic(self,p,omega,T):
return self.alpha_s*self.C_A*self.m_D(T)**2/(4.*omega*omega)
def dGamma_domega(self, p,omega,T):
return self.K_factor_fct_inel(T)*self.dGamma_domega_inel(p,omega,T)+self.K_factor_fct_elastic(T)*self.dGamma_domega_elastic(p,omega,T)
#plt.figure()
#plt.axes().xaxis.set_minor_locator(AutoMinorLocator())
#plt.gca().yaxis.set_ticks_position('both')
#plt.axes().yaxis.set_minor_locator(AutoMinorLocator())
#
#plt.xscale('log')
#plt.yscale('log')
#plt.xlim(3e-4,5e-1)
#plt.ylim(1e-6,1e-1)
#
#plt.xlabel(r'$\omega/p$')
#plt.ylabel(r"$d\Gamma/d\omega$")
#
## Plot the exact solution
#z_range=np.array([10.**x for x in np.arange(-4,0,.1)])
#p=1000.
#T=0.3
#omega_range=z_range*p
#dGamma_list=[dGamma_domega_inel(p,omega,T) for omega in omega_range]
#
#plt.plot(z_range, dGamma_list,"-",color='red',label="", linewidth=3)
#
#plt.legend()
#plt.tight_layout()
##plt.savefig("T_profile.pdf")
#plt.show()