-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsetstrongcoupl.f
254 lines (245 loc) · 7.94 KB
/
setstrongcoupl.f
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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
subroutine setscalesbtilde
implicit none
include 'pwhg_math.h'
include 'pwhg_st.h'
include 'pwhg_flg.h'
real * 8 pwhg_alphas
external pwhg_alphas
real * 8 muf,mur
c signal we will begin by computing Born type contributions
flg_btildepart='b'
call set_fac_ren_scales(muf,mur)
st_mufact2= muf**2*st_facfact**2
st_muren2 = mur**2*st_renfact**2
st_alpha = pwhg_alphas(st_muren2,st_lambda5MSB,st_nlight)
end
subroutine setscalesbtlreal
implicit none
include 'pwhg_math.h'
include 'pwhg_st.h'
include 'pwhg_flg.h'
real * 8 pwhg_alphas
external pwhg_alphas
real * 8 muf,mur
logical ini
data ini/.true./
save ini
real * 8 powheginput
external powheginput
if(ini) then
if(powheginput("#btlscalereal").eq.1d0) then
flg_btlscalereal=.true.
else
flg_btlscalereal=.false.
endif
ini=.false.
endif
if(flg_btlscalereal) then
c if this is active we may compute scales that depends upon
c the real kinematics; the user routine set_fac_ren_scales
c should test the flag flg_btildepart to see if this is the case
flg_btildepart='r'
call set_fac_ren_scales(muf,mur)
st_mufact2= muf**2*st_facfact**2
st_muren2 = mur**2*st_renfact**2
st_alpha = pwhg_alphas(st_muren2,st_lambda5MSB,st_nlight)
endif
end
subroutine setscalesbtlct
implicit none
include 'pwhg_math.h'
include 'pwhg_st.h'
include 'pwhg_flg.h'
real * 8 pwhg_alphas
external pwhg_alphas
real * 8 muf,mur
logical ini
data ini/.true./
save ini
real * 8 powheginput
external powheginput
if(ini) then
if(powheginput("#btlscalereal").eq.1d0) then
flg_btlscalereal=.true.
else
flg_btlscalereal=.false.
endif
if(powheginput("#btlscalect").eq.1d0) then
flg_btlscalect=.true.
else
flg_btlscalect=.false.
endif
endif
if(flg_btlscalereal.and.flg_btlscalect) then
c signal we will begin by computing counterterm contributions, in cases
c when it is desirable to have the scales of the counterterm differ from
c those of the real term (the token btlscalect selects this case)
c The user routine should test the flag flg_btildepart to see if
c we are in a counterterm.
flg_btildepart='c'
call set_fac_ren_scales(muf,mur)
st_mufact2= muf**2*st_facfact**2
st_muren2 = mur**2*st_renfact**2
st_alpha = pwhg_alphas(st_muren2,st_lambda5MSB,st_nlight)
endif
end
subroutine set_rad_scales(ptsq)
implicit none
real * 8 ptsq
include 'pwhg_math.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_st.h'
include 'pwhg_rad.h'
include 'pwhg_pdf.h'
real * 8 pwhg_alphas
integer nf
external pwhg_alphas
character * 3 whichpdfpk
external whichpdfpk
integer ini,mem
data ini/0/
save ini
if (ini.eq.0) then
if( whichpdfpk().eq.'lha') then
continue
elseif( whichpdfpk().eq.'mlm') then
continue
else
write(*,*) ' unimplemented pdf package',whichpdfpk()
stop
endif
ini=1
endif
st_mufact2=max(pdf_q2min,ptsq)
cccccccccccccccccccccccccccccccccc
c In case of final-state radiation, Born and real PDF's
c should always cancel out in the ratio R/B. If the radiation scale
c is too low, this cancellation can be spoilt because PDF's can vanish,
c typically when a heavy flavour is present as initial state.
c To prevent this, we use a scale higher than the heavy-flavour
c threshold, so that PDF's are evaluated with a safe value for
c mufact (50 is an arbitrary choice).
if(rad_kinreg.ge.2) st_mufact2=50.**2
st_muren2=ptsq
st_alpha = pwhg_alphas(st_muren2,st_lambda5MSB,-1)
if(st_muren2.lt.rad_charmthr2) then
nf=3
elseif(st_muren2.lt.rad_bottomthr2) then
nf=4
else
nf=5
endif
st_alpha = st_alpha *
# (1+st_alpha/(2*pi)*((67d0/18-pi**2/6)*ca-5d0/9*nf))
end
subroutine init_rad_lambda
implicit none
include 'pwhg_math.h'
include 'pwhg_st.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_rad.h'
real * 8 b0,mu0sq,as,pwhg_alphas
external pwhg_alphas
b0=(33-2*5)/(12*pi)
mu0sq=(2*st_lambda5MSB)**2
c running value of alpha at initial scale (see notes: running_coupling)
as=pwhg_alphas(mu0sq,st_lambda5MSB,-1)
c for better NLL accuracy (FNO2006, (4.32) and corresponding references)
as=as*(1+as/(2*pi)*((67d0/18-pi**2/6)*ca-5d0/9*3))
rad_lamll=sqrt(exp(-1/(b0*as))*mu0sq)
end
function pwhg_alphas0(q2,xlam,inf)
implicit none
real * 8 pwhg_alphas0,q2,xlam
integer inf
real * 8 pi
parameter (pi=3.141592653589793d0)
real * 8 b0
b0=(33-2*inf)/(12*pi)
pwhg_alphas0=1/(b0*log(q2/xlam**2))
end
C----------------------------------------------------------------------------
C-------------------------------------------------------------------
C------- ALPHA QCD -------------------------------------
c Program to calculate alfa strong with nf flavours,
c as a function of lambda with 5 flavors.
c The value of alfa is matched at the thresholds q = mq.
c When invoked with nf < 0 it chooses nf as the number of
c flavors with mass less then q.
c
function pwhg_alphas(q2,xlam,inf)
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_rad.h'
real * 8 pwhg_alphas,q2,xlam
integer inf
real * 8 pi
parameter (pi=3.141592653589793d0)
real * 8 olam,b5,bp5,b4,bp4,b3,bp3,xlc,xlb,xllc,xllb,c45,c35,
# xmc,xmb
real * 8 q,xlq,xllq
integer nf
data olam/0.d0/
save olam,b5,bp5,b4,bp4,b3,bp3,xlc,xlb,xllc,xllb,c45,c35,xmc,xmb
c logical ini
c data ini/.true./
c save ini
c pwhg_alphas = 30d0
c if (ini) then
c write(*,*) '****************************************'
c write(*,*) '****************************************'
c write(*,*) ' RETURN alpha_s = ',pwhg_alphas
c write(*,*) '****************************************'
c write(*,*) '****************************************'
c ini = .false.
c endif
c return
if(xlam.ne.olam) then
olam = xlam
xmc=sqrt(rad_charmthr2)
xmb=sqrt(rad_bottomthr2)
b5 = (33-2*5)/pi/12
bp5 = (153 - 19*5) / pi / 2 / (33 - 2*5)
b4 = (33-2*4)/pi/12
bp4 = (153 - 19*4) / pi / 2 / (33 - 2*4)
b3 = (33-2*3)/pi/12
bp3 = (153 - 19*3) / pi / 2 / (33 - 2*3)
xlc = 2 * log(xmc/xlam)
xlb = 2 * log(xmb/xlam)
xllc = log(xlc)
xllb = log(xlb)
c45 = 1/( 1/(b5 * xlb) - xllb*bp5/(b5 * xlb)**2 )
# - 1/( 1/(b4 * xlb) - xllb*bp4/(b4 * xlb)**2 )
c35 = 1/( 1/(b4 * xlc) - xllc*bp4/(b4 * xlc)**2 )
# - 1/( 1/(b3 * xlc) - xllc*bp3/(b3 * xlc)**2 ) + c45
endif
q = sqrt(q2)
xlq = 2 * log( q/xlam )
xllq = log( xlq )
nf = inf
if( nf .lt. 0) then
if( q .gt. xmb ) then
nf = 5
elseif( q .gt. xmc ) then
nf = 4
else
nf = 3
endif
endif
if ( nf .eq. 5 ) then
pwhg_alphas = 1/(b5 * xlq) - bp5/(b5 * xlq)**2 * xllq
elseif( nf .eq. 4 ) then
pwhg_alphas =
# 1/( 1/(1/(b4 * xlq) - bp4/(b4 * xlq)**2 * xllq) + c45 )
elseif( nf .eq. 3 ) then
pwhg_alphas =
# 1/( 1/(1/(b3 * xlq) - bp3/(b3 * xlq)**2 * xllq) + c35 )
else
print *,'error in alfa: unimplemented # of light flavours',nf
call exit(1)
endif
return
end