← Back to branch summary

~maddevelopers/mg5amcnlo/WWW_tg

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_427003/PROC_427003/SubProcesses/reweight.f

  • Committer: Martin Delcourt
  • Date: 2013-03-25 20:56:51 UTC
  • Revision ID: mdel@perso.be-20130325205651-f3bisj4z2jlfyjh3
added test et test2

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_427003/PROC_427003/SubProcesses/reweight.f
1
 
      double precision function gamma(q0)
2
 
c**************************************************
3
 
c   calculates the branching probability
4
 
c**************************************************
5
 
      implicit none
6
 
      include 'nexternal.inc'
7
 
      include 'message.inc'
8
 
      include 'maxamps.inc'
9
 
      include 'cluster.inc'
10
 
      include 'sudakov.inc'
11
 
      include 'run.inc'
12
 
      integer i
13
 
      double precision q0, val, add, add2
14
 
      double precision qr,lf
15
 
      double precision alphas
16
 
      external alphas
17
 
      double precision pi
18
 
      parameter (pi=3.141592654d0)
19
 
 
20
 
      gamma=0.0d0
21
 
 
22
 
      if (Q1<m_qmass(iipdg)) return
23
 
      m_lastas=Alphas(alpsfact*q0)
24
 
      val=2d0*m_colfac(iipdg)*m_lastas/PI/q0
25
 
c   if (m_mode & bpm::power_corrs) then
26
 
      qr=q0/Q1
27
 
      if(m_pca(iipdg,iimode).eq.0)then
28
 
        lf=log(1d0/qr-1d0)
29
 
      else 
30
 
        lf=log(1d0/qr)
31
 
      endif
32
 
      val=val*(m_dlog(iipdg)*(1d0+m_kfac*m_lastas/(2d0*PI))*lf+m_slog(iipdg)
33
 
     $   +qr*(m_power(iipdg,1,iimode)+qr*(m_power(iipdg,2,iimode)
34
 
     $   +qr*m_power(iipdg,3,iimode))))
35
 
c   else
36
 
c   val=val*m_dlog*(1d0+m_kfac*m_lastas/(2d0*PI))*log(Q1/q0)+m_slog;
37
 
c   endif
38
 
      if(m_qmass(iipdg).gt.0d0)then
39
 
        val=val+m_colfac(iipdg)*m_lastas/PI/q0*(0.5-q0/m_qmass(iipdg)*
40
 
     $     atan(m_qmass(iipdg)/q0)-
41
 
     $     (1.0-0.5*(q0/m_qmass(iipdg))**2)*log(1.0+(m_qmass(iipdg)/q0)**2))
42
 
      endif
43
 
      val=max(val,0d0)
44
 
      if (iipdg.eq.21) then
45
 
        add=0d0
46
 
        do i=-6,-1
47
 
          if(m_qmass(abs(i)).gt.0d0)then
48
 
            add2=m_colfac(i)*m_lastas/PI/q0/
49
 
     $         (1.0+(m_qmass(abs(i))/q0)**2)*
50
 
     $         (1.0-1.0/3.0/(1.0+(m_qmass(abs(i))/q0)**2))
51
 
          else
52
 
            add2=2d0*m_colfac(i)*m_lastas/PI/q0*(m_slog(i)
53
 
     $         +qr*(m_power(i,1,iimode)+qr*(m_power(i,2,iimode)
54
 
     $         +qr*m_power(i,3,iimode))))
55
 
          endif
56
 
          add=add+max(add2,0d0)
57
 
        enddo
58
 
        val=val+add
59
 
      endif
60
 
      
61
 
      gamma = max(val,0d0)
62
 
 
63
 
      if (btest(mlevel,6)) then
64
 
        write(*,*)'       \\Delta^I_{',iipdg,'}(',
65
 
     &     q0,',',q1,') -> ',gamma
66
 
        write(*,*) val,m_lastas,m_dlog(iipdg),m_slog(iipdg)
67
 
        write(*,*) m_power(iipdg,1,iimode),m_power(iipdg,2,iimode),m_power(iipdg,3,iimode)
68
 
      endif
69
 
 
70
 
      return
71
 
      end
72
 
 
73
 
      double precision function sud(q0,Q11,ipdg,imode)
74
 
c**************************************************
75
 
c   actually calculates is sudakov weight
76
 
c**************************************************
77
 
      implicit none
78
 
      include 'message.inc'
79
 
      include 'nexternal.inc'
80
 
      include 'maxamps.inc'
81
 
      include 'cluster.inc'      
82
 
      integer ipdg,imode
83
 
      double precision q0, Q11
84
 
      double precision gamma,DGAUSS
85
 
      external gamma,DGAUSS
86
 
      double precision eps
87
 
      parameter (eps=1d-5)
88
 
      
89
 
      sud=0.0d0
90
 
 
91
 
      Q1=Q11
92
 
      iipdg=iabs(ipdg)
93
 
      iimode=imode
94
 
 
95
 
      sud=exp(-DGAUSS(gamma,q0,Q1,eps))
96
 
 
97
 
      if (btest(mlevel,6)) then
98
 
        write(*,*)'       \\Delta^',imode,'_{',ipdg,'}(',
99
 
     &     2*log10(q0/q1),') -> ',sud
100
 
      endif
101
 
 
102
 
      return
103
 
      end
104
 
 
105
 
      double precision function sudwgt(q0,q1,q2,ipdg,imode)
106
 
c**************************************************
107
 
c   calculates is sudakov weight
108
 
c**************************************************
109
 
      implicit none
110
 
      include 'message.inc'
111
 
      integer ipdg,imode
112
 
      double precision q0, q1, q2
113
 
      double precision sud
114
 
      external sud
115
 
      
116
 
      sudwgt=1.0d0
117
 
 
118
 
      if(q2.le.q1)then
119
 
         if(q2.lt.q1.and.btest(mlevel,4))
120
 
     $        write(*,*)'Warning! q2 < q1 in sudwgt. Return 1.'
121
 
         return
122
 
      endif
123
 
 
124
 
      sudwgt=sud(q0,q2,ipdg,imode)/sud(q0,q1,ipdg,imode)
125
 
 
126
 
      if (btest(mlevel,5)) then
127
 
        write(*,*)'       \\Delta^',imode,'_{',ipdg,'}(',
128
 
     &     q0,',',q1,',',q2,') -> ',sudwgt
129
 
      endif
130
 
 
131
 
      return
132
 
      end
133
 
 
134
 
      logical function isqcd(ipdg)
135
 
c**************************************************
136
 
c   determines whether particle is qcd particle
137
 
c**************************************************
138
 
      implicit none
139
 
      integer ipdg, irfl
140
 
      integer get_color
141
 
 
142
 
      isqcd=(iabs(get_color(ipdg)).gt.1)
143
 
 
144
 
      return
145
 
      end
146
 
 
147
 
      logical function isjet(ipdg)
148
 
c**************************************************
149
 
c   determines whether particle is qcd jet particle
150
 
c**************************************************
151
 
      implicit none
152
 
 
153
 
      include 'cuts.inc'
154
 
 
155
 
      integer ipdg, irfl
156
 
 
157
 
      isjet=.true.
158
 
 
159
 
      irfl=abs(ipdg)
160
 
      if (irfl.gt.maxjetflavor.and.irfl.ne.21) isjet=.false.
161
 
