1
double precision function gamma(q0)
2
c**************************************************
3
c calculates the branching probability
4
c**************************************************
6
include 'nexternal.inc'
12
double precision q0, val, add, add2
13
double precision qr,lf
14
double precision alphas
17
parameter (pi=3.141592654d0)
21
if (Q1<m_qmass(iipdg)) return
22
m_lastas=Alphas(alpsfact*q0)
23
val=2d0*m_colfac(iipdg)*m_lastas/PI/q0
24
c if (m_mode & bpm::power_corrs) then
26
if(m_pca(iipdg,iimode).eq.0)then
31
val=val*(m_dlog(iipdg)*(1d0+m_kfac*m_lastas/(2d0*PI))*lf+m_slog(iipdg)
32
$ +qr*(m_power(iipdg,1,iimode)+qr*(m_power(iipdg,2,iimode)
33
$ +qr*m_power(iipdg,3,iimode))))
35
c val=val*m_dlog*(1d0+m_kfac*m_lastas/(2d0*PI))*log(Q1/q0)+m_slog;
37
if(m_qmass(iipdg).gt.0d0)then
38
val=val+m_colfac(iipdg)*m_lastas/PI/q0*(0.5-q0/m_qmass(iipdg)*
39
$ atan(m_qmass(iipdg)/q0)-
40
$ (1.0-0.5*(q0/m_qmass(iipdg))**2)*log(1.0+(m_qmass(iipdg)/q0)**2))
46
if(m_qmass(abs(i)).gt.0d0)then
47
add2=m_colfac(i)*m_lastas/PI/q0/
48
$ (1.0+(m_qmass(abs(i))/q0)**2)*
49
$ (1.0-1.0/3.0/(1.0+(m_qmass(abs(i))/q0)**2))
51
add2=2d0*m_colfac(i)*m_lastas/PI/q0*(m_slog(i)
52
$ +qr*(m_power(i,1,iimode)+qr*(m_power(i,2,iimode)
53
$ +qr*m_power(i,3,iimode))))
62
if (btest(mlevel,6)) then
63
write(*,*)' \\Delta^I_{',iipdg,'}(',
64
& q0,',',q1,') -> ',gamma
65
write(*,*) val,m_lastas,m_dlog(iipdg),m_slog(iipdg)
66
write(*,*) m_power(iipdg,1,iimode),m_power(iipdg,2,iimode),m_power(iipdg,3,iimode)
72
double precision function sud(q0,Q11,ipdg,imode)
73
c**************************************************
74
c actually calculates is sudakov weight
75
c**************************************************
78
include 'nexternal.inc'
81
double precision q0, Q11
82
double precision gamma,DGAUSS
93
sud=exp(-DGAUSS(gamma,q0,Q1,eps))
95
if (btest(mlevel,6)) then
96
write(*,*)' \\Delta^',imode,'_{',ipdg,'}(',
97
& 2*log10(q0/q1),') -> ',sud
103
double precision function sudwgt(q0,q1,q2,ipdg,imode)
104
c**************************************************
105
c calculates is sudakov weight
106
c**************************************************
108
include 'message.inc'
110
double precision q0, q1, q2
117
if(q2.lt.q1.and.btest(mlevel,4))
118
$ write(*,*)'Warning! q2 < q1 in sudwgt. Return 1.'
122
sudwgt=sud(q0,q2,ipdg,imode)/sud(q0,q1,ipdg,imode)
124
if (btest(mlevel,5)) then
125
write(*,*)' \\Delta^',imode,'_{',ipdg,'}(',
126
& q0,',',q1,',',q2,') -> ',sudwgt
132
logical function isqcd(ipdg)
133
c**************************************************
134
c determines whether particle is qcd particle
135
c**************************************************
141
irfl=iand(abs(ipdg),63)
142
if (irfl.gt.8.and.irfl.ne.21) isqcd=.false.
143
c write(*,*)'iqcd? pdg = ',ipdg,' -> ',irfl,' -> ',isqcd
148
logical function isjet(ipdg)
149
c**************************************************
150
c determines whether particle is qcd jet particle
151
c**************************************************
161
if (irfl.gt.maxjetflavor.and.irfl.ne.21) isjet=.false.
162
c write(*,*)'isjet? pdg = ',ipdg,' -> ',irfl,' -> ',isjet
167
logical function isparton(ipdg)
168
c**************************************************
169
c determines whether particle is qcd jet particle
170
c**************************************************
180
if (irfl.gt.5.and.irfl.ne.21) isparton=.false.
181
c write(*,*)'isparton? pdg = ',ipdg,' -> ',irfl,' -> ',isparton
187
subroutine ipartupdate(p,imo,ida1,ida2,ipdg,ipart)
188
c**************************************************
189
c Traces particle lines according to CKKW rules
190
c**************************************************
194
include 'nexternal.inc'
195
include 'message.inc'
197
double precision p(0:3,nexternal)
198
integer imo,ida1,ida2,i,idmo,idda1,idda2
199
integer ipdg(n_max_cl),ipart(2,n_max_cl)
209
if (btest(mlevel,1)) then
210
write(*,*) ' updating ipart for: ',ida1,ida2,imo
213
if (btest(mlevel,1)) then
214
write(*,*) ' daughters: ',(ipart(i,ida1),i=1,2),(ipart(i,ida2),i=1,2)
217
c IS clustering - just transmit info on incoming line
218
if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
219
$ (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
220
if(ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2)
221
$ ipart(1,imo)=ipart(1,ida2)
222
if(ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2)
223
$ ipart(1,imo)=ipart(1,ida1)
224
if (btest(mlevel,1)) then
225
write(*,*) ' -> ',(ipart(i,imo),i=1,2)
226
c Set intermediate particle identity
227
if(iabs(idmo).lt.6)then
228
if(iabs(idda1).lt.6) ipdg(imo)=-idda1
229
if(iabs(idda2).lt.6) ipdg(imo)=-idda2
231
if (btest(mlevel,1)) then
232
write(*,*) ' particle identities: ',idda1,idda2,idmo
240
if(idmo.eq.21.and.idda1.eq.21.and.idda2.eq.21)then
241
c gluon -> 2 gluon splitting: Choose hardest gluon
242
if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
243
$ p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
244
ipart(1,imo)=ipart(1,ida1)
245
ipart(2,imo)=ipart(2,ida1)
247
ipart(1,imo)=ipart(1,ida2)
248
ipart(2,imo)=ipart(2,ida2)
250
else if(idmo.eq.21.and.idda1.eq.-idda2)then
251
c gluon -> quark anti-quark: use both, but take hardest as 1
252
if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
253
$ p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
254
ipart(1,imo)=ipart(1,ida1)
255
ipart(2,imo)=ipart(1,ida2)
257
ipart(1,imo)=ipart(1,ida2)
258
ipart(2,imo)=ipart(1,ida1)
260
else if(idmo.eq.idda1.or.idmo.eq.idda1+sign(1,idda2))then
261
c quark -> quark-gluon or quark-Z or quark-h or quark-W
262
ipart(1,imo)=ipart(1,ida1)
263
else if(idmo.eq.idda2.or.idmo.eq.idda2+sign(1,idda1))then
264
c quark -> gluon-quark or Z-quark or h-quark or W-quark
265
ipart(1,imo)=ipart(1,ida2)
268
if (btest(mlevel,1)) then
269
write(*,*) ' -> ',(ipart(i,imo),i=1,2)
272
c Set intermediate particle identity
273
if(iabs(idmo).lt.6)then
274
if(iabs(idda1).lt.6) ipdg(imo)=idda1
275
if(iabs(idda2).lt.6) ipdg(imo)=idda2
277
if (btest(mlevel,1)) then
278
write(*,*) ' particle identities: ',idda1,idda2,idmo
285
logical function isjetvx(imo,ida1,ida2,ipdg,ipart)
286
c***************************************************
287
c Checks if a qcd vertex generates a jet
288
c***************************************************
292
include 'nexternal.inc'
294
integer imo,ida1,ida2,idmo,idda1,idda2,i
295
integer ipdg(n_max_cl),ipart(2,n_max_cl)
304
if(.not.isqcd(idmo).or..not.isqcd(idda1).or.
305
& .not.isqcd(idda2)) then
311
if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
312
$ (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
313
c Check if ida1 is outgoing parton or ida2 is outgoing parton
314
if(ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2.and.isjet(idda1).or.
315
$ ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2.and.isjet(idda2))then
324
if(isjet(idda1).or.isjet(idda2))then
333
logical function setclscales(p)
334
c**************************************************
335
c reweight the hard me according to ckkw
336
c employing the information in common/cl_val/
337
c**************************************************
340
include 'message.inc'
342
include 'nexternal.inc'
343
include 'cluster.inc'
349
DOUBLE PRECISION P(0:3,NEXTERNAL)
352
integer i, j, idi, idj
354
parameter( PI = 3.14159265358979323846d0 )
356
integer mapconfig(0:lmaxconfigs), this_config
357
integer iforest(2,-max_branch:-1,lmaxconfigs)
358
integer sprop(-max_branch:-1,lmaxconfigs)
359
integer tprid(-max_branch:-1,lmaxconfigs)
360
include 'configs.inc'
361
real*8 xptj,xptb,xpta,xptl,xmtc
362
real*8 xetamin,xqcut,deltaeta
363
common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
366
include 'maxamps.inc'
367
integer idup(nexternal,maxproc)
368
integer mothup(2,nexternal,maxproc)
369
integer icolup(2,nexternal,maxflow)
370
include 'leshouche.inc'
371
double precision asref, pt2prev(n_max_cl),pt2min
372
integer n, icmp, ibeam(2), iqcd(0:2)!, ilast(0:nexternal)
373
integer idfl, ipdg(n_max_cl), idmap(-nexternal:nexternal)
374
integer ipart(2,n_max_cl)
375
double precision xnow(2)
376
integer jlast(2),jfirst(2)
377
logical qcdline(2),qcdrad(2)
381
logical isqcd,isjet,isparton,isjetvx,cluster
382
double precision alphas
383
external isqcd, isjet, isparton, isjetvx, cluster, alphas
387
if(ickkw.le.0.and.xqcut.le.0d0.and.q2fact(1).gt.0.and.scale.gt.0) return
390
c Cluster the configuration
393
if (.not.cluster(p(0,1))) then
394
c if (xqcut.gt.0d0) then
396
cc if(pt2ijcl(1).lt.xqcut**2) failed=.true.
398
c if (btest(mlevel,3)) then
399
c write(*,*)'q_min = ',pt2ijcl(1),' < ',xqcut**2
401
c setclscales=.false.
406
write(*,*)'setclscales: Error. Clustering failed.'
411
c Preparing graph particle information
413
ipart(1,ishft(1,i))=i
414
ipart(2,ishft(1,i))=0
418
if (btest(mlevel,1)) then
419
write(*,*)'setclscales: identified tree {'
421
write(*,*)' ',i,': ',idacl(i,1),'&',idacl(i,2),
422
& ' -> ',imocl(i),', ptij = ',dsqrt(pt2ijcl(i))
424
write(*,*)' graphs (',igscl(0),'):',(igscl(i),i=1,igscl(0))
427
c fill particle information
428
icmp=ishft(1,nexternal+1)-2
431
ipdg(idmap(i))=idup(i,1)
433
$ write(*,*) i,' got id ',idmap(i),' -> ',ipdg(idmap(i))
436
idi=iforest(1,-i,igscl(1))
437
idj=iforest(2,-i,igscl(1))
438
idmap(-i)=idmap(idi)+idmap(idj)
439
idfl=sprop(-i,igscl(1))
442
ipdg(icmp-idmap(-i))=idfl
444
idfl=tprid(-i,igscl(1))
447
ipdg(icmp-idmap(-i))=idfl
449
c write(*,*) -i,' (',idi,',',idj,') got id ',idmap(-i),
450
c & ' -> ',ipdg(idmap(-i))
454
cc Set factorization scale as for the MLM case
456
c if(xqcut.gt.0) then
457
cc Using last clustering value
458
c if(pt2ijcl(nexternal-2).lt.max(4d0,xqcut**2))then
459
c setclscales=.false.
463
c If last clustering is s-channel QCD (e.g. ttbar) use mt2last instead
464
c (i.e. geom. average of transverse mass of t and t~)
465
if(mt2last.gt.4d0 .and. nexternal.gt.3 .and. isqcd(ipdg(idacl(nexternal-3,1)))
466
$ .and. isqcd(ipdg(idacl(nexternal-3,2)))
467
$ .and. isqcd(ipdg(imocl(nexternal-3))))then
468
mt2ij(nexternal-2)=mt2last
469
mt2ij(nexternal-3)=mt2last
470
if (btest(mlevel,3)) then
471
write(*,*)' setclscales: set last vertices to mtlast: ',sqrt(mt2last)
475
C If we have fixed factorization scale, for ickkw>0 means central
476
C scale, i.e. last two scales (ren. scale for these vertices are
477
C anyway already set by "scale" above
479
if(fixed_fac_scale.and.first)then
483
else if(fixed_fac_scale) then
498
c Go through clusterings and set factorization scales for use in dsig
500
c Update particle tree map
501
call ipartupdate(p,imocl(n),idacl(n,1),idacl(n,2),
505
if (isqcd(ipdg(idacl(n,i)))) then
507
if ((isparton(ipdg(idacl(n,i))).and.idacl(n,i).eq.ibeam(j).or.
508
$ isparton(ipdg(imocl(n))).and.imocl(n).eq.ibeam(j))
509
$ .and.qcdrad(j)) then
510
c is emission - this is what we want
511
c Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
512
c f1(x1,pt2E) is given by DSIG, just need to set scale.
514
if(ickkw.eq.0.or.jfirst(j).eq.0) jfirst(j)=n
516
qcdline(j)=isqcd(ipdg(imocl(n)))
517
if(n.lt.nexternal-2)then
518
qcdrad(j)=isqcd(ipdg(idacl(n,3-i)))
527
$ write(*,*) 'jfirst is ',jfirst(1),jfirst(2),
528
$ ' and jlast is ',jlast(1),jlast(2)
530
c Set central scale to mT2 and multiply with scalefact
531
if(jlast(1).gt.0.and.mt2ij(jlast(1)).gt.0d0)
532
$ pt2ijcl(jlast(1))=mt2ij(jlast(1))
533
if(jlast(2).gt.0.and.mt2ij(jlast(2)).gt.0d0)
534
$ pt2ijcl(jlast(2))=mt2ij(jlast(2))
535
if(qcdline(1).and.qcdline(2).and.jlast(1).ne.jlast(2)) then
536
c If not WBF or similar, set uniform scale to be geom. average
537
pt2ijcl(jlast(1))=sqrt(pt2ijcl(jlast(1))*pt2ijcl(jlast(2)))
538
pt2ijcl(jlast(2))=pt2ijcl(jlast(1))
540
if(jlast(1).gt.0) pt2ijcl(jlast(1))=scalefact**2*pt2ijcl(jlast(1))
541
if(jlast(2).gt.0) pt2ijcl(jlast(2))=scalefact**2*pt2ijcl(jlast(2))
542
if(lpp(1).eq.0.and.lpp(2).eq.0)then
543
pt2ijcl(nexternal-2)=scalefact**2*pt2ijcl(nexternal-2)
544
pt2ijcl(nexternal-3)=pt2ijcl(nexternal-2)
547
c Check xqcut for vertices with jet daughters only
550
if (n.lt.nexternal-2.and.n.ne.jlast(1).and.n.ne.jlast(2).and.
551
$ (isjet(ipdg(idacl(n,1))).or.isjet(ipdg(idacl(n,2)))).and.
552
$ sqrt(pt2ijcl(n)).lt.xqcut)then
559
c JA: Check xmtc cut for central process
560
if(pt2ijcl(jlast(1)).lt.xmtc**2.or.pt2ijcl(jlast(2)).lt.xmtc**2)then
562
if(btest(mlevel,3)) write(*,*)'Failed xmtc cut ',
563
$ sqrt(pt2ijcl(jlast(1))),sqrt(pt2ijcl(jlast(1))),' < ',xmtc
567
if(ickkw.eq.0.and.(fixed_fac_scale.or.q2fact(1).gt.0).and.
568
$ (fixed_ren_scale.or.scale.gt.0)) return
570
c Set renormalization scale to largest factorization scale
571
if(scale.eq.0d0) then
572
if(jlast(1).gt.0.and.jlast(2).gt.0) then
573
scale=(pt2ijcl(jlast(1))*pt2ijcl(jlast(2)))**0.25d0
574
elseif(jlast(1).gt.0) then
575
scale=sqrt(pt2ijcl(jlast(1)))
576
elseif(jlast(2).gt.0) then
577
scale=sqrt(pt2ijcl(jlast(2)))
578
elseif(lpp(1).eq.0.and.lpp(2).eq.0) then
579
scale=sqrt(pt2ijcl(nexternal-2))
582
$ G = SQRT(4d0*PI*ALPHAS(scale))
585
$ write(*,*) 'Set ren scale to ',scale
587
if(ickkw.gt.0.and.q2fact(1).gt.0) then
588
c Use the fixed or previously set scale for central scale
589
if(jlast(1).gt.0) pt2ijcl(jlast(1))=q2fact(1)
590
if(jlast(2).gt.0) pt2ijcl(jlast(2))=q2fact(2)
593
if(lpp(1).eq.0.and.lpp(2).eq.0)then
594
if(q2fact(1).gt.0)then
595
pt2ijcl(nexternal-2)=q2fact(1)
596
pt2ijcl(nexternal-3)=q2fact(1)
598
q2fact(1)=pt2ijcl(nexternal-2)
601
elseif(ickkw.eq.2.or.pdfwgt)then
602
c Use the minimum scale found for fact scale in ME
603
if(jlast(1).gt.0) q2fact(1)=min(pt2ijcl(jfirst(1)),pt2ijcl(jlast(1)))
604
if(jlast(2).gt.0) q2fact(2)=min(pt2ijcl(jfirst(2)),pt2ijcl(jlast(2)))
605
else if(q2fact(1).eq.0d0) then
606
if(jlast(1).gt.0) q2fact(1)=pt2ijcl(jlast(1))
607
if(jlast(2).gt.0) q2fact(2)=pt2ijcl(jlast(2))
610
if(nexternal.eq.3.and.nincoming.eq.2) then
612
$ q2fact(1)=pt2ijcl(nexternal-2)
614
$ q2fact(2)=pt2ijcl(nexternal-2)
618
$ write(*,*) 'Set fact scales to ',sqrt(q2fact(1)),sqrt(q2fact(2))
622
double precision function rewgt(p)
623
c**************************************************
624
c reweight the hard me according to ckkw
625
c employing the information in common/cl_val/
626
c**************************************************
629
include 'message.inc'
631
include 'nexternal.inc'
632
include 'cluster.inc'
638
DOUBLE PRECISION P(0:3,NEXTERNAL)
642
DOUBLE PRECISION PD(0:MAXPROC)
643
COMMON /SubProc/ PD, IPROC
646
integer i, j, idi, idj
648
parameter( PI = 3.14159265358979323846d0 )
650
integer mapconfig(0:lmaxconfigs), this_config
651
integer iforest(2,-max_branch:-1,lmaxconfigs)
652
integer sprop(-max_branch:-1,lmaxconfigs)
653
integer tprid(-max_branch:-1,lmaxconfigs)
654
include 'configs.inc'
655
real*8 xptj,xptb,xpta,xptl,xmtc
656
real*8 xetamin,xqcut,deltaeta
657
common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
658
include 'maxamps.inc'
659
integer idup(nexternal,maxproc)
660
integer mothup(2,nexternal,maxproc)
661
integer icolup(2,nexternal,maxflow)
662
include 'leshouche.inc'
663
double precision asref, pt2prev(n_max_cl),pt2pdf(n_max_cl),pt2min
664
integer n, icmp, ibeam(2), iqcd(0:2)!, ilast(0:nexternal)
665
integer idfl, ipdg(n_max_cl), idmap(-nexternal:nexternal)
666
integer ipart(2,n_max_cl)
667
double precision xnow(2)
668
double precision xtarget, tmp
673
logical isqcd,isjet,isparton,isjetvx
674
double precision alphas,getissud,pdg2pdf, sudwgt
676
external isqcd,isjet,isparton
677
external alphas, isjetvx, getissud, pdg2pdf, xran1, sudwgt
681
if(ickkw.le.0) return
683
if(.not.clustered)then
684
write(*,*)'Error: No clustering done when calling rewgt!'
689
c Set mimimum kt scale, depending on highest mult or not
690
if(hmult.or.ickkw.eq.1)then
696
$ write(*,*) 'pt2min set to ',pt2min
698
c Since we use pdf reweighting, need to know particle identities
701
xtarget=xran1(iseed)*pd(np)
703
do while (pd(iprocset) .lt. xtarget .and. iprocset .lt. np)
706
if (btest(mlevel,1)) then
707
write(*,*) 'Set process number ',iprocset
710
c Preparing graph particle information
712
c ilast(i)=ishft(1,i)
714
pt2prev(ishft(1,i))=max(pt2min,p(0,i)**2-p(1,i)**2-p(2,i)**2-p(3,i)**2)
716
pt2prev(ishft(1,i))=0d0
718
pt2pdf(ishft(1,i))=pt2prev(ishft(1,i))
719
ptclus(i)=sqrt(pt2prev(ishft(1,i)))
720
ipart(1,ishft(1,i))=i
721
ipart(2,ishft(1,i))=0
726
if (btest(mlevel,1)) then
727
write(*,*)'rewgt: identified tree {'
729
write(*,*)' ',i,': ',idacl(i,1),'&',idacl(i,2),
730
& ' -> ',imocl(i),', ptij = ',dsqrt(pt2ijcl(i))
732
write(*,*)' graphs (',igscl(0),'):',(igscl(i),i=1,igscl(0))
735
c fill particle information
736
icmp=ishft(1,nexternal+1)-2
739
ipdg(idmap(i))=idup(i,iprocset)
741
$ write(*,*) i,' got id ',idmap(i),' -> ',ipdg(idmap(i))
744
idi=iforest(1,-i,igscl(1))
745
idj=iforest(2,-i,igscl(1))
746
idmap(-i)=idmap(idi)+idmap(idj)
747
idfl=sprop(-i,igscl(1))
750
ipdg(icmp-idmap(-i))=idfl
752
idfl=tprid(-i,igscl(1))
755
ipdg(icmp-idmap(-i))=idfl
757
c write(*,*) -i,' (',idi,',',idj,') got id ',idmap(-i),
758
c & ' -> ',ipdg(idmap(-i))
760
c Set x values for the two sides, for IS Sudakovs
764
if(btest(mlevel,3))then
765
write(*,*) 'Set x values to ',xnow(1),xnow(2)
768
c Prepare for resetting q2fact based on PDF reweighting
769
if(ickkw.eq.2.or.pdfwgt)then
774
c Set strong coupling used
778
c Perform alpha_s reweighting based on type of vertex
780
if (btest(mlevel,3)) then
781
write(*,*)' ',n,': ',idacl(n,1),'(',ipdg(idacl(n,1)),
782
& ')&',idacl(n,2),'(',ipdg(idacl(n,2)),') -> ',
783
& imocl(n),'(',ipdg(imocl(n)),'), ptij = ',
786
c perform alpha_s reweighting only for vertices where a jet is produced
787
c and not for the last clustering (use non-fixed ren. scale for these)
788
if (n.lt.nexternal-2.and.
789
$ isjetvx(imocl(n),idacl(n,1),idacl(n,2),ipdg,ipart)) then
791
rewgt=rewgt*alphas(alpsfact*sqrt(pt2ijcl(n)))/asref
792
if (btest(mlevel,3)) then
793
write(*,*)' reweight vertex: ',ipdg(imocl(n)),ipdg(idacl(n,1)),ipdg(idacl(n,2))
794
write(*,*)' as: ',alphas(alpsfact*dsqrt(pt2ijcl(n))),
795
& '/',asref,' -> ',alphas(alpsfact*dsqrt(pt2ijcl(n)))/asref
796
write(*,*)' and G=',SQRT(4d0*PI*ALPHAS(scale))
799
c Update starting values for FS parton showering
802
if(ipart(j,idacl(n,i)).gt.0)then
803
ptclus(ipart(j,idacl(n,i)))=dsqrt(pt2ijcl(n))
807
c Update particle tree map
808
call ipartupdate(p,imocl(n),idacl(n,1),idacl(n,2),ipdg,ipart)
809
if(ickkw.eq.2.or.pdfwgt) then
810
c Perform PDF and, if ickkw=2, Sudakov reweighting
813
c write(*,*)'weight ',idacl(n,i),', ptij=',pt2prev(idacl(n,i))
814
if (isqcd(ipdg(idacl(n,i)))) then
815
if(pt2min.eq.0d0) then
818
$ write(*,*) 'pt2min set to ',pt2min
820
if(pt2prev(idacl(n,i)).eq.0d0) pt2prev(idacl(n,i))=
821
$ max(pt2min,p(0,i)**2-p(1,i)**2-p(2,i)**2-p(3,i)**2)
823
if (isparton(ipdg(idacl(n,i))).and.idacl(n,i).eq.ibeam(j)) then
824
c is sudakov weight - calculate only once for each parton
825
c line where parton line ends with change of parton id or
826
c non-radiation vertex
829
if(pt2pdf(idacl(n,i)).eq.0d0) pt2pdf(idacl(n,i))=pt2prev(idacl(n,i))
830
if(ickkw.eq.2.and.(ipdg(idacl(n,i)).ne.ipdg(imocl(n)).or.
831
$ .not.isjetvx(imocl(n),idacl(n,1),idacl(n,2),ipdg,ipart)).and.
832
$ pt2prev(idacl(n,i)).lt.pt2ijcl(n).and.zcl(n).gt.1d-20)then
833
tmp=min(1d0,max(getissud(ibeam(j),ipdg(idacl(n,i)),xnow(j),xnow(3-j),pt2ijcl(n)),1d-20)/
834
$ max(getissud(ibeam(j),ipdg(idacl(n,i)),xnow(j),xnow(3-j),pt2prev(idacl(n,i))),1d-20))
836
pt2prev(imocl(n))=pt2ijcl(n)
837
if (btest(mlevel,3)) then
838
write(*,*)' reweight line: ',ipdg(idacl(n,i)), idacl(n,i)
839
write(*,*)' pt2prev, pt2new, x1, x2: ',pt2prev(idacl(n,i)),pt2ijcl(n),xnow(j),xnow(3-j)
840
write(*,*)' Sud: ',tmp
841
write(*,*)' -> rewgt: ',rewgt
844
pt2prev(imocl(n))=pt2prev(idacl(n,i))
846
c Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
847
c f1(x1,pt2E) is given by DSIG, already set scale for that
848
xnow(j)=xnow(j)*zcl(n)
849
if(q2fact(j).eq.0d0.and.ickkw.eq.2)then
850
q2fact(j)=pt2min ! Starting scale for PS
851
if (btest(mlevel,3)) then
852
write(*,*)' reweight: set fact scale ',j,' for PS scale to: ',q2fact(j)
854
else if(q2fact(j).eq.0d0)then
856
else if(pt2pdf(idacl(n,i)).lt.pt2ijcl(n).and.zcl(n).gt.1d-20) then
857
if(ickkw.eq.1) q2fact(j)=pt2ijcl(n)
858
rewgt=rewgt*max(pdg2pdf(abs(ibeam(j)),ipdg(idacl(n,i))*sign(1,ibeam(j)),xnow(j),sqrt(pt2ijcl(n))),1d-20)/
859
$ max(pdg2pdf(abs(ibeam(j)),ipdg(idacl(n,i))*sign(1,ibeam(j)),xnow(j),sqrt(pt2pdf(idacl(n,i)))),1d-20)
860
if (btest(mlevel,3)) then
861
write(*,*)' reweight ',ipdg(idacl(n,i)),' by pdfs: '
862
write(*,*)' x, pt2prev, ptnew: ',xnow(j),pt2pdf(idacl(n,i)),pt2ijcl(n)
864
$ pdg2pdf(abs(ibeam(j)),ipdg(idacl(n,i))*sign(1,ibeam(j)),xnow(j),sqrt(pt2ijcl(n))),' / ',
865
$ pdg2pdf(abs(ibeam(j)),ipdg(idacl(n,i))*sign(1,ibeam(j)),xnow(j),sqrt(pt2pdf(idacl(n,i))))
866
write(*,*)' -> rewgt: ',rewgt
867
c write(*,*)' (compare for glue: ',
868
c $ pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i)))),' / ',
869
c $ pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2ijcl(n)))
870
c write(*,*)' = ',pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i))))/
871
c $ pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2ijcl(n)))
872
c write(*,*)' -> ',pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i))))/
873
c $ pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2ijcl(n)))*rewgt,' )'
876
c End both Sudakov and pdf reweighting when we reach a
877
c non-radiation vertex
878
if(isjetvx(imocl(n),idacl(n,1),idacl(n,2),ipdg,ipart)) then
879
pt2pdf(imocl(n))=pt2ijcl(n)
881
pt2pdf(imocl(n))=1d30
882
pt2prev(imocl(n))=1d30
883
if (btest(mlevel,3)) then
884
write(*,*)' rewgt: for vertex ',idacl(n,1),idacl(n,2),imocl(n),
885
$ ' with ids ',ipdg(idacl(n,1)),ipdg(idacl(n,2)),ipdg(imocl(n))
886
write(*,*)' set pt2prev, pt2pdf: ',pt2prev(imocl(n)),pt2pdf(imocl(n))
893
if(ickkw.eq.2.and.pt2prev(idacl(n,i)).lt.pt2ijcl(n).and.
894
$ (isvx.or.ipdg(idacl(n,i)).ne.ipdg(imocl(n)).or.
895
$ (ipdg(idacl(n,i)).ne.ipdg(idacl(n,3-i)).and.
896
$ pt2prev(idacl(n,i)).gt.pt2prev(idacl(n,3-i))))) then
897
tmp=sudwgt(sqrt(pt2min),sqrt(pt2prev(idacl(n,i))),
898
& dsqrt(pt2ijcl(n)),ipdg(idacl(n,i)),1)
900
if (btest(mlevel,3)) then
901
write(*,*)' reweight fs line: ',ipdg(idacl(n,i)), idacl(n,i)
902
write(*,*)' pt2prev, pt2new: ',pt2prev(idacl(n,i)),pt2ijcl(n)
903
write(*,*)' Sud: ',tmp
904
write(*,*)' -> rewgt: ',rewgt
906
pt2prev(imocl(n))=pt2ijcl(n)
908
pt2prev(imocl(n))=pt2prev(idacl(n,i))
913
if (ickkw.eq.2.and.n.eq.nexternal-2.and.isqcd(ipdg(imocl(n))).and.
914
$ pt2prev(imocl(n)).lt.pt2ijcl(n)) then
915
tmp=sudwgt(sqrt(pt2min),sqrt(pt2prev(imocl(n))),
916
& dsqrt(pt2ijcl(n)),ipdg(imocl(n)),1)
918
if (btest(mlevel,3)) then
919
write(*,*)' reweight last fs line: ',ipdg(imocl(n)), imocl(n)
920
write(*,*)' pt2prev, pt2new: ',pt2prev(imocl(n)),pt2ijcl(n)
921
write(*,*)' Sud: ',tmp
922
write(*,*)' -> rewgt: ',rewgt
928
if((ickkw.eq.2.or.pdfwgt).and.lpp(1).eq.0.and.lpp(2).eq.0)then
933
if (btest(mlevel,3)) then
934
write(*,*)'} -> w = ',rewgt
added
removed
Template/SubProcesses/reweight.f
