« back to all changes in this revision
Viewing changes to setlocalscales.f
- browse files at revision 116
- compare with another revision
- download diff
- download tarball
- view history from revision 116
- Committer: Rikkert Frederix
- Date: 2017-06-19 18:34:45 UTC
- mfrom: (112.1.26 PWHG_STJ)
- Revision ID: frederix@physik.uzh.ch-20170619183445-5832xjj3q30gdr9f
merge with latest PWHG_STJ
- files added:
- auxiliary-progress.f
- reweight
- reweight/obscorrelations
- files modified:
- Born_phsp.f
added
removed
setlocalscales.f
1
1
2
subroutine setlocalscales(iuborn,imode,rescfac)
2
3
c returns the rescaling factor including sudakov form factors and
3
4
c coupling rescaling, for born (imode=1) and NLO corrections (imode=2)
60
61
endif
61
62
end
62
63
64
C at the moment momenta routine works for 5 incoming momenta (not 6)
63
65
subroutine setlocalscales0
64
66
1 (flav,pin,basicfac,bornfac,nlofac)
65
67
use auxiliary
72
74
real * 8 pin(0:3,nlegborn),basicfac,bornfac,nlofac
73
75
integer onem
74
76
parameter (onem=1000000)
75
real * 8 scales(nlegborn),p(0:3,nlegborn),ptot(0:3),ptot2(0:3),
77
real * 8 scales(nlegborn),p(0:3,nlegborn),ptot(0:3),
76
78
1 lscalej,lscalek
77
79
integer j,k,l,lflav(nlegborn),jmerge,kmerge,inlofac,npart
78
80
integer mergedfl
79
real * 8 q2merge,q2merge0,renfac2,facfact2,alphas,mu2,muf2,q2om2
81
real * 8 q2merge,renfac2,facfact2,alphas,mu2,muf2,q2om2
80
82
real * 8 sudakov,expsudakov,pwhg_alphas,b0,powheginput
81
83
$ ,sudakov_massive,expsudakov_massive,sudakov_exp
82
84
external sudakov,expsudakov,pwhg_alphas,powheginput
83
85
$ ,sudakov_massive,expsudakov_massive
84
real * 8 q2mergeMAX,single_top_scale(4)
85
logical dijetflag
86
common/cdijetflag/dijetflag
87
logical raisingscales,ini
88
save raisingscales,ini
89
data ini/.true./
86
real * 8 q2mergeMAX,q2hard
90
87
logical pwhg_isfinite
91
integer n_badq2merge,ixx
92
data n_badq2merge/0/
93
save n_badq2merge
88
real*8 qMom(0:3),qprMom(0:3),bMom(0:3),tMom(0:3),WMom(0:3),jmom(0:3)
89
real*8 rad1Mom(0:3),rad2Mom(0:3),totalMom(0:3)
90
integer nUnClusteredPart,nClusteredPart
91
real*8 m2top,Q2W,two_pt_dot_pb
92
real *8 q2merge_new,bMom_new(0:3),tMom_new(0:3),WMom_new(0:3)
94
93
95
if(ini) then
96
if(powheginput("#raisingscales").eq.0) then
97
raisingscales = .false.
98
else
99
raisingscales = .true.
100
endif
101
ini = .false.
102
endif
94
C - Initialization of dynamical parameters
103
95
renfac2=st_renfact**2
104
96
facfact2=st_facfact**2
105
lflav=flav
106
p=pin
107
97
scales=0
98
sudakov_exp=0d0
99
100
C - Make sure momenta and flavour arrays have no residual junk in them
101
p=0d0
102
lflav=0
103
104
C - Initialize kinematics
108
105
q2mergeMAX=-1d10
109
sudakov_exp=0d0
110
do l=1,nlegborn
111
call findNearestNeighbours(p,lflav,jmerge,kmerge,mergedfl,
112
$ q2merge)
113
if(.not.pwhg_isfinite(q2merge)) then
114
q2merge=1d10
115
n_badq2merge=n_badq2merge+1
116
if((n_badq2merge.le.100).or.mod(n_badq2merge,100).eq.0) then
117
write(*,*) 'setlocalscales.f:'
118
write(*,*) 'bad q2merge from findNearestNeighbours'
119
write(*,*) 'resetting q2merge=1d10'
120
write(*,*) 'n_badq2merge = ',n_badq2merge
121
write(*,*) 'q2merge = ', q2merge
122
write(*,*) 'lflav = ',lflav
123
write(*,*) 'jmerge = ',jmerge
124
write(*,*) 'kmerge = ',kmerge
125
write(*,*) 'mergedfl = ',mergedfl
126
write(*,*) 'momenta = '
127
do ixx=3,nlegborn
128
write(*,*) ixx,p(:,ixx)
129
enddo
130
endif
131
if(n_badq2merge.eq.100) then
132
write(*,*) 'Further bad q2 merge messages are suppressed'
133
endif
134
endif
135
if(q2merge.lt.1d10) then
136
c perform the merging
137
if(q2merge.gt.q2mergeMAX) q2mergeMAX=q2merge
138
lscalej=scales(jmerge)
139
lscalek=scales(kmerge)
140
scales(jmerge)=q2merge
141
if(lscalej.eq.0) then
142
c This is the first merge; it sets the low scale for
143
c all partons; no Sudakov factor or reweighting is introduced
144
do j=1,nlegborn
145
scales(j)=q2merge
146
enddo
147
c save this scale; it is the Q_0 scale that appears in all Sudakovs
148
q2merge0=q2merge
149
bornfac=0
150
c Provide alpha_S reweighting for the first merge
151
alphas=pwhg_alphas(max(q2merge*renfac2,1d0),
152
1 st_lambda5MSB,st_nlight)
153
basicfac=alphas/st_alpha
154
nlofac=basicfac
155
mu2=max(q2merge*renfac2,1d0)
156
muf2=max(q2merge*facfact2,1d0)
157
inlofac=1
158
else
159
c provide Sudakov
160
sudakov_exp=sudakov_exp+
161
1 sudakov(q2merge0,q2merge,lscalej,lflav(jmerge))
162
sudakov_exp=sudakov_exp+
163
1 sudakov(q2merge0,q2merge,lscalek,lflav(kmerge))
164
bornfac=bornfac+
165
1 expsudakov(q2merge0,q2merge,lscalej,lflav(jmerge))
166
bornfac=bornfac+
167
1 expsudakov(q2merge0,q2merge,lscalek,lflav(kmerge))
168
c provide alpha_S reweighting
169
alphas=pwhg_alphas(max(q2merge*renfac2,1d0),
170
1 st_lambda5MSB,st_nlight)
171
basicfac=basicfac*alphas/st_alpha
172
mu2=mu2*max(q2merge*renfac2,1d0)
173
nlofac=nlofac+alphas/st_alpha
174
inlofac=inlofac+1
175
endif
176
if(jmerge.gt.2) then
177
p(:,jmerge)=p(:,jmerge)+p(:,kmerge)
178
else
179
p(3,jmerge)=p(3,jmerge)-p(3,kmerge)
180
p(0,jmerge)=abs(p(3,jmerge))
181
endif
182
lflav(kmerge)=onem
183
lflav(jmerge)=mergedfl
184
else
185
goto 99
186
endif
187
enddo
188
99 continue
189
190
c No more merging is possible.
191
192
if(.not.dijetflag) then
193
c Define the initial scale as
194
c the invariant mass of the remaining system
195
ptot=0
196
do j=3,nlegborn
197
if(lflav(j).ne.onem) then
198
ptot=ptot+p(:,j)
199
endif
200
enddo
201
if(raisingscales) then
202
q2merge=max(q2mergeMAX,
203
$ ptot(0)**2-ptot(1)**2-ptot(2)**2-ptot(3)**2)
204
else
205
q2merge=ptot(0)**2-ptot(1)**2-ptot(2)**2-ptot(3)**2
206
endif
207
208
c SINGLE TOP SPECIAL: first compute -Q^2 (W-boson t-channel)
209
npart=0
210
do j=1,nlegborn
211
if (lflav(j).ne.onem) npart=npart+1
212
enddo
213
if (npart.eq.4) then
214
if (abs(lflav(1)).eq.5) then
215
ptot=p(:,1)
216
ptot2=p(:,1)
217
else
218
ptot=p(:,2)
219
ptot2=p(:,2)
220
endif
221
do j=3,nlegborn
222
if (abs(lflav(j)).eq.6) then
223
ptot=ptot-p(:,j)
224
ptot2=ptot2+p(:,j)
225
endif
226
enddo
227
elseif (npart.eq.5) then ! must have initial state gluon
228
if (lflav(1).ne.0 .and. lflav(2).ne.0) then
229
write (*,*) 'ERROR #2 in setting cluster scales',npart
230
write (*,*) lflav
231
stop
232
endif
233
if (abs(lflav(1)).ne.5 .and. abs(lflav(2)).ne.5) then
234
c if b is final state, gluon is attached to b-top line
235
if (lflav(1).eq.0) then
236
ptot=p(:,1)
237
ptot2=p(:,1)
238
else
239
ptot=p(:,2)
240
ptot2=p(:,2)
241
endif
242
do j=3,nlegborn
243
if (abs(lflav(j)).eq.5 .or. abs(lflav(j)).eq.6) then
244
ptot=ptot-p(:,j)
245
ptot2=ptot2+p(:,j)
246
endif
247
enddo
248
else
249
c if b is initial state, gluon is attached to massless quark line
250
if (lflav(1).ne.0) then
251
ptot=p(:,1)
252
ptot2=p(:,1)
253
else
254
ptot=p(:,2)
255
ptot2=p(:,2)
256
endif
257
do j=3,nlegborn
258
if (abs(lflav(j)).eq.6) then
259
ptot=ptot-p(:,j)
260
ptot2=ptot2+p(:,j)
261
endif
262
enddo
263
endif
264
else
265
write (*,*) 'ERROR #1 in setting cluster scales',npart
266
write (*,*) lflav
267
stop
268
endif
269
270
271
! top mass squared:
272
do j=3,nlegborn
273
if (abs(lflav(j)).eq.6) then
274
single_top_scale(4)=
275
& p(0,j)**2-p(1,j)**2-p(2,j)**2-p(3,j)**2
276
endif
277
enddo
278
! typical hard scale for massive quark line:
279
C p_w^2
280
single_top_scale(1)=
281
& abs(ptot(0)**2-ptot(1)**2-ptot(2)**2-ptot(3)**2)+
282
& single_top_scale(4)
283
! typical hard scale for massless quark line:
284
single_top_scale(2)=
285
& abs(ptot(0)**2-ptot(1)**2-ptot(2)**2-ptot(3)**2)
286
! Hard scale in mass logarithm:
287
single_top_scale(3)=single_top_scale(1)
288
! typical scale for this process:
289
q2merge=(single_top_scale(1)+single_top_scale(2))/2d0
290
else
291
c Dijet case: use the scalar sum of the pt of the two partons
292
q2merge = 0
293
do j=3,nlegborn
294
if(lflav(j).ne.onem) then
295
q2merge=sqrt(p(1,j)**2+p(2,j)**2)+q2merge
296
endif
297
enddo
298
q2merge=q2merge**2
299
if(raisingscales) then
300
q2merge=max(q2mergeMAX,q2merge)
301
endif
302
endif
303
304
if(scales(1).gt.0) then
305
do j=1,nlegborn
306
if(abs(lflav(j)).lt.st_nlight) then
307
sudakov_exp=sudakov_exp+sudakov(q2merge0
308
$ ,single_top_scale(2),scales(j),lflav(j))
309
bornfac=bornfac+ expsudakov(q2merge0,single_top_scale(2)
310
$ ,scales(j),lflav(j))
311
elseif(abs(lflav(j)).eq.st_nlight) then
312
sudakov_exp=sudakov_exp+sudakov(q2merge0
313
$ ,single_top_scale(1),scales(j),lflav(j))
314
bornfac=bornfac+ expsudakov(q2merge0,single_top_scale(1)
315
$ ,scales(j),lflav(j))
316
elseif(abs(lflav(j)).gt.st_nlight .and. abs(lflav(j)).le.6)
317
$ then
318
C q2om2 = pw^2/mtop^2
319
q2om2=single_top_scale(3)/single_top_scale(4)
320
sudakov_exp=sudakov_exp+sudakov_massive(q2merge0
321
$ ,single_top_scale(1),scales(j),q2om2)
322
bornfac=bornfac+ expsudakov_massive(q2merge0
323
$ ,single_top_scale(1),scales(j),q2om2)
324
endif
325
enddo
326
else
106
p=pin
107
lflav=flav
108
nUnClusteredPart=nlegborn
109
110
C --- ANALYSIS
111
C
112
Call ClusterEvent(p,lflav, q2merge,tMom,bMom,WMom,jmom)
113
114
C - End of kinematic analysis and merging, now compute dynamical factors below
115
C --- END OF ANALYSIS
116
117
118
C hard scale, average of the scales in the Sudakovs
119
m2top=abs(tMom(0)**2-tMom(1)**2-tMom(2)**2-tMom(3)**2)
120
Q2W=abs(WMom(0)**2-WMom(1)**2-WMom(2)**2-WMom(3)**2)
121
two_pt_dot_pb=Q2W+m2top
122
q2hard = (Q2W + two_pt_dot_pb)/2d0
123
124
if(q2merge.eq.1d10) then
327
125
c If scales(1)=0 no merge has taken place: no sudakovs.
328
mu2=1
329
muf2=max(q2merge*facfact2,1d0)
330
inlofac=0
126
C muf2=max(q2mergeMAX*facfact2,1d0)
127
muf2=max(q2hard*facfact2,1d0)
331
128
bornfac=0
332
basicfac=1
333
nlofac=0
334
endif
335
if(st_bornorder.gt.inlofac) then
336
alphas=pwhg_alphas(max(q2merge*renfac2,1d0),
129
alphas=pwhg_alphas(max(q2hard*renfac2,1d0),
337
130
1 st_lambda5MSB,st_nlight)
338
c do j=inlofac+1,st_bornorder
339
c mu2=mu2*max(q2merge*renfac2,1d0)
340
c nlofac=nlofac+alphas/st_alpha
341
c basicfac=basicfac*alphas/st_alpha
342
c enddo
343
mu2=mu2*max(q2merge*renfac2,1d0)**(st_bornorder-inlofac)
344
nlofac=nlofac+alphas/st_alpha*(st_bornorder-inlofac)
345
basicfac=basicfac*(alphas/st_alpha)**(st_bornorder-inlofac)
131
mu2=max(q2hard*renfac2,1d0)
132
nlofac=alphas/st_alpha
133
basicfac=(alphas/st_alpha)
346
134
inlofac=st_bornorder
135
136
elseif(q2merge.gt.0d0.and.q2merge.lt.1d10) then
137
bornfac=0
138
alphas=pwhg_alphas(max(q2merge*renfac2,1d0),
139
1 st_lambda5MSB,st_nlight)
140
basicfac=alphas/st_alpha
141
nlofac=basicfac
142
mu2=max(q2merge*renfac2,1d0)
143
muf2=max(q2merge*facfac2,1d0)
144
inlofac=1
145
146
do j=1,nlegborn
147
if(abs(lflav(j)).lt.5) then
148
C light Sudakovs
149
sudakov_exp=sudakov_exp+sudakov(q2merge
150
$ ,Q2W,q2merge,lflav(j))
151
bornfac=bornfac+ expsudakov(q2merge,q2W
152
$ ,q2merge,lflav(j))
153
elseif(abs(lflav(j)).eq.5) then
154
C bottom Sudakov
155
sudakov_exp=sudakov_exp+sudakov(q2merge
156
$ ,two_pt_dot_pb,q2merge,lflav(j))
157
bornfac=bornfac+ expsudakov(q2merge
158
$ ,two_pt_dot_pb,q2merge,lflav(j))
159
elseif(abs(lflav(j)).eq.6) then
160
C top Sudakov
161
q2om2=two_pt_dot_pb/m2top
162
sudakov_exp=sudakov_exp+sudakov_massive(q2merge
163
$ ,two_pt_dot_pb,q2merge,q2om2)
164
bornfac=bornfac+ expsudakov_massive(q2merge
165
$ ,two_pt_dot_pb,q2merge,q2om2)
166
endif
167
enddo
168
else
169
write(*,*) ''
170
write(*,*) 'setlocalscales.f'
171
write(*,*) 'q2merge = ',q2merge
172
write(*,*) 'Value of q2merge not allowed!'
347
173
endif
348
nlofac=nlofac/inlofac
349
mu2=mu2**(1d0/inlofac)
350
174
b0=(33-2*st_nlight)/(12*pi)
351
175
bornfac=1+st_alpha*nlofac*
352
176
1 (bornfac+st_bornorder*b0*log(mu2/st_muren2))
353
177
st_mufact2=muf2
354
178
basicfac=basicfac*exp(sudakov_exp)
355
c$$$ write (*,*) sqrt(q2merge),sqrt(mu2),sqrt(muf2),inlofac
356
c$$$ $ ,st_bornorder,alphas,scales(1),basicfac,bornfac
357
179
end
358
180
359
C $$$C ------------------------------------------------ C
360
C $$$C - Inputs: - C
361
C $$$C - ******* - C
362
C $$$C - p - Underlying born momenta - C
363
C $$$C - lflav - Flavour list derived from - C
364
C $$$C - flst_born by subjecting it to - C
365
C $$$C - repeated QCD clusterings. - C
366
C $$$C - - C
367
C $$$C - Outputs: - C
368
C $$$C - ******** - C
369
C $$$C - jmerge - Index in lflav of one of the two - C
370
C $$$C - closest partons. - C
371
C $$$C - kmerge - Index in lflav of the - C
372
C $$$C - corresponding parton. - C
373
C $$$C - mergedfl - Flavour of parton resulting from - C
374
C $$$C - combination. - C
375
C $$$C - q2merge - pT^2 scale associated to the - C
376
C $$$C - merging of jmerge and kmerge. - C
377
C $$$C - - C
378
C $$$C - checked 24/03/12 - C
379
C $$$C ------------------------------------------------ C
380
C $$$ subroutine findNearestNeighbours(p,lflav,jmerge,kmerge,mergedfl,
381
C $$$ $ q2merge)
382
C $$$ use auxiliary
383
C $$$ implicit none
384
C $$$ include 'nlegborn.h'
385
C $$$ include 'pwhg_st.h'
386
C $$$ include 'pwhg_math.h'
387
C $$$C - Input / output:
388
C $$$ real * 8 p(0:3,nlegborn)
389
C $$$ integer lflav(nlegborn)
390
C $$$ integer jmerge,kmerge,mergedfl
391
C $$$ real * 8 q2merge
392
C $$$C - Local variables:
393
C $$$ real * 8 ycm
394
C $$$ integer onem
395
C $$$ parameter (onem=1000000)
396
C $$$ integer j,k
397
C $$$ integer fl1,fl2,fl
398
C $$$ integer npartons,nparticles
399
C $$$ real * 8 yj,phij,q2j
400
C $$$ real * 8 yk,phik,q2k
401
C $$$ real * 8 dphi
402
C $$$ real * 8 q2
403
C $$$ logical dijetflag
404
C $$$ common/cdijetflag/dijetflag
405
C $$$
406
C $$$ q2merge=1d10
407
C $$$ ycm=log(p(0,1)/p(0,2))/2
408
C $$$ mergedfl=onem
409
C $$$
410
C $$$c Count particles and partons in the final state.
411
C $$$c If we have two particles and two partons, it
412
C $$$c is the dijet case, return with no merging.
413
C $$$
414
C $$$ npartons = 0
415
C $$$ nparticles = 0
416
C $$$ do j=3,nlegborn
417
C $$$ if(abs(lflav(j)).le.st_nlight) npartons = npartons+1
418
C $$$ if(lflav(j).ne.mergedfl) nparticles = nparticles+1
419
C $$$ enddo
420
C $$$ if(npartons.eq.nparticles.and.npartons.eq.2) then
421
C $$$ dijetflag = .true.
422
C $$$ return
423
C $$$ else
424
C $$$ dijetflag = .false.
425
C $$$ endif
426
C $$$c
427
C $$$
428
C $$$ do j=3,nlegborn
429
C $$$ if(abs(lflav(j)).gt.st_nlight) goto 11
430
C $$$ yj=0.5d0*log((p(0,j)+p(3,j))/(p(0,j)-p(3,j)))
431
C $$$ if(yj.gt.ycm) then
432
C $$$ call aux_validmergeisr(lflav,1,j,fl1)
433
C $$$ if(fl1.ne.onem) then
434
C $$$ q2j = p(1,j)**2+p(2,j)**2
435
C $$$ if(q2j.lt.q2merge) then
436
C $$$ q2merge=q2j
437
C $$$ jmerge=1
438
C $$$ kmerge=j
439
C $$$ mergedfl=fl1
440
C $$$ endif
441
C $$$ endif
442
C $$$ else
443
C $$$ call aux_validmergeisr(lflav,2,j,fl2)
444
C $$$ if(fl2.ne.onem) then
445
C $$$ q2j = p(1,j)**2+p(2,j)**2
446
C $$$ if(q2j.lt.q2merge) then
447
C $$$ q2merge=q2j
448
C $$$ jmerge=2
449
C $$$ kmerge=j
450
C $$$ mergedfl=fl2
451
C $$$ endif
452
C $$$ endif
453
C $$$ endif
454
C $$$ do k=j+1,nlegborn
455
C $$$ if(abs(lflav(k)).gt.st_nlight) goto 12
456
C $$$ call aux_validmergefsr(lflav,j,k,fl)
457
C $$$ if(fl.ne.onem) then
458
C $$$ yk=0.5d0*log((p(0,k)+p(3,k))/(p(0,k)-p(3,k)))
459
C $$$ call aux_phipt2(p(:,k),phik,q2k)
460
C $$$ call aux_phipt2(p(:,j),phij,q2j)
461
C $$$ dphi=abs(phik-phij)
462
C $$$ if(dphi.gt.2*pi) dphi=dphi-2*pi
463
C $$$ if(dphi.gt.pi) dphi=2*pi-dphi
464
C $$$ q2=((yk-yj)**2+dphi**2)*min(q2k,q2j)
465
C $$$ if(q2.lt.q2merge) then
466
C $$$ q2merge=q2
467
C $$$ jmerge=j
468
C $$$ kmerge=k
469
C $$$ mergedfl=fl
470
C $$$ endif
471
C $$$ endif
472
C $$$ 12 continue
473
C $$$ enddo
474
C $$$ 11 continue
475
C $$$ enddo
476
C $$$ end
477
478
181
C ------------------------------------------------ C
479
182
C - Inputs: - C
480
183
C - ******* - C
505
208
endif
506
209
lam2=st_lambda5MSB**2
507
210
if(q20.le.lam2.or.q2l.lt.lam2.or.q2h.lt.lam2) then
211
sudakov=-99d99
212
goto 999
213
endif
214
if(q2l.ge.q2h.or.q2h.le.q20) then
508
215
sudakov=0
509
216
goto 999
510
217
endif
511
if(q2l.ge.q2h.or.q2h.le.q20) then
512
sudakov=1
513
goto 999
514
endif
515
218
if(flav.eq.0) then
516
219
isQuark=.false.
517
220
else
537
240
include 'pwhg_st.h'
538
241
include 'pwhg_math.h'
539
242
real * 8 lam2
540
real * 8 theExponentN,theExponentD
243
real * 8 theExponentN,theExponentD,theExponentNQC,theExponentDQC
541
244
lam2=st_lambda5MSB**2
542
245
if(q20.le.lam2.or.q2l.lt.lam2.or.q2h.lt.lam2) then
246
sudakov_massive=-99d99
247
goto 999
248
endif
249
if(q2l.ge.q2h.or.q2h.le.q20) then
543
250
sudakov_massive=0
544
251
goto 999
545
252
endif
546
if(q2l.ge.q2h.or.q2h.le.q20) then
547
sudakov_massive=1
548
goto 999
549
endif
550
253
if(q2l.le.q20) then
551
254
call sudakov_exponent_massive(q20,q2h,q2om2,theExponentN)
552
sudakov_massive=theExponentN
255
call sudakov_exponent_quasi_col(q20,q2h,q2om2,theExponentNQC)
256
sudakov_massive=theExponentN+theExponentNQC
553
257
else
554
258
call sudakov_exponent_massive(q20,q2h,q2om2,theExponentN)
555
259
call sudakov_exponent_massive(q20,q2l,q2om2,theExponentD)
556
sudakov_massive=theExponentN-theExponentD
260
call sudakov_exponent_quasi_col(q20,q2h,q2om2,theExponentNQC)
261
call sudakov_exponent_quasi_col(q20,q2l,q2om2,theExponentDQC)
262
sudakov_massive=theExponentN-theExponentD+theExponentNQC
263
$ -theExponentDQC
557
264
endif
558
265
999 continue
559
266
end
674
381
675
382
function expsudakov_massive(q20,q2h,q2l,q2om2)
676
383
implicit none
677
real * 8 expsudakov_massive,q2h,q2l,q20,q2om2
384
real * 8 expsudakov_massive,q2h,q2l,q20,q2om2,mt2,mt,qh,ql,q0
678
385
include 'pwhg_st.h'
679
386
include 'pwhg_math.h'
680
387
include 'pwhg_flg.h'
681
real * 8 b,lam2
388
real * 8 b,lam2,ddilog
389
external ddilog
682
390
lam2=st_lambda5MSB**2
683
391
if(q20.le.lam2.or.q2l.lt.lam2.or.q2h.lt.lam2) then
684
392
c in this case everything is zero, irrelevant
698
406
1 1d0/pi*( - b*log(q2h/q20))
699
407
2 - 1d0/pi*( - b*log(q2l/q20))
700
408
endif
409
410
c add the quasi-collinear stuff
411
mt2=q2h/q2om2
412
mt=sqrt(mt2)
413
qh=sqrt(q2h)
414
ql=sqrt(q2l)
415
q0=sqrt(q20)
416
if(q2l.le.q20) then
417
expsudakov_massive=expsudakov_massive+4d0/3d0 * 1d0/(2d0*pi)
418
$ *((-4*mt*qh*ATan(mt/qh) + 4*mt*q0*ATan(mt/q0) + q2h
419
$ *Log(1 + mt2/q2h) - mt2*Log(mt2 + q2h) + q20 *Log(q20) +
420
$ (mt2 - q20)*Log(mt2 + q20) + 2*mt2*(ddilog(-(q2h/mt2)) -
421
$ ddilog(-(q20/mt2))))/(2.*mt2))
422
else
423
expsudakov_massive=expsudakov_massive+4d0/3d0 * 1d0/(2d0*pi)
424
$ *((-4*mt*qh*ATan(mt/qh) + 4*mt*q0*ATan(mt/q0) + q2h
425
$ *Log(1 + mt2/q2h) - mt2*Log(mt2 + q2h) + q20 *Log(q20) +
426
$ (mt2 - q20)*Log(mt2 + q20) + 2*mt2*(ddilog(-(q2h/mt2)) -
427
$ ddilog(-(q20/mt2))))/(2.*mt2))
428
expsudakov_massive=expsudakov_massive-4d0/3d0 * 1d0/(2d0*pi)
429
$ *((-4*mt*ql*ATan(mt/ql) + 4*mt*q0*ATan(mt/q0) + q2l
430
$ *Log(1 + mt2/q2l) - mt2*Log(mt2 + q2l) + q20 *Log(q20) +
431
$ (mt2 - q20)*Log(mt2 + q20) + 2*mt2*(ddilog(-(q2l/mt2)) -
432
$ ddilog(-(q20/mt2))))/(2.*mt2))
433
endif
701
434
end
702
435
703
C $$$C ---------------------------------------------- C
704
C $$$C - Inputs: - C
705
C $$$C - ******* - C
706
C $$$C - flav - flavour list derived from flst_born - C
707
C $$$C - by subjecting it to repeated QCD - C
708
C $$$C - compatible clusterings. - C
709
C $$$C - i - index of i-th initial-state partON - C
710
C $$$C - in flav: hence i = 1 or 2 only. - C
711
C $$$C - j - index of j-th final-state partICLE - C
712
C $$$C - particle in flav. - C
713
C $$$C - - C
714
C $$$C - Outputs: - C
715
C $$$C - ******** - C
716
C $$$C - fl - Would-be PDG code of spacelike - C
717
C $$$C - "mother" parton obtained by merging - C
718
C $$$C - (~on-shell) incoming parton i with - C
719
C $$$C - outgoing particle j: - C
720
C $$$C - i -> fl + j - C
721
C $$$C - Note gluons have id=0 in Powheg-Box - C
722
C $$$C - instead of 21. If the splitting is - C
723
C $$$C - not possible in QCD, fl=1000000 ; - C
724
C $$$C - this setting signals to the rest of - C
725
C $$$C - the algorithm that this is not a - C
726
C $$$C - candidate pair for combination. - C
727
C $$$C - - C
728
C $$$C ---------------------------------------------- C
729
C $$$ subroutine validmergeisr(flav,i,j,fl)
730
C $$$ implicit none
731
C $$$ include 'nlegborn.h'
732
C $$$ include 'pwhg_flst.h'
733
C $$$ include 'pwhg_st.h'
734
C $$$ integer onem
735
C $$$ parameter (onem=1000000)
736
C $$$ integer flav(nlegborn),i,j,fl
737
C $$$ integer lflav(nlegborn)
738
C $$$C logical validflav
739
C $$$C external validflav
740
C $$$ interface
741
C $$$ logical function validflav(flav)
742
C $$$ integer flav(:)
743
C $$$ end function
744
C $$$ end interface
745
C $$$ if(i.gt.2.or.j.le.2) then ! Remove when development is finished.
746
C $$$ write(*,*) 'validmergeisr: fatal error'
747
C $$$ write(*,*) 'Routine demands an i.s. and f.s. particle'
748
C $$$ write(*,*) 'index for the 2nd and 3rd input values '
749
C $$$ write(*,*) 'respectively. Quitting.'
750
C $$$ call exit(-1)
751
C $$$ endif
752
C $$$ if(abs(flav(i)).gt.st_nlight.or.abs(flav(j)).gt.st_nlight) then
753
C $$$ fl=onem
754
C $$$ return
755
C $$$ endif
756
C $$$ if(flav(i).eq.flav(j)) then
757
C $$$c g -> g g or q -> g q
758
C $$$ fl=0
759
C $$$ goto 999
760
C $$$ endif
761
C $$$ if(flav(j).eq.0) then
762
C $$$c q -> q g
763
C $$$ fl=flav(i)
764
C $$$ goto 999
765
C $$$ endif
766
C $$$ if(flav(i).eq.0) then
767
C $$$c g -> qbar q
768
C $$$ fl=-flav(j)
769
C $$$ goto 999
770
C $$$ endif
771
C $$$ fl=onem
772
C $$$ return
773
C $$$ 999 continue
774
C $$$C - Check that the flavour list that results from the merging
775
C $$$C - is acceptable e.g. check that for HJJ you don't get back to
776
C $$$C - qqbar->H; if you do then set fl to 1000000, as if the
777
C $$$C - branching were not possible in QCD s.t. it will be neglected
778
C $$$C - as a candidate for clustering.
779
C $$$ lflav=flav
780
C $$$ lflav(j)=onem
781
C $$$ lflav(i)=fl
782
C $$$ if(.not.validflav(lflav)) then
783
C $$$ fl=onem
784
C $$$ endif
785
C $$$ end
786
C $$$
787
C $$$
788
C $$$C ---------------------------------------------- C
789
C $$$C - Inputs: - C
790
C $$$C - ******* - C
791
C $$$C - flav - flavour list derived from flst_born - C
792
C $$$C - by subjecting it to repeated QCD - C
793
C $$$C - compatible clusterings. - C
794
C $$$C - i - index of i-th final-state partICLE - C
795
C $$$C - in flav: hence i = 1 or 2 only. - C
796
C $$$C - j - index of j-th final-state partICLE - C
797
C $$$C - particle in flav. - C
798
C $$$C - - C
799
C $$$C - Outputs: - C
800
C $$$C - ******** - C
801
C $$$C - fl - Would-be PDG code of timelike - C
802
C $$$C - "mother" parton obtained by merging - C
803
C $$$C - outgoing particles i and j: - C
804
C $$$C - fl -> i + j - C
805
C $$$C - Note gluons have id=0 in Powheg-Box - C
806
C $$$C - instead of 21. If the splitting is - C
807
C $$$C - not possible in QCD, fl=1000000 ; - C
808
C $$$C - this setting signals to the rest of - C
809
C $$$C - the algorithm that this is not a - C
810
C $$$C - candidate pair for combination. - C
811
C $$$C - - C
812
C $$$C ---------------------------------------------- C
813
C $$$ subroutine validmergefsr(flav,i,j,fl)
814
C $$$ implicit none
815
C $$$ include 'nlegborn.h'
816
C $$$ include 'pwhg_flst.h'
817
C $$$ include 'pwhg_st.h'
818
C $$$ integer onem
819
C $$$ parameter (onem=1000000)
820
C $$$ integer flav(nlegborn),i,j,fl
821
C $$$ integer lflav(nlegborn)
822
C $$$C logical validflav
823
C $$$C external validflav
824
C $$$ interface
825
C $$$ logical function validflav(flav)
826
C $$$ integer flav(:)
827
C $$$ end function
828
C $$$ end interface
829
C $$$
830
C $$$ if(i.le.2.or.j.le.2) then ! Remove when development is finished.
831
C $$$ write(*,*) 'validmergefsr: fatal error'
832
C $$$ write(*,*) 'Routine demands an f.s. and f.s. particle'
833
C $$$ write(*,*) 'index for the 2nd and 3rd input values '
834
C $$$ write(*,*) 'respectively. Quitting.'
835
C $$$ call exit(-1)
836
C $$$ endif
837
C $$$ if(abs(flav(i)).gt.st_nlight.or.abs(flav(j)).gt.st_nlight) then
838
C $$$ fl=onem
839
C $$$ return
840
C $$$ endif
841
C $$$ if(flav(i).eq.-flav(j)) then
842
C $$$c g -> g g or g -> q qbar
843
C $$$ fl=0
844
C $$$ goto 999
845
C $$$ endif
846
C $$$ if(flav(j).eq.0) then
847
C $$$c q -> q g
848
C $$$ fl=flav(i)
849
C $$$ goto 999
850
C $$$ endif
851
C $$$ if(flav(i).eq.0) then
852
C $$$c q -> g q
853
C $$$ fl=flav(j)
854
C $$$ goto 999
855
C $$$ endif
856
C $$$ fl=onem
857
C $$$ return
858
C $$$ 999 continue
859
C $$$C - Check that the flavour list that results from the merging
860
C $$$C - is acceptable e.g. check that for HJJ you don't get back to
861
C $$$C - qqbar->H; if you do then set fl to 1000000, as if the
862
C $$$C - branching were not possible in QCD s.t. it will be neglected
863
C $$$C - as a candidate for clustering.
864
C $$$ lflav=flav
865
C $$$ lflav(j)=onem
866
C $$$ lflav(i)=fl
867
C $$$ if(.not.validflav(lflav)) then
868
C $$$ fl=onem
869
C $$$ endif
870
C $$$ end
871
C $$$
872
C $$$
873
C $$$
874
C $$$C ---------------------------------------------- C
875
C $$$C - Inputs: - C
876
C $$$C - ******* - C
877
C $$$C - p - p(0) = Energy, p(3) = p_Z - C
878
C $$$C - - C
879
C $$$C - Outputs: - C
880
C $$$C - ******** - C
881
C $$$C - y - Rapidity - C
882
C $$$C - phi - phi - C
883
C $$$C - q2 - pT^2 w.r.t the beam - C
884
C $$$C - - C
885
C $$$C ---------------------------------------------- C
886
C $$$ subroutine phipt2(p,phi,q2)
887
C $$$ implicit none
888
C $$$ real * 8 p(0:3),phi,q2
889
C $$$ q2=p(1)**2+p(2)**2
890
C $$$ phi=atan2(p(2),p(1))
891
C $$$ end
892
893
436
894
437
C ********* DDT / Ellis-Veseli / Nason-Ridolfi Sudakov ************ C
895
438
C - - C
1205
748
aSbar = pwhg_alphas(q2l,st_lambda5MSB,-1)/2/Pi
1206
749
1207
750
B1 = CF*(log(q2om2)-1d0)
1208
K = (67d0/18-Pi**2/6)*CA-5d0/9*nf
1209
B2 = B1*K
1210
751
K = (67d0/18-Pi**2/6)*CA-5d0/9*nf
752
B2 = B1*K
1211
753
1212
754
theB1coeff =
1213
755
$ - Log(Lq2h/Lq2l)/(2.*bnf*Pi)
1215
757
$ + Lq2l*Log(Lq2h)
1216
758
$ )/(2.*bnf**2*Lq2l*Lq2h*Pi)
1217
759
1218
theB2coeff =
760
theB2coeff =
1219
761
$ - (1/Lq2l - 1/Lq2h)/(4.*bnf**2*Pi**2)
1220
762
$ - bpnf*( Lq2l**2 - Lq2h**2 - 2*Lq2h**2*Log(Lq2l)
1221
763
$ + 2*Lq2l**2*Log(Lq2h)
1226
768
$ )/(108.*bnf**4*Lq2l**3*Lq2h**3*Pi**2)
1227
769
1228
770
1229
C theExponent = B1*theB1coeff
1230
theExponent = B1*theB1coeff + B2*theB2coeff
1231
1232
if(.not.pwhg_isfinite(theExponent)) then
1233
write(6,*) ' '
1234
write(6,*) 'Warning: sudakov_exponent is weird.'
1235
write(6,*) 'theExponent = ',theExponent
1236
write(6,*) 'exp(theExponent) = ',exp(theExponent)
1237
write(6,*) 'q_low = ',sqrt(q2l)
1238
write(6,*) 'q_hi = ',sqrt(q2h)
1239
write(6,*) 'q/m = ',sqrt(q2om2)
1240
endif
1241
1242
end
1243
771
C theExponent = B1*theB1coeff
772
theExponent = B1*theB1coeff + B2*theB2coeff
773
774
if(.not.pwhg_isfinite(theExponent)) then
775
write(6,*) ' '
776
write(6,*) 'Warning: sudakov_exponent is weird.'
777
write(6,*) 'theExponent = ',theExponent
778
write(6,*) 'exp(theExponent) = ',exp(theExponent)
779
write(6,*) 'q_low = ',sqrt(q2l)
780
write(6,*) 'q_hi = ',sqrt(q2h)
781
write(6,*) 'q/m = ',sqrt(q2om2)
782
endif
783
784
end
785
786
787
subroutine sudakov_exponent_quasi_col(q2l,q2h,q2om2,theExponent)
788
implicit none
789
real * 8 q2l,q2h,q2om2,theExponent,eps
790
logical pwhg_isfinite
791
external pwhg_isfinite
792
real * 8 m2_common
793
common/sudakov_integral_quasi_col/m2_common
794
real * 8 dgauss,sudakov_exponent_integrand_quasi_col
795
external dgauss,sudakov_exponent_integrand_quasi_col
796
797
eps = 1d-6
798
m2_common=q2h/q2om2
799
theExponent=dgauss(sudakov_exponent_integrand_quasi_col,q2l,q2h,eps)
800
801
if(.not.pwhg_isfinite(theExponent)) then
802
write(6,*) ' '
803
write(6,*) 'Warning: sudakov_exponent is weird.'
804
write(6,*) 'theExponent = ',theExponent
805
write(6,*) 'exp(theExponent) = ',exp(theExponent)
806
write(6,*) 'q_low = ',sqrt(q2l)
807
write(6,*) 'q_hi = ',sqrt(q2h)
808
write(6,*) 'q/m = ',sqrt(q2om2)
809
endif
810
811
end
812
813
double precision function sudakov_exponent_integrand_quasi_col(qt2)
814
include 'nlegborn.h'
815
include 'pwhg_st.h'
816
include 'pwhg_flst.h'
817
include 'pwhg_rad.h'
818
include 'pwhg_math.h'
819
real * 8 qt2,m2
820
real * 8 m2_common
821
common/sudakov_integral_quasi_col/m2_common
822
real * 8 aSbar
823
integer nf
824
real * 8 bnf,bpnf,K
825
real * 8 pwhg_alphas
826
external pwhg_alphas
827
nf=5
828
c The max(qt2,1d0) in the argument of alpha-s is to prevent numerical
829
c inaccuracies. This approximation is completely justified, since this
830
c integrand tents to zero for qt2<<mt2 by construction.
831
aSbar = pwhg_alphas(max(qt2,1d0),st_lambda5MSB,nf)/2/Pi
832
bnf = (11d0*CA-2d0*nf)/12/Pi
833
bpnf = (153 - 19d0*nf) / Pi / 2 / (33 - 2*nf)
834
K = (67d0/18-Pi**2/6)*CA-5d0/9*nf
835
m2=m2_common
836
sudakov_exponent_integrand_quasi_col=
837
& -Cf*(aSbar+K*aSbar**2)*(
838
& log(m2/qt2)-
839
& sqrt(qt2/m2)*atan(sqrt(m2/qt2))-
840
& (1d0-qt2/(2*m2))*log(1+m2/qt2)
841
& )/qt2
842
return
843
end
844
1244
845
1245
846
C ***************************************************************** C
1246
847
C - The integrand in Sudakov exponent times q^2 i.e. we effectiv - C
Loggerhead is a web-based interface for Breezy
Version: 2.0.1