c      write(*,*)'isjet? pdg = ',ipdg,' -> ',irfl,' -> ',isjet
162
 
 
163
 
      return
164
 
      end
165
 
 
166
 
      logical function isparton(ipdg)
167
 
c**************************************************
168
 
c   determines whether particle is qcd jet particle
169
 
c**************************************************
170
 
      implicit none
171
 
 
172
 
      include 'cuts.inc'
173
 
      include 'run.inc'
174
 
 
175
 
      integer ipdg, irfl
176
 
 
177
 
      isparton=.true.
178
 
 
179
 
      irfl=abs(ipdg)
180
 
      if (irfl.gt.max(asrwgtflavor,maxjetflavor).and.irfl.ne.21)
181
 
     $     isparton=.false.
182
 
c      write(*,*)'isparton? pdg = ',ipdg,' -> ',irfl,' -> ',isparton
183
 
 
184
 
      return
185
 
      end
186
 
 
187
 
 
188
 
      subroutine ipartupdate(p,imo,ida1,ida2,ipdg,ipart)
189
 
c**************************************************
190
 
c   Traces particle lines according to CKKW rules
191
 
c**************************************************
192
 
      implicit none
193
 
 
194
 
      include 'ncombs.inc'
195
 
      include 'nexternal.inc'
196
 
      include 'message.inc'
197
 
 
198
 
      double precision p(0:3,nexternal)
199
 
      integer imo,ida1,ida2,i,idmo,idda1,idda2
200
 
      integer ipdg(n_max_cl),ipart(2,n_max_cl)
201
 
 
202
 
      idmo=ipdg(imo)
203
 
      idda1=ipdg(ida1)
204
 
      idda2=ipdg(ida2)
205
 
 
206
 
      if (btest(mlevel,4)) then
207
 
        write(*,*) ' updating ipart for: ',ida1,ida2,' -> ',imo
208
 
      endif
209
 
 
210
 
      if (btest(mlevel,4)) then
211
 
         write(*,*) ' daughters: ',(ipart(i,ida1),i=1,2),(ipart(i,ida2),i=1,2)
212
 
      endif
213
 
 
214
 
c     IS clustering - just transmit info on incoming line
215
 
      if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
