-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathsetstrongcoupl.f
270 lines (256 loc) · 8.8 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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
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
if( whichpdfpk().eq.'lha') then
continue
elseif( whichpdfpk().eq.'mlm') then
continue
else
write(*,*) ' unimplemented pdf package',whichpdfpk()
call exit(-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'
include 'pwhg_pdf.h'
real * 8 b0,mu0sq,as,pwhg_alphas
external pwhg_alphas
b0=(33-2*5)/(12*pi)
mu0sq=(2*st_lambda5MSB)**2
ccccccccccccc
c !: 20-05-2016: Improvement: rather than freezing CMW alphas, in
c order to avoid that aCMW/alphas0 exceeds 1 (which creates an
c upper-bound violation when generating ISR), it's enough to just
c increase a bit the scale at which aCMW (computed starting from the
c running of LHAPDF) is matched to alphas0. This scale is mu0sq.
c Using 4*lambda rather than 2*lambda was found empirically. Notice
c that this should not affect physics result, since the cutoff on
c radiation is above mu0sq, i.e. rad_ptsqmin > mu0sq.
if(pdf_alphas_from_pdf) then
mu0sq=(4*st_lambda5MSB)**2
else
mu0sq=(2*st_lambda5MSB)**2
endif
c Moreover, pwhg_alphas0 can be called always with nlc=5, so
c changing nlc in this file and in gen_radiation.f is not needed
c anymore. Recall that alphas0 is just used as a function that
c should be bigger than aCMW through all the pt region that can be
c probed when generating ISR (BOX paper, E.2). For single top, I
c have kept these changes, but there was no real reason to do so.
c Notice that in this way we can reproduce exactly what we run for the
c WWJ-MiNLO paper, as well as for the WW@NNLOPS paper
ccccccccccc
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'
include 'pwhg_pdf.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
real *8 alphasfrompdf
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(pdf_alphas_from_pdf) then
pwhg_alphas=alphasfrompdf(q)
else
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
endif
return
end