72
74
real * 8 pin(0:3,nlegborn),basicfac,bornfac,nlofac
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),
77
79
integer j,k,l,lflav(nlegborn),jmerge,kmerge,inlofac,npart
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)
86
common/cdijetflag/dijetflag
87
logical raisingscales,ini
88
save raisingscales,ini
86
real * 8 q2mergeMAX,q2hard
90
87
logical pwhg_isfinite
91
integer n_badq2merge,ixx
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)
96
if(powheginput("#raisingscales").eq.0) then
97
raisingscales = .false.
99
raisingscales = .true.
94
C - Initialization of dynamical parameters
103
95
renfac2=st_renfact**2
104
96
facfact2=st_facfact**2
100
C - Make sure momenta and flavour arrays have no residual junk in them
104
C - Initialize kinematics
111
call findNearestNeighbours(p,lflav,jmerge,kmerge,mergedfl,
113
if(.not.pwhg_isfinite(q2merge)) then
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 = '
128
write(*,*) ixx,p(:,ixx)
131
if(n_badq2merge.eq.100) then
132
write(*,*) 'Further bad q2 merge messages are suppressed'
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
147
c save this scale; it is the Q_0 scale that appears in all Sudakovs
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
155
mu2=max(q2merge*renfac2,1d0)
156
muf2=max(q2merge*facfact2,1d0)
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))
165
1 expsudakov(q2merge0,q2merge,lscalej,lflav(jmerge))
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
177
p(:,jmerge)=p(:,jmerge)+p(:,kmerge)
179
p(3,jmerge)=p(3,jmerge)-p(3,kmerge)
180
p(0,jmerge)=abs(p(3,jmerge))
183
lflav(jmerge)=mergedfl
190
c No more merging is possible.
192
if(.not.dijetflag) then
193
c Define the initial scale as
194
c the invariant mass of the remaining system
197
if(lflav(j).ne.onem) then
201
if(raisingscales) then
202
q2merge=max(q2mergeMAX,
203
$ ptot(0)**2-ptot(1)**2-ptot(2)**2-ptot(3)**2)
205
q2merge=ptot(0)**2-ptot(1)**2-ptot(2)**2-ptot(3)**2
208
c SINGLE TOP SPECIAL: first compute -Q^2 (W-boson t-channel)
211
if (lflav(j).ne.onem) npart=npart+1
214
if (abs(lflav(1)).eq.5) then
222
if (abs(lflav(j)).eq.6) then
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
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
243
if (abs(lflav(j)).eq.5 .or. abs(lflav(j)).eq.6) then
249
c if b is initial state, gluon is attached to massless quark line
250
if (lflav(1).ne.0) then
258
if (abs(lflav(j)).eq.6) then
265
write (*,*) 'ERROR #1 in setting cluster scales',npart
273
if (abs(lflav(j)).eq.6) then
275
& p(0,j)**2-p(1,j)**2-p(2,j)**2-p(3,j)**2
278
! typical hard scale for massive quark line:
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:
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
291
c Dijet case: use the scalar sum of the pt of the two partons
294
if(lflav(j).ne.onem) then
295
q2merge=sqrt(p(1,j)**2+p(2,j)**2)+q2merge
299
if(raisingscales) then
300
q2merge=max(q2mergeMAX,q2merge)
304
if(scales(1).gt.0) then
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)
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)
108
nUnClusteredPart=nlegborn
112
Call ClusterEvent(p,lflav, q2merge,tMom,bMom,WMom,jmom)
114
C - End of kinematic analysis and merging, now compute dynamical factors below
115
C --- END OF ANALYSIS
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
124
if(q2merge.eq.1d10) then
327
125
c If scales(1)=0 no merge has taken place: no sudakovs.
329
muf2=max(q2merge*facfact2,1d0)
126
C muf2=max(q2mergeMAX*facfact2,1d0)
127
muf2=max(q2hard*facfact2,1d0)
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
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
136
elseif(q2merge.gt.0d0.and.q2merge.lt.1d10) then
138
alphas=pwhg_alphas(max(q2merge*renfac2,1d0),
139
1 st_lambda5MSB,st_nlight)
140
basicfac=alphas/st_alpha
142
mu2=max(q2merge*renfac2,1d0)
143
muf2=max(q2merge*facfac2,1d0)
147
if(abs(lflav(j)).lt.5) then
149
sudakov_exp=sudakov_exp+sudakov(q2merge
150
$ ,Q2W,q2merge,lflav(j))
151
bornfac=bornfac+ expsudakov(q2merge,q2W
153
elseif(abs(lflav(j)).eq.5) then
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
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)
170
write(*,*) 'setlocalscales.f'
171
write(*,*) 'q2merge = ',q2merge
172
write(*,*) 'Value of q2merge not allowed!'
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))
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
359
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
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
378
C $$$C - checked 24/03/12 - C
379
C $$$C ------------------------------------------------ C
380
C $$$ subroutine findNearestNeighbours(p,lflav,jmerge,kmerge,mergedfl,
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:
395
C $$$ parameter (onem=1000000)
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
403
C $$$ logical dijetflag
404
C $$$ common/cdijetflag/dijetflag
407
C $$$ ycm=log(p(0,1)/p(0,2))/2
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.
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
420
C $$$ if(npartons.eq.nparticles.and.npartons.eq.2) then
421
C $$$ dijetflag = .true.
424
C $$$ dijetflag = .false.
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
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
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
478
181
C ------------------------------------------------ C
698
406
1 1d0/pi*( - b*log(q2h/q20))
699
407
2 - 1d0/pi*( - b*log(q2l/q20))
410
c add the quasi-collinear stuff
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))
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))
703
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
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
728
C $$$C ---------------------------------------------- C
729
C $$$ subroutine validmergeisr(flav,i,j,fl)
731
C $$$ include 'nlegborn.h'
732
C $$$ include 'pwhg_flst.h'
733
C $$$ include 'pwhg_st.h'
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
741
C $$$ logical function validflav(flav)
742
C $$$ integer flav(:)
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.'
752
C $$$ if(abs(flav(i)).gt.st_nlight.or.abs(flav(j)).gt.st_nlight) then
756
C $$$ if(flav(i).eq.flav(j)) then
757
C $$$c g -> g g or q -> g q
761
C $$$ if(flav(j).eq.0) then
766
C $$$ if(flav(i).eq.0) then
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.
782
C $$$ if(.not.validflav(lflav)) then
788
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
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
812
C $$$C ---------------------------------------------- C
813
C $$$ subroutine validmergefsr(flav,i,j,fl)
815
C $$$ include 'nlegborn.h'
816
C $$$ include 'pwhg_flst.h'
817
C $$$ include 'pwhg_st.h'
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
825
C $$$ logical function validflav(flav)
826
C $$$ integer flav(:)
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.'
837
C $$$ if(abs(flav(i)).gt.st_nlight.or.abs(flav(j)).gt.st_nlight) then
841
C $$$ if(flav(i).eq.-flav(j)) then
842
C $$$c g -> g g or g -> q qbar
846
C $$$ if(flav(j).eq.0) then
851
C $$$ if(flav(i).eq.0) then
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.
867
C $$$ if(.not.validflav(lflav)) then
874
C $$$C ---------------------------------------------- C
877
C $$$C - p - p(0) = Energy, p(3) = p_Z - 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
885
C $$$C ---------------------------------------------- C
886
C $$$ subroutine phipt2(p,phi,q2)
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))
894
437
C ********* DDT / Ellis-Veseli / Nason-Ridolfi Sudakov ************ C
1226
768
$ )/(108.*bnf**4*Lq2l**3*Lq2h**3*Pi**2)
1229
C theExponent = B1*theB1coeff
1230
theExponent = B1*theB1coeff + B2*theB2coeff
1232
if(.not.pwhg_isfinite(theExponent)) then
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)
771
C theExponent = B1*theB1coeff
772
theExponent = B1*theB1coeff + B2*theB2coeff
774
if(.not.pwhg_isfinite(theExponent)) then
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)
787
subroutine sudakov_exponent_quasi_col(q2l,q2h,q2om2,theExponent)
789
real * 8 q2l,q2h,q2om2,theExponent,eps
790
logical pwhg_isfinite
791
external pwhg_isfinite
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
799
theExponent=dgauss(sudakov_exponent_integrand_quasi_col,q2l,q2h,eps)
801
if(.not.pwhg_isfinite(theExponent)) then
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)
813
double precision function sudakov_exponent_integrand_quasi_col(qt2)
816
include 'pwhg_flst.h'
818
include 'pwhg_math.h'
821
common/sudakov_integral_quasi_col/m2_common
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
836
sudakov_exponent_integrand_quasi_col=
837
& -Cf*(aSbar+K*aSbar**2)*(
839
& sqrt(qt2/m2)*atan(sqrt(m2/qt2))-
840
& (1d0-qt2/(2*m2))*log(1+m2/qt2)
1245
846
C ***************************************************************** C
1246
847
C - The integrand in Sudakov exponent times q^2 i.e. we effectiv - C