216
 
     $   (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
217
 
         ipart(2,imo)=0
218
 
         if(ipart(1,ida1).le.2.and.ipart(1,ida2).le.2)then
219
 
c           This is last clustering - keep mother ipart
220
 
            ipart(1,imo)=ipart(1,imo)
221
 
         elseif(ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2)then
222
 
            ipart(1,imo)=ipart(1,ida2)        
223
 
         elseif(ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2)then
224
 
            ipart(1,imo)=ipart(1,ida1)
225
 
         endif
226
 
         if (btest(mlevel,4))
227
 
     $        write(*,*) ' -> ',(ipart(i,imo),i=1,2)
228
 
         return
229
 
      endif        
230
 
 
231
 
c     FS clustering
232
 
      if(idmo.eq.21.and.idda1.eq.21.and.idda2.eq.21)then
233
 
c     gluon -> 2 gluon splitting: Choose hardest gluon
234
 
        if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
235
 
     $     p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
236
 
          ipart(1,imo)=ipart(1,ida1)
237
 
          ipart(2,imo)=ipart(2,ida1)
238
 
        else
239
 
          ipart(1,imo)=ipart(1,ida2)
240
 
          ipart(2,imo)=ipart(2,ida2)
241
 
        endif
242
 
      else if(idmo.eq.21.and.idda1.eq.-idda2)then
243
 
c     gluon -> quark anti-quark: use both, but take hardest as 1
244
 
        if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
245
 
     $     p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
246
 
          ipart(1,imo)=ipart(1,ida1)
247
 
          ipart(2,imo)=ipart(1,ida2)
248
 
        else
249
 
          ipart(1,imo)=ipart(1,ida2)
250
 
          ipart(2,imo)=ipart(1,ida1)
251
 
        endif
252
 
      else if(idmo.eq.idda1.or.idmo.eq.idda1+sign(1,idda2))then
253
 
c     quark -> quark-gluon or quark-Z or quark-h or quark-W
254
 
        ipart(1,imo)=ipart(1,ida1)
255
 
        ipart(2,imo)=0
256
 
      else if(idmo.eq.idda2.or.idmo.eq.idda2+sign(1,idda1))then
257
 
c     quark -> gluon-quark or Z-quark or h-quark or W-quark
258
 
        ipart(1,imo)=ipart(1,ida2)
259
 
        ipart(2,imo)=0
260
 
      else
261
 
c     Color singlet
262
 
         ipart(1,imo)=ipart(1,ida1)
263
 
         ipart(2,imo)=ipart(1,ida2)
264
 
      endif
265
 
      
266
 
      if (btest(mlevel,4)) then
267
 
        write(*,*) ' -> ',(ipart(i,imo),i=1,2)
268
 
      endif
269
 
 
270
 
      return
271
 
      end
272
 
      
273
 
      logical function isjetvx(imo,ida1,ida2,ipdg,ipart,islast)
274
 
c***************************************************
275
 
c   Checks if a qcd vertex generates a jet
276
 
c***************************************************
277
 
      implicit none
278
 
 
279
 
      include 'ncombs.inc'
280
 
      include 'nexternal.inc'
281
 
 
282
 
      integer imo,ida1,ida2,idmo,idda1,idda2,i
283
 
      integer ipdg(n_max_cl),ipart(2,n_max_cl)
284
 
      logical isqcd,isjet,islast
285
 
      external isqcd,isjet
286
 
 
287
 
      idmo=ipdg(imo)
288
 
      idda1=ipdg(ida1)
289
 
      idda2=ipdg(ida2)
290
 
 
291
 
c     Check QCD vertex
292
 
      if(islast.or..not.isqcd(idmo).or..not.isqcd(idda1).or.
293
 
     &     .not.isqcd(idda2)) then
294
 
         isjetvx = .false.
295
 
         return
296
 
      endif
297
 
 
298
 
c     IS clustering
299
 
      if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
300
 
     $   (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
301
 
c     Check if ida1 is outgoing parton or ida2 is outgoing parton
302
 
         if(ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2.and.isjet(idda1).or.
303
 
     $        ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2.and.isjet(idda2))then
304
 
           isjetvx=.true.
305
 
        else
306
 
           isjetvx=.false.
307
 
        endif
308
 
        return
309
 
      endif        
310
 
 
311
 
c     FS clustering
312
 
      if(isjet(idda1).or.isjet(idda2))then
313
 
         isjetvx=.true.
314
 
      else
315
 
         isjetvx=.false.
316
 
      endif
317
 
      
318
 
      return
319
 
      end
320
 
 
321
 
      logical function ispartonvx(imo,ida1,ida2,ipdg,ipart,islast)
322
 
c***************************************************
323
 
c   Checks if a qcd vertex generates a jet
324
 
c***************************************************
325
 
      implicit none
326
 
 
327
 
      include 'ncombs.inc'
328
 
      include 'nexternal.inc'
329
 
 
330
 
      integer imo,ida1,ida2,idmo,idda1,idda2,i
331
 
      integer ipdg(n_max_cl),ipart(2,n_max_cl)
332
 
      logical isqcd,isparton,islast
333
 
      external isqcd,isparton
334
 
 
335
 
      idmo=ipdg(imo)
336
 
      idda1=ipdg(ida1)
337
 
      idda2=ipdg(ida2)
338
 
 
339
 
c     Check QCD vertex
340
 
      if(.not.isqcd(idmo).or..not.isqcd(idda1).or.
341
 
     &     .not.isqcd(idda2)) then
342
 
         ispartonvx = .false.
343
 
         return
344
 
      endif
345
 
 
346
 
c     IS clustering
347
 
      if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
348
 
     $   (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
349
 
c     Check if ida1 is outgoing parton or ida2 is outgoing parton
350
 
         if(.not.islast.and.ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2.and.isparton(idda1).or.
351
 
     $        ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2.and.isparton(idda2))then
352
 
           ispartonvx=.true.
353
 
        else
354
 
           ispartonvx=.false.
355
 
        endif
356
 
        return
357
 
      endif        
358
 
 
359
 
c     FS clustering
360
 
      if(isparton(idda1).or.isparton(idda2))then
361
 
         ispartonvx=.true.
362
 
      else
363
 
         ispartonvx=.false.
364
 
      endif
365
 
      
366
 
      return
367
 
      end
368
 
 
369
 
      logical function setclscales(p)
370
 
c**************************************************
371
 
c     Calculate dynamic scales based on clustering
372
 
c     Also perform xqcut and xmtc cuts
373
 
c**************************************************
374
 
      implicit none
375
 
 
376
 
      include 'message.inc'
377
 
      include 'genps.inc'
378
 
      include 'maxconfigs.inc'
379
 
      include 'nexternal.inc'
380
 
      include 'maxamps.inc'
381
 
      include 'cluster.inc'
382
 
      include 'run.inc'
383
 
      include 'coupl.inc'
384
 
C   
385
 
C   ARGUMENTS 
386
 
C   
387
 
      DOUBLE PRECISION P(0:3,NEXTERNAL)
388
 
 
389
 
C   global variables
390
 
C     Present process number
391
 
      INTEGER IMIRROR,IPROC
392
 
      COMMON/TO_MIRROR/IMIRROR, IPROC
393
 
 
394
 
C   local variables
395
 
      integer i, j, idi, idj
396
 
      real*8 PI
397
 
      parameter( PI = 3.14159265358979323846d0 )
398
 
 
399
 
      integer mapconfig(0:lmaxconfigs), this_config
400
 
      integer iforest(2,-max_branch:-1,lmaxconfigs)
401
 
      real*8 xptj,xptb,xpta,xptl,xmtc
402
 
      real*8 xetamin,xqcut,deltaeta
403
 
      common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
404
 
c     q2bck holds the central q2fact scales
405
 
      real*8 q2bck(2)
406
 
      common /to_q2bck/q2bck
407
 
      double precision asref, pt2prev(n_max_cl),pt2min
408
 
      integer n, ibeam(2), iqcd(0:2)!, ilast(0:nexternal)
409
 
      integer idfl, idmap(-nexternal:nexternal)
410
 
      integer ipart(2,n_max_cl)
411
 
      double precision xnow(2)
412
 
      integer jlast(2),jfirst(2),jcentral(2),nwarning
413
 
      logical qcdline(2),partonline(2)
414
 
      logical failed,first
415
 
      data first/.true./
416
 
      data nwarning/0/
417
 
 
418
 
      logical isqcd,isjet,isparton,cluster,isjetvx
419
 
      double precision alphas
420
 
      external isqcd, isjet, isparton, cluster, isjetvx, alphas
421
 
 
422
 
      setclscales=.true.
423
 
 
424
 
      if(ickkw.le.0.and.xqcut.le.0d0.and.q2fact(1).gt.0.and.scale.gt.0) return
425
 
 
426
 
c   
427
 
c   Cluster the configuration
428
 
c   
429
 
      
430
 
c      if (.not.clustered) then
431
 
         clustered = cluster(p(0,1))
432
 
         if(.not.clustered) then
433
 
            open(unit=26,file='../../../error',status='unknown',err=999)
434
 
            write(26,*) 'Error: Clustering failed in cluster.f.'
435
 
            write(*,*) 'Error: Clustering failed in cluster.f.'
436
 
            stop
437
 
 999        write(*,*) 'error'
438
 
            setclscales=.false.
439
 
            clustered = .false.
440
 
            return
441
 
         endif
442
 
c      endif
443
 
 
444
 
      if (btest(mlevel,1)) then
445
 
        write(*,*)'setclscales: identified tree {'
446
 
        do i=1,nexternal-2
447
 
          write(*,*)'  ',i,': ',idacl(i,1),'(',ipdgcl(idacl(i,1),igraphs(1),iproc),')',
448
 
     $       '&',idacl(i,2),'(',ipdgcl(idacl(i,2),igraphs(1),iproc),')',
449
 
     $       ' -> ',imocl(i),'(',ipdgcl(imocl(i),igraphs(1),iproc),')',
450
 
     $       ', ptij = ',dsqrt(pt2ijcl(i))
451
 
        enddo
452
 
        write(*,*)'  process: ',iproc
453
 
        write(*,*)'  graphs (',igraphs(0),'):',(igraphs(i),i=1,igraphs(0))
454
 
        write(*,*)'}'
455
 
      endif
456
 
 
457
 
c     If last clustering is s-channel QCD (e.g. ttbar) use mt2last instead
458
 
c     (i.e. geom. average of transverse mass of t and t~)
459
 
        if(mt2last.gt.4d0 .and. nexternal.gt.3) then
460
 
           if(isqcd(ipdgcl(idacl(nexternal-3,1),igraphs(1),iproc))
461
 
     $      .and. isqcd(ipdgcl(idacl(nexternal-3,2),igraphs(1),iproc))
462
 
     $      .and. isqcd(ipdgcl(imocl(nexternal-3),igraphs(1),iproc)))then
463
 
              mt2ij(nexternal-2)=mt2last
464
 
              mt2ij(nexternal-3)=mt2last
465
 
              if (btest(mlevel,3)) then
466
 
                 write(*,*)' setclscales: set last vertices to mtlast: ',sqrt(mt2last)
467
 
              endif
468
 
           endif
469
 
        endif
470
 
 
471
 
C   If we have fixed factorization scale, for ickkw>0 means central
472
 
C   scale, i.e. last two scales (ren. scale for these vertices are
473
 
C   anyway already set by "scale" above)
474
 
      if(ickkw.gt.0) then
475
 
         if(fixed_fac_scale.and.first)then
476
 
            q2bck(1)=q2fact(1)
477
 
            q2bck(2)=q2fact(2)
478
 
            first=.false.
479
 
         else if(fixed_fac_scale) then
480
 
            q2fact(1)=q2bck(1)
481
 
            q2fact(2)=q2bck(2)
482
 
         endif
483
 
      endif
484
 
 
485
 
      do i=1,2
486
 
         ibeam(i)=ishft(1,i-1)
487
 
         jfirst(i)=0
488
 
         jlast(i)=0
489
 
         jcentral(i)=0
490
 
         partonline(i)=isparton(ipdgcl(ibeam(i),igraphs(1),iproc))
491
 
         qcdline(i)=isqcd(ipdgcl(ibeam(i),igraphs(1),iproc))
492
 
      enddo
493
 
 
494
 
c   Go through clusterings and set factorization scales for use in dsig
495
 
      if (nexternal.eq.3) goto 10
496
 
      do n=1,nexternal-2
497
 
        do i=1,2
498
 
            do j=1,2
499
 
              if (isqcd(ipdgcl(idacl(n,i),igraphs(1),iproc)).and.
500
 
     $            idacl(n,i).eq.ibeam(j).and.qcdline(j)) then
501
 
c             is emission - this is what we want
502
 
c             Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
503
 
c             f1(x1,pt2E) is given by DSIG, just need to set scale.
504
 
                 ibeam(j)=imocl(n)
505
 
                 if(jfirst(j).eq.0)then
506
 
                    if(isjetvx(imocl(n),idacl(n,1),idacl(n,2),
507
 
     $                 ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)) then
508
 
                       jfirst(j)=n
509
 
                    else
510
 
                       jfirst(j)=-1
511
 
                    endif
512
 
                 endif
513
 
                 if(partonline(j))then
514
 
c                   Stop fact scale where parton line stops
515
 
                    jlast(j)=n
516
 
                    partonline(j)=isjet(ipdgcl(imocl(n),igraphs(1),iproc))
517
 
                 endif
518
 
c                Trace QCD line through event
519
 
                 jcentral(j)=n
520
 
                 qcdline(j)=isqcd(ipdgcl(imocl(n),igraphs(1),iproc))
521
 
              endif
522
 
            enddo
523
 
        enddo
524
 
      enddo
525
 
 
526
 
 10   if(jfirst(1).le.0) jfirst(1)=jlast(1)
527
 
      if(jfirst(2).le.0) jfirst(2)=jlast(2)
528
 
 
529
 
      if (btest(mlevel,3))
530
 
     $     write(*,*) 'jfirst is ',jfirst(1),jfirst(2),
531
 
     $     ' jlast is ',jlast(1),jlast(2),
532
 
     $     ' and jcentral is ',jcentral(1),jcentral(2)
533
 
 
534
 
c     Set central scale to mT2
535
 
      if(jcentral(1).gt.0) then
536
 
         if(mt2ij(jcentral(1)).gt.0d0)
537
 
     $        pt2ijcl(jcentral(1))=mt2ij(jcentral(1))
538
 
      endif
539
 
      if(jcentral(2).gt.0)then
540
 
         if(mt2ij(jcentral(2)).gt.0d0)
541
 
     $     pt2ijcl(jcentral(2))=mt2ij(jcentral(2))
542
 
      endif
543
 
      if(btest(mlevel,4))then
544
 
         print *,'jlast, jcentral: ',(jlast(i),i=1,2),(jcentral(i),i=1,2)
545
 
         if(jlast(1).gt.0) write(*,*)'pt(jlast 1): ', sqrt(pt2ijcl(jlast(1)))
546
 
         if(jlast(2).gt.0) write(*,*)'pt(jlast 2): ', sqrt(pt2ijcl(jlast(2)))
547
 
         if(jcentral(1).gt.0) write(*,*)'pt(jcentral 1): ', sqrt(pt2ijcl(jcentral(1)))
548
 
         if(jcentral(2).gt.0) write(*,*)'pt(jcentral 2): ', sqrt(pt2ijcl(jcentral(2)))
549
 
      endif
550
 
c     Check xqcut for vertices with jet daughters only
551
 
      ibeam(1)=ishft(1,0)
552
 
      ibeam(2)=ishft(1,1)
553
 
      if(xqcut.gt.0) then
554
 
         do n=1,nexternal-2
555
 
            do j=1,2
556
 
               if (btest(mlevel,4))
557
 
     $              write(*,*) 'ibeam(1,2), daughter, mother:',
558
 
     $              ibeam(1),ibeam(2),idacl(n,j),imocl(n)
559
 
               if (n.lt.nexternal-2) then
560
 
                  if(n.ne.jlast(1).and.n.ne.jlast(2).and.
561
 
     $              isjet(ipdgcl(idacl(n,j),igraphs(1),iproc)).and.
562
 
     $              (idacl(n,3-j).eq.ibeam(1).or.
563
 
     $              idacl(n,3-j).eq.ibeam(2)).and.
564
 
     $              sqrt(pt2ijcl(n)).lt.xqcut)then
565
 
c                   ISR
566
 
                     if (btest(mlevel,3))
567
 
     $                    write(*,*) 'Failed xqcut: ',n, ipdgcl(idacl(n,1),igraphs(1),iproc),
568
 
     $                    ipdgcl(idacl(n,2),igraphs(1),iproc), xqcut
569
 
                     setclscales=.false.
570
 
                     clustered = .false.
571
 
                     return
572
 
                  endif
573
 
               endif
574
 
               if (idacl(n,j).eq.ibeam(1).and.imocl(n).ne.ibeam(2))
575
 
     $              ibeam(1)=imocl(n)
576
 
               if (idacl(n,j).eq.ibeam(2).and.imocl(n).ne.ibeam(1))
577
 
     $              ibeam(2)=imocl(n)
578
 
            enddo
579
 
            if (n.lt.nexternal-2) then
580
 
               if(n.ne.jlast(1).and.n.ne.jlast(2).and.
581
 
     $           isjet(ipdgcl(idacl(n,1),igraphs(1),iproc)).and.
582
 
     $           isjet(ipdgcl(idacl(n,2),igraphs(1),iproc)).and.
583
 
     $           idacl(n,1).ne.ibeam(1).and.idacl(n,1).ne.ibeam(2).and.
584
 
     $           idacl(n,2).ne.ibeam(1).and.idacl(n,2).ne.ibeam(2).and.
585
 
     $           sqrt(pt2ijcl(n)).lt.xqcut)then
586
 
c           FSR
587
 
                  if (btest(mlevel,3))
588
 
     $                 write(*,*) 'Failed xqcut: ',n, ipdgcl(idacl(n,1),igraphs(1),iproc),
589
 
     $                 ipdgcl(idacl(n,2),igraphs(1),iproc), xqcut
590
 
                  setclscales=.false.
591
 
                  clustered = .false.
592
 
                  return
593
 
               endif
594
 
            endif
595
 
         enddo
596
 
      endif
597
 
 
598
 
c     JA: Check xmtc cut for central process
599
 
      if(xmtc**2.gt.0) then
600
 
         if(jcentral(1).gt.0.and.pt2ijcl(jcentral(1)).lt.xmtc**2
601
 
     $      .or.jcentral(2).gt.0.and.pt2ijcl(jcentral(2)).lt.xmtc**2)then
602
 
            setclscales=.false.
603
 
            clustered = .false.
604
 
            if(btest(mlevel,3)) write(*,*)'Failed xmtc cut ',
605
 
     $           sqrt(pt2ijcl(jcentral(1))),sqrt(pt2ijcl(jcentral(1))),
606
 
     $           ' < ',xmtc
607
 
            return
608
 
         endif
609
 
      endif
610
 
      
611
 
      if(ickkw.eq.0.and.(fixed_fac_scale.or.q2fact(1).gt.0).and.
612
 
     $     (fixed_ren_scale.or.scale.gt.0)) return
613
 
 
614
 
c     Ensure that last scales are at least as big as first scales
615
 
      if(jlast(1).gt.0)
616
 
     $     pt2ijcl(jlast(1))=max(pt2ijcl(jlast(1)),pt2ijcl(jfirst(1)))
617
 
      if(jlast(2).gt.0)
618
 
     $     pt2ijcl(jlast(2))=max(pt2ijcl(jlast(2)),pt2ijcl(jfirst(2)))
619
 
 
620
 
c     Set renormalization scale to geom. aver. of central scales
621
 
c     if both beams are qcd
622
 
      if(scale.eq.0d0) then
623
 
         if(jcentral(1).gt.0.and.jcentral(2).gt.0) then
624
 
            scale=(pt2ijcl(jcentral(1))*pt2ijcl(jcentral(2)))**0.25d0
625
 
         elseif(jcentral(1).gt.0) then
626
 
            scale=sqrt(pt2ijcl(jcentral(1)))
627
 
         elseif(jcentral(2).gt.0) then
628
 
            scale=sqrt(pt2ijcl(jcentral(2)))
629
 
         else
630
 
            scale=sqrt(pt2ijcl(nexternal-2))
631
 
         endif
632
 
         scale = scalefact*scale
633
 
         if(scale.gt.0)
634
 
     $        G = SQRT(4d0*PI*ALPHAS(scale))
635
 
      endif
636
 
      if (btest(mlevel,3))
637
 
     $     write(*,*) 'Set ren scale to ',scale
638
 
 
639
 
      if(ickkw.gt.0.and.q2fact(1).gt.0) then
640
 
c     Use the fixed or previously set scale for central scale
641
 
         if(jcentral(1).gt.0) pt2ijcl(jcentral(1))=q2fact(1)
642
 
         if(jcentral(2).gt.0.and.jcentral(2).ne.jcentral(1))
643
 
     $        pt2ijcl(jcentral(2))=q2fact(2)
644
 
      endif
645
 
 
646
 
      if(nexternal.eq.3.and.nincoming.eq.2.and.q2fact(1).eq.0) then
647
 
         q2fact(1)=pt2ijcl(nexternal-2)
648
 
         q2fact(2)=pt2ijcl(nexternal-2)
649
 
      endif
650
 
 
651
 
      if(q2fact(1).eq.0d0) then
652
 
c     Use the geom. average of central scale and first non-radiation vertex
653
 
         if(jlast(1).gt.0) q2fact(1)=sqrt(pt2ijcl(jlast(1))*pt2ijcl(jcentral(1)))
654
 
         if(jlast(2).gt.0) q2fact(2)=sqrt(pt2ijcl(jlast(2))*pt2ijcl(jcentral(2)))
655
 
         if(jcentral(1).gt.0.and.jcentral(1).eq.jcentral(2))then
656
 
c     We have a qcd line going through the whole event, use single scale
657
 
            q2fact(1)=max(q2fact(1),q2fact(2))
658
 
            q2fact(2)=q2fact(1)
659
 
         endif
660
 
      endif
661
 
      if(.not. fixed_fac_scale) then
662
 
         q2fact(1)=scalefact**2*q2fact(1)
663
 
         q2fact(2)=scalefact**2*q2fact(2)
664
 
         q2bck(1)=q2fact(1)
665
 
         q2bck(2)=q2fact(2)
666
 
         if (btest(mlevel,3))
667
 
     $      write(*,*) 'Set central fact scales to ',sqrt(q2bck(1)),sqrt(q2bck(2))
668
 
      endif
669
 
         
670
 
      if(jcentral(1).eq.0.and.jcentral(2).eq.0)then
671
 
         if(q2fact(1).gt.0)then
672
 
            pt2ijcl(nexternal-2)=q2fact(1)
673
 
            pt2ijcl(nexternal-3)=q2fact(1)
674
 
         else
675
 
            q2fact(1)=pt2ijcl(nexternal-2)
676
 
            q2fact(2)=q2fact(1)
677
 
         endif
678
 
      elseif(ickkw.eq.2.or.pdfwgt)then
679
 
c     Use the minimum scale found for fact scale in ME
680
 
         if(jlast(1).gt.0.and.jfirst(1).lt.jlast(1))
681
 
     $        q2fact(1)=min(pt2ijcl(jfirst(1)),q2fact(1))
682
 
         if(jlast(2).gt.0.and.jfirst(2).lt.jlast(2))
683
 
     $        q2fact(2)=min(pt2ijcl(jfirst(2)),q2fact(2))
684
 
      endif
685
 
 
686
 
c     Check that factorization scale is >= 2 GeV
687
 
      if(lpp(1).ne.0.and.q2fact(1).lt.4d0.or.
688
 
     $   lpp(2).ne.0.and.q2fact(2).lt.4d0)then
689
 
         if(nwarning.le.10) then
690
 
             nwarning=nwarning+1
691
 
             write(*,*) 'Warning: Too low fact scales: ',
692
 
     $            sqrt(q2fact(1)), sqrt(q2fact(2))
693
 
          endif
694
 
         if(nwarning.eq.11) then
695
 
             nwarning=nwarning+1
696
 
             write(*,*) 'No more warnings written out this run.'
697
 
          endif
698
 
         setclscales=.false.
699
 
         clustered = .false.
700
 
         return
701
 
      endif
702
 
 
703
 
      if (btest(mlevel,3))
704
 
     $     write(*,*) 'Set fact scales to ',sqrt(q2fact(1)),sqrt(q2fact(2))
705
 
      return
706
 
      end
707
 
      
708
 
      double precision function rewgt(p)
709
 
c**************************************************
710
 
c   reweight the hard me according to ckkw
711
 
c   employing the information in common/cl_val/
712
 
c**************************************************
713
 
      implicit none
714
 
 
715
 
      include 'message.inc'
716
 
      include 'genps.inc'
717
 
      include 'maxconfigs.inc'
718
 
      include 'nexternal.inc'
719
 
      include 'maxamps.inc'
720
 
      include 'cluster.inc'
721
 
      include 'run.inc'
722
 
      include 'coupl.inc'
723
 
C   
724
 
C   ARGUMENTS 
725
 
C   
726
 
      DOUBLE PRECISION P(0:3,NEXTERNAL)
727
 
 
728
 
C   global variables
729
 
C     Present process number
730
 
      INTEGER IMIRROR,IPROC
731
 
      COMMON/TO_MIRROR/IMIRROR, IPROC
732
 
      integer              IPSEL
733
 
      COMMON /SubProc/ IPSEL
734
 
      INTEGER SUBDIAG(MAXSPROC),IB(2)
735
 
      COMMON/TO_SUB_DIAG/SUBDIAG,IB
736
 
      data IB/1,2/
737
 
c     q2bck holds the central q2fact scales
738
 
      real*8 q2bck(2)
739
 
      common /to_q2bck/q2bck
740
 
      double precision stot,m1,m2
741
 
      common/to_stot/stot,m1,m2
742
 
 
743
 
C   local variables
744
 
      integer i, j, idi, idj
745
 
      real*8 PI
746
 
      parameter( PI = 3.14159265358979323846d0 )
747
 
 
748
 
      integer mapconfig(0:lmaxconfigs), this_config
749
 
      integer iforest(2,-max_branch:-1,lmaxconfigs)
750
 
      integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
751
 
      integer tprid(-max_branch:-1,lmaxconfigs)
752
 
      include 'configs.inc'
753
 
      real*8 xptj,xptb,xpta,xptl,xmtc
754
 
      real*8 xetamin,xqcut,deltaeta
755
 
      common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
756
 
      double precision asref, pt2prev(n_max_cl),pt2pdf(n_max_cl),pt2min
757
 
      integer n, ibeam(2), iqcd(0:2)!, ilast(0:nexternal)
758
 
      integer idfl, idmap(-nexternal:nexternal)
759
 
c     ipart gives external particle number chain
760
 
      integer ipart(2,n_max_cl)
761
 
      double precision xnow(2)
762
 
      double precision xtarget,tmp,pdfj1,pdfj2,q2now,etot
763
 
      integer iseed,np
764
 
      data iseed/0/
765
 
      logical isvx
766
 
 
767
 
      logical isqcd,isjet,isparton,isjetvx,ispartonvx
768
 
      double precision alphas,getissud,pdg2pdf, sudwgt
769
 
      real xran1
770
 
      external isqcd,isjet,isparton,ispartonvx
771
 
      external alphas, isjetvx, getissud, pdg2pdf, xran1,  sudwgt
772
 
 
773
 
      rewgt=1.0d0
774
 
      clustered=.false.
775
 
 
776
 
      if(ickkw.le.0) return
777
 
 
778
 
c   Set mimimum kt scale, depending on highest mult or not
779
 
      if(hmult.or.ickkw.eq.1)then
780
 
        pt2min=0
781
 
      else
782
 
        pt2min=xqcut**2
783
 
      endif
784
 
      if (btest(mlevel,3))
785
 
     $     write(*,*) 'pt2min set to ',pt2min
786
 
 
787
 
c   Set etot, used for non-radiating partons
788
 
      etot=sqrt(stot)
789
 
 
790
 
c   Since we use pdf reweighting, need to know particle identities
791
 
      if (btest(mlevel,1)) then
792
 
         write(*,*) 'Set process number ',ipsel
793
 
      endif
794
 
 
795
 
c   Preparing graph particle information (ipart, needed to keep track of
796
 
c   external particle clustering scales)
797
 
      do i=1,nexternal
798
 
c        ilast(i)=ishft(1,i)
799
 
         pt2prev(ishft(1,i-1))=0d0
800
 
         if (ickkw.eq.2) then
801
 
            if(pt2min.gt.0)then
802
 
               pt2prev(ishft(1,i-1))=
803
 
     $              max(pt2min,p(0,i)**2-p(1,i)**2-p(2,i)**2-p(3,i)**2)
804
 
            endif
805
 
            pt2pdf(ishft(1,i-1))=pt2prev(ishft(1,i-1))
806
 
         else if(pdfwgt) then
807
 
            pt2pdf(ishft(1,i-1))=0d0
808
 
         endif
809
 
         ptclus(i)=sqrt(pt2prev(ishft(1,i-1)))
810
 
         if (btest(mlevel,3))
811
 
     $        write(*,*) 'Set ptclus for ',i,' to ', ptclus(i)
812
 
         ipart(1,ishft(1,i-1))=i
813
 
         ipart(2,ishft(1,i-1))=0
814
 
         if (btest(mlevel,4))
815
 
     $        write(*,*) 'Set ipart for ',ishft(1,i-1),' to ',
816
 
     $        ipart(1,ishft(1,i-1)),ipart(2,ishft(1,i-1))
817
 
      enddo
818
 
c      ilast(0)=nexternal
819
 
      ibeam(1)=ishft(1,0)
820
 
      ibeam(2)=ishft(1,1)
821
 
      if (btest(mlevel,1)) then
822
 
        write(*,*)'rewgt: identified tree {'
823
 
        do i=1,nexternal-2
824
 
          write(*,*)'  ',i,': ',idacl(i,1),'(',ipdgcl(idacl(i,1),igraphs(1),iproc),')',
825
 
     $       '&',idacl(i,2),'(',ipdgcl(idacl(i,2),igraphs(1),iproc),')',
826
 
     $       ' -> ',imocl(i),'(',ipdgcl(imocl(i),igraphs(1),iproc),')',
827
 
     $       ', ptij = ',dsqrt(pt2ijcl(i))
828
 
        enddo
829
 
        write(*,*)'  graphs (',igraphs(0),'):',(igraphs(i),i=1,igraphs(0))
830
 
        write(*,*)'}'
831
 
      endif
832
 
c     Set x values for the two sides, for IS Sudakovs
833
 
      do i=1,2
834
 
        xnow(i)=xbk(ib(i))
835
 
      enddo
836
 
      if(btest(mlevel,3))then
837
 
        write(*,*) 'Set x values to ',xnow(1),xnow(2)
838
 
      endif
839
 
 
840
 
c     Prepare for resetting q2fact based on PDF reweighting
841
 
      if(ickkw.eq.2)then
842
 
         q2fact(1)=0d0
843
 
         q2fact(2)=0d0
844
 
      endif
845
 
c   
846
 
c   Set strong coupling used
847
 
c   
848
 
      asref=G**2/(4d0*PI)
849
 
 
850
 
c   Perform alpha_s reweighting based on type of vertex
851
 
      do n=1,nexternal-2
852
 
c       scale for alpha_s reweighting
853
 
        q2now=pt2ijcl(n)
854
 
        if(n.eq.nexternal-2) then
855
 
           q2now = scale**2
856
 
        endif
857
 
        if (btest(mlevel,3)) then
858
 
          write(*,*)'  ',n,': ',idacl(n,1),'(',ipdgcl(idacl(n,1),igraphs(1),iproc),
859
 
     &       ')&',idacl(n,2),'(',ipdgcl(idacl(n,2),igraphs(1),iproc),
860
 
     &       ') -> ',imocl(n),'(',ipdgcl(imocl(n),igraphs(1),iproc),
861
 
     &       '), ptij = ',dsqrt(q2now) 
862
 
        endif
863
 
c     perform alpha_s reweighting only for vertices where a parton is produced
864
 
c     and not for the last clustering (use non-fixed ren. scale for these)
865
 
        if (n.lt.nexternal-2)then
866
 
           if(ispartonvx(imocl(n),idacl(n,1),idacl(n,2),
867
 
     $       ipdgcl(1,igraphs(1),iproc),ipart,.false.)) then
868
 
c       alpha_s weight
869
 
              rewgt=rewgt*alphas(alpsfact*sqrt(q2now))/asref
870
 
              if (btest(mlevel,3)) then
871
 
                 write(*,*)' reweight vertex: ',ipdgcl(imocl(n),igraphs(1),iproc),
872
 
     $                ipdgcl(idacl(n,1),igraphs(1),iproc),ipdgcl(idacl(n,2),igraphs(1),iproc)
873
 
                 write(*,*)'       as: ',alphas(alpsfact*dsqrt(q2now)),
874
 
     &                '/',asref,' -> ',alphas(alpsfact*dsqrt(q2now))/asref
875
 
                 write(*,*)' and G=',SQRT(4d0*PI*ALPHAS(scale))
876
 
              endif
877
 
           endif
878
 
        endif
879
 
c   Update starting values for FS parton showering
880
 
        do i=1,2
881
 
          do j=1,2
882
 
            if(ipart(j,idacl(n,i)).gt.0.and.ipart(j,idacl(n,i)).gt.2)then
883
 
              ptclus(ipart(j,idacl(n,i)))=
884
 
     $              max(ptclus(ipart(j,idacl(n,i))),dsqrt(pt2ijcl(n)))
885
 
              if(ickkw.ne.2.and.
886
 
     $             (.not.isqcd(ipdgcl(imocl(n),igraphs(1),iproc)).or.
887
 
     $             ipart(1,idacl(n,3-i)).le.2.and.
888
 
     $             .not.isqcd(ipdgcl(idacl(n,3-i),igraphs(1),iproc)).or.
889
 
     $             isbw(imocl(n))))then
890
 
c             For particles originating in non-qcd t-channel vertices or decay vertices,
891
 
c             set origination scale to machine energy since we don't want these
892
 
c             to be included in matching.
893
 
                 ptclus(ipart(j,idacl(n,i)))=etot
894
 
              endif
895
 
              if (btest(mlevel,3))
896
 
     $             write(*,*) 'Set ptclus for ',ipart(j,idacl(n,i)),
897
 
     $             ' to ', ptclus(ipart(j,idacl(n,i))),
898
 
     $             ipdgcl(imocl(n),igraphs(1),iproc),
899
 
     $             isqcd(ipdgcl(imocl(n),igraphs(1),iproc)),isbw(imocl(n))
900
 
            endif
901
 
          enddo
902
 
        enddo
903
 
c     Special case for last 1,2->i vertex
904
 
        if(n.eq.nexternal-2)then
905
 
           ptclus(ipart(1,imocl(n)))=
906
 
     $              max(ptclus(ipart(1,imocl(n))),dsqrt(pt2ijcl(n)))
907
 
              if(ickkw.ne.2.and.
908
 
     $             (.not.isqcd(ipdgcl(idacl(n,1),igraphs(1),iproc)).or.
909
 
     $             .not.isqcd(ipdgcl(idacl(n,2),igraphs(1),iproc))))then
910
 
c             For particles originating in non-qcd vertices or decay vertices,
911
 
c             set origination scale to machine energy since we don't want these
912
 
c             to be included in matching.
913
 
                 ptclus(ipart(1,imocl(n)))=etot
914
 
              endif
915
 
              if (btest(mlevel,3))
916
 
     $             write(*,*) 'Set ptclus for ',ipart(1,imocl(n)),
917
 
     $             ' to ', ptclus(ipart(1,imocl(n))),
918
 
     $             ipdgcl(idacl(n,1),igraphs(1),iproc),
919
 
     $             ipdgcl(idacl(n,2),igraphs(1),iproc)
920
 
        endif
921
 
c   Update particle tree map
922
 
        call ipartupdate(p,imocl(n),idacl(n,1),idacl(n,2),
923
 
     $       ipdgcl(1,igraphs(1),iproc),ipart)
924
 
        if(ickkw.eq.2.or.pdfwgt) then
925
 
c       Perform PDF and, if ickkw=2, Sudakov reweighting
926
 
          isvx=.false.
927
 
          do i=1,2
928
 
c         write(*,*)'weight ',idacl(n,i),', ptij=',pt2prev(idacl(n,i))
929
 
            if (isqcd(ipdgcl(idacl(n,i),igraphs(1),iproc))) then
930
 
               if(ickkw.eq.2.and.pt2min.eq.0d0) then
931
 
                  pt2min=pt2ijcl(n)
932
 
                  if (btest(mlevel,3))
933
 
     $                 write(*,*) 'pt2min set to ',pt2min
934
 
               endif
935
 
               if(ickkw.eq.2.and.pt2prev(idacl(n,i)).eq.0d0)
936
 
     $              pt2prev(idacl(n,i))=
937
 
     $              max(pt2min,p(0,i)**2-p(1,i)**2-p(2,i)**2-p(3,i)**2)
938
 
               do j=1,2
939
 
                  if (isparton(ipdgcl(idacl(n,i),igraphs(1),iproc)).and
940
 
     $                 .idacl(n,i).eq.ibeam(j)) then
941
 
c               is sudakov weight - calculate only once for each parton
942
 
c               line where parton line ends with change of parton id or
943
 
c               non-radiation vertex
944
 
                     isvx=.true.
945
 
                     ibeam(j)=imocl(n)
946
 
c                    Perform Sudakov reweighting if ickkw=2
947
 
                     if(ickkw.eq.2.and.(ipdgcl(idacl(n,i),igraphs(1),iproc).ne.
948
 
     $                    ipdgcl(imocl(n),igraphs(1),iproc).or.
949
 
     $                    .not.isjetvx(imocl(n),idacl(n,1),idacl(n,2),
950
 
     $                    ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)).and.
951
 
     $                    pt2prev(idacl(n,i)).lt.pt2ijcl(n))then
952
 
                        tmp=min(1d0,max(getissud(ibeam(j),ipdgcl(idacl(n,i),
953
 
     $                       igraphs(1),iproc),xnow(j),xnow(3-j),pt2ijcl(n)),1d-20)/
954
 
     $                       max(getissud(ibeam(j),ipdgcl(idacl(n,i),
955
 
     $                       igraphs(1),iproc),xnow(j),xnow(3-j),pt2prev(idacl(n,i))),1d-20))
956
 
                        rewgt=rewgt*tmp
957
 
                        pt2prev(imocl(n))=pt2ijcl(n)
958
 
                        if (btest(mlevel,3)) then
959
 
                           write(*,*)' reweight line: ',ipdgcl(idacl(n,i),igraphs(1),iproc), idacl(n,i)
960
 
                           write(*,*)'     pt2prev, pt2new, x1, x2: ',pt2prev(idacl(n,i)),pt2ijcl(n),xnow(j),xnow(3-j)
961
 
                           write(*,*)'           Sud: ',tmp
962
 
                           write(*,*)'        -> rewgt: ',rewgt
963
 
                        endif
964
 
                     else if(ickkw.eq.2) then
965
 
                        pt2prev(imocl(n))=pt2prev(idacl(n,i))
966
 
                     endif
967
 
c                 End Sudakov reweighting when we reach a non-radiation vertex
968
 
                     if(ickkw.eq.2.and..not.
969
 
     $                    ispartonvx(imocl(n),idacl(n,1),idacl(n,2),
970
 
     $                    ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)) then
971
 
                        pt2prev(imocl(n))=1d30
972
 
                        if (btest(mlevel,3)) then
973
 
                          write(*,*)' rewgt: ending reweighting for vx ',
974
 
     $                          idacl(n,1),idacl(n,2),imocl(n),
975
 
     $                          ' with ids ',ipdgcl(idacl(n,1),igraphs(1),iproc),
976
 
     $                          ipdgcl(idacl(n,2),igraphs(1),iproc),ipdgcl(imocl(n),igraphs(1),iproc)
977
 
                        endif
978
 
                     endif
979
 
c               PDF reweighting
980
 
c               Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
981
 
c               f1(x1,pt2E) is given by DSIG, already set scale for that
982
 
                     if (zcl(n).gt.0d0.and.zcl(n).lt.1d0) then
983
 
                        xnow(j)=xnow(j)*zcl(n)
984
 
                     endif
985
 
c                    PDF scale
986
 
                     q2now=min(pt2ijcl(n), q2bck(j))
987
 
c                    Set PDF scale to central factorization scale
988
 
c                    if non-radiating vertex or last 2->2
989
 
                     if(.not.isjetvx(imocl(n),idacl(n,1),idacl(n,2),
990
 
     $                    ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)) then
991
 
                        q2now=q2bck(j)
992
 
                     endif
993
 
                     if (btest(mlevel,3))
994
 
     $                    write(*,*)' set q2now for pdf to ',sqrt(q2now)
995
 
                     if(q2fact(j).eq.0d0.and.ickkw.eq.2)then
996
 
                        q2fact(j)=pt2min ! Starting scale for PS
997
 
                        pt2pdf(imocl(n))=q2now
998
 
                        if (btest(mlevel,3))
999
 
     $                       write(*,*)' set fact scale ',j,
1000
 
     $                          ' for PS scale to: ',sqrt(q2fact(j))
1001
 
                     else if(pt2pdf(idacl(n,i)).eq.0d0)then
1002
 
                        pt2pdf(imocl(n))=q2now
1003
 
                        if (btest(mlevel,3))
1004
 
     $                       write(*,*)' set pt2pdf for ',imocl(n),
1005
 
     $                          ' to: ',sqrt(pt2pdf(imocl(n)))
1006
 
                     else if(pt2pdf(idacl(n,i)).lt.q2now
1007
 
     $                       .and.isjet(ipdgcl(idacl(n,i),igraphs(1),iproc))) then
1008
 
                        pdfj1=pdg2pdf(abs(lpp(IB(j))),ipdgcl(idacl(n,i),
1009
 
     $                       igraphs(1),iproc)*sign(1,lpp(IB(j))),
1010
 
     $                       xnow(j),sqrt(q2now))
1011
 
                        pdfj2=pdg2pdf(abs(lpp(IB(j))),ipdgcl(idacl(n,i),
1012
 
     $                       igraphs(1),iproc)*sign(1,lpp(IB(j))),
1013
 
     $                       xnow(j),sqrt(pt2pdf(idacl(n,i))))
1014
 
                        if(pdfj2.lt.1d-10)then
1015
 
c                          Scale too low for heavy quark
1016
 
                           rewgt=0d0
1017
 
                           if (btest(mlevel,3))
1018
 
     $                        write(*,*) 'Too low scale for quark pdf: ',
1019
 
     $                        sqrt(pt2pdf(idacl(n,i))),pdfj2,pdfj1
1020
 
                           return
1021
 
                        endif
1022
 
                        rewgt=rewgt*pdfj1/pdfj2
1023
 
                        if (btest(mlevel,3)) then
1024
 
                           write(*,*)' reweight ',n,i,ipdgcl(idacl(n,i),igraphs(1),iproc),' by pdfs: '
1025
 
                           write(*,*)'     x, ptprev, ptnew: ',xnow(j),
1026
 
     $                          sqrt(pt2pdf(idacl(n,i))),sqrt(q2now)
1027
 
                           write(*,*)'           PDF: ',pdfj1,' / ',pdfj2
1028
 
                           write(*,*)'        -> rewgt: ',rewgt
1029
 
c                           write(*,*)'  (compare for glue: ',
1030
 
c     $                          pdg2pdf(lpp(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i)))),' / ',
1031
 
c     $                          pdg2pdf(lpp(j),21,xbk(j),sqrt(pt2ijcl(n)))
1032
 
c                           write(*,*)'       = ',pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i))))/
1033
 
