-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathALLM.py
More file actions
160 lines (103 loc) · 3.74 KB
/
ALLM.py
File metadata and controls
160 lines (103 loc) · 3.74 KB
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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
# Final Version -- October 2023 -- Hamzeh Khanpour
import math
import scipy.integrate as integ
# units in GeV
pmass = 0.938272081
pi0mass = 0.1349768
# ALLM parameters -- arXiv:hep-ph/9712415
Mass2_0 = 0.31985
Mass2_P = 49.457
Mass2_R = 0.15052
Q2_0 = 0.52544
Lambda2 = 0.06527
Ccp = (0.28067, 0.22291, 2.1979)
Cap = (-0.0808, -0.44812, 1.1709)
Cbp = (0.36292, 1.8917, 1.8439)
Ccr = (0.80107, 0.97307, 3.4942)
Car = (0.58400, 0.37888, 2.6063)
Cbr = (0.01147, 3.7582, 0.49338)
# --------------------------------------------------------------
def tvalue(Q2):
return math.log \
((math.log((Q2 + Q2_0) / Lambda2) / math.log(Q2_0 / Lambda2)))
def xP(xbj, Q2):
if xbj == 0:
print("xbj zero")
return -1.
xPinv = 1. + Q2 / (Q2 + Mass2_P) * (1. / xbj - 1.)
return 1. / xPinv
def xR(xbj, Q2):
if xbj == 0:
print("xbj zero")
return -1.
xPinv = 1. + Q2 / (Q2 + Mass2_R) * (1. / xbj - 1.)
return 1. / xPinv
def type1(tval, tuple1):
return tuple1[0] + tuple1[1] * (tval ** tuple1[2])
def type2(tval, tuple1):
return tuple1[0] +\
(tuple1[0] - tuple1[1]) * (1. / (1. + tval ** tuple1[2]) - 1.)
def aP(tval):
return type2(tval, Cap)
def bP(tval):
return type1(tval, Cbp)
def cP(tval):
return type2(tval, Ccp)
def aR(tval):
return type1(tval, Car)
def bR(tval):
return type1(tval, Cbr)
def cR(tval):
return type1(tval, Ccr)
def allm_f2P(xbj, Q2):
tval = tvalue(Q2)
return cP(tval) * (xP(xbj, Q2) ** aP(tval)) * ((1. - xbj) ** bP(tval))
def allm_f2R(xbj, Q2):
tval = tvalue(Q2)
return cR(tval) * (xR(xbj, Q2) ** aR(tval)) * ((1. - xbj) ** bR(tval))
def allm_f2(xbj, Q2):
return Q2 / (Q2 + Mass2_0) * (allm_f2P(xbj, Q2) + allm_f2R(xbj, Q2))
# --------------------------------------------------------------
def allm_f2divx_mN(mN, Q2, yp):
A2 = mN*mN - pmass * pmass # Hamzeh
mqdiff = mN*mN - pmass * pmass + Q2 # Hamzeh
if mqdiff < 0:
print('mN*mN, Q2:', mN*mN, Q2)
return 0.
xbj = Q2 / mqdiff
minM2 = (pmass + pi0mass) * (pmass + pi0mass)
qmin2 = (mN*mN / (1.0 - yp) - pmass * pmass) * yp
if xbj < 0:
print('xbj: ', xbj)
return 0.
else:
# 27 Jul 2021: adding Qmin2
if qmin2 < Q2:
return allm_f2(xbj, Q2) / Q2**0.0 * 2.0 * mN * mqdiff # Hamzeh: It should be Q2**2.0 in Syy200.py
else:
return 0.
def allm_formM_mN2(Q2, yp, mMin2, mNmax):
return integ.quad(allm_f2divx_mN, mMin2, mNmax, args=(Q2, yp),
epsrel=1.e-2)
# --------------------------------------------------------------
def allm_xf2_mN(mN, Q2, yp):
A2 = mN*mN - pmass * pmass # Hamzeh
mqdiff = mN*mN - pmass * pmass + Q2
if mqdiff < 0:
print('mN*mN, Q2:', mN*mN, Q2)
return 0.
xbj = Q2 / mqdiff
minM2 = (pmass + pi0mass) * (pmass + pi0mass)
qmin2 = (mN*mN / (1.0 - yp) - pmass * pmass) * yp
if xbj < 0:
print('xbj: ', xbj)
return 0.
else:
if qmin2 < Q2:
return allm_f2(xbj, Q2) / Q2**0.0 * 2.0 * mN * (1.0/mqdiff) * ( 1.0 - qmin2 / Q2 )
else:
return 0.
def allm_formE_qmin2(Q2, yp, mMin2, mNmax):
return integ.quad(allm_xf2_mN, mMin2, mNmax, args=(Q2, yp),
epsrel=1.e-2)
# --------------------------------------------------------------