c     $                          pdg2pdf(lpp(j),21,xbk(j),sqrt(pt2ijcl(n)))
1034
 
c                           write(*,*)'       -> ',pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i))))/
1035
 
c     $                          pdg2pdf(lpp(j),21,xbk(j),sqrt(pt2ijcl(n)))*rewgt,' )'
1036
 
                        endif
1037
 
c                       Set scale for mother as this scale
1038
 
                        pt2pdf(imocl(n))=q2now                           
1039
 
                     else if(pt2pdf(idacl(n,i)).ge.q2now) then
1040
 
c                    If no reweighting, just copy daughter scale for mother
1041
 
                        pt2pdf(imocl(n))=pt2pdf(idacl(n,i))
1042
 
                     endif
1043
 
                     goto 10
1044
 
                  endif
1045
 
               enddo
1046
 
c           fs sudakov weight
1047
 
               if(ickkw.eq.2.and.pt2prev(idacl(n,i)).lt.pt2ijcl(n).and.
1048
 
     $              (isvx.or.ipdgcl(idacl(n,i),igraphs(1),iproc).ne.ipdgcl(imocl(n),igraphs(1),iproc).or.
1049
 
     $              (ipdgcl(idacl(n,i),igraphs(1),iproc).ne.
1050
 
     $              ipdgcl(idacl(n,3-i),igraphs(1),iproc).and.
1051
 
     $              pt2prev(idacl(n,i)).gt.pt2prev(idacl(n,3-i))))) then
1052
 
                  tmp=sudwgt(sqrt(pt2min),sqrt(pt2prev(idacl(n,i))),
1053
 
     &                 dsqrt(pt2ijcl(n)),ipdgcl(idacl(n,i),igraphs(1),iproc),1)
1054
 
                  rewgt=rewgt*tmp
1055
 
                  if (btest(mlevel,3)) then
1056
 
                     write(*,*)' reweight fs line: ',ipdgcl(idacl(n,i),igraphs(1),iproc), idacl(n,i)
1057
 
                     write(*,*)'     pt2prev, pt2new: ',pt2prev(idacl(n,i)),pt2ijcl(n)
1058
 
                     write(*,*)'           Sud: ',tmp
1059
 
                     write(*,*)'        -> rewgt: ',rewgt
1060
 
                  endif
1061
 
                  pt2prev(imocl(n))=pt2ijcl(n)
1062
 
               else
1063
 
                  pt2prev(imocl(n))=pt2prev(idacl(n,i))
1064
 
               endif 
1065
 
            endif
1066
 
 10         continue
1067
 
          enddo
1068
 
          if (ickkw.eq.2.and.n.eq.nexternal-2.and.isqcd(ipdgcl(imocl(n),igraphs(1),iproc)).and.
1069
 
     $         pt2prev(imocl(n)).lt.pt2ijcl(n)) then
1070
 
             tmp=sudwgt(sqrt(pt2min),sqrt(pt2prev(imocl(n))),
1071
 
     &            dsqrt(pt2ijcl(n)),ipdgcl(imocl(n),igraphs(1),iproc),1)
1072
 
             rewgt=rewgt*tmp
1073
 
             if (btest(mlevel,3)) then
1074
 
                write(*,*)' reweight last fs line: ',ipdgcl(imocl(n),igraphs(1),iproc), imocl(n)
1075
 
                write(*,*)'     pt2prev, pt2new: ',pt2prev(imocl(n)),pt2ijcl(n)
1076
 
                write(*,*)'           Sud: ',tmp
1077
 
                write(*,*)'        -> rewgt: ',rewgt
1078
 
             endif
1079
 
          endif
1080
 
        endif
1081
 
      enddo
1082
 
 
1083
 
      if(ickkw.eq.2.and.lpp(1).eq.0.and.lpp(2).eq.0)then
1084
 
         q2fact(1)=pt2min
1085
 
         q2fact(2)=q2fact(1)
1086
 
      else if (ickkw.eq.1.and.pdfwgt) then
1087
 
         q2fact(1)=q2bck(1)
1088
 
         q2fact(2)=q2bck(2)         
1089
 
         if (btest(mlevel,3))
1090
 
     $        write(*,*)' set fact scales for PS to ',
1091
 
     $        sqrt(q2fact(1)),sqrt(q2fact(2))
1092
 
      endif
1093
 
 
1094
 
      if (btest(mlevel,3)) then
1095
 
        write(*,*)'} ->  w = ',rewgt
1096
 
      endif
1097
 
      return
1098
 
      end
1099