← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_141512/PROC_141512/SubProcesses/cluster.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_141512/PROC_141512/SubProcesses/cluster.f
 
1
      subroutine crossp(p1,p2,p)
 
2
c**************************************************************************
 
3
c     input:
 
4
c            p1, p2    vectors to cross
 
5
c**************************************************************************
 
6
      implicit none
 
7
      real*8 p1(0:3), p2(0:3), p(0:3)
 
8
 
 
9
      p(0)=0d0
 
10
      p(1)=p1(2)*p2(3)-p1(3)*p2(2)
 
11
      p(2)=p1(3)*p2(1)-p1(1)*p2(3)
 
12
      p(3)=p1(1)*p2(2)-p1(2)*p2(1)
 
13
 
 
14
      return 
 
15
      end
 
16
 
 
17
 
 
18
      subroutine rotate(p1,p2,n,nn2,ct,st,d)
 
19
c**************************************************************************
 
20
c     input:
 
21
c            p1        vector to be rotated
 
22
c            n         vector perpendicular to plane of rotation
 
23
c            nn2       squared norm of n to improve numerics
 
24
c            ct, st    cos/sin theta of rotation in plane 
 
25
c            d         direction: 1 there / -1 back
 
26
c     output:
 
27
c            p2        p1 rotated using defined rotation
 
28
c**************************************************************************
 
29
      implicit none
 
30
      real*8 p1(0:3), p2(0:3), n(0:3), at(0:3), ap(0:3), cr(0:3)
 
31
      double precision nn2, ct, st, na, nn
 
32
      integer d, i
 
33
 
 
34
      if (nn2.eq.0d0) then
 
35
         do i=0,3
 
36
            p2(i)=p1(i)
 
37
         enddo   
 
38
         return
 
39
      endif
 
40
      nn=dsqrt(nn2)
 
41
      na=(n(1)*p1(1)+n(2)*p1(2)+n(3)*p1(3))/nn2
 
42
      do i=1,3
 
43
         at(i)=n(i)*na
 
44
         ap(i)=p1(i)-at(i)
 
45
      enddo
 
46
c      write(*,*)'nn2 ',nn2,' ',nn,' ',na
 
47
c      write(*,*)'ap ',ap(1),',',ap(2),',',ap(3)
 
48
c      write(*,*)'at ',at(1),',',at(2),',',at(3)
 
49
      p2(0)=p1(0)
 
50
      call crossp(n,ap,cr)
 
51
c      write(*,*)'cr ',cr(1),',',cr(2),',',cr(3)
 
52
      do i=1,3
 
53
         if (d.ge.0) then
 
54
            p2(i)=at(i)+ct*ap(i)+st/nn*cr(i)
 
55
         else 
 
56
            p2(i)=at(i)+ct*ap(i)-st/nn*cr(i)
 
57
         endif
 
58
      enddo
 
59
      
 
60
      return 
 
61
      end
 
62
 
 
63
 
 
64
      subroutine constr(p1,p2,n,nn2,ct,st)
 
65
c**************************************************************************
 
66
c     input:
 
67
c            p1, p2    p1 rotated onto p2 defines plane of rotation
 
68
c     output:
 
69
c            n         vector perpendicular to plane of rotation
 
70
c            nn2       squared norm of n to improve numerics
 
71
c            ct, st    cos/sin theta of rotation in plane 
 
72
c**************************************************************************
 
73
      implicit none
 
74
      real*8 p1(0:3), p2(0:3), n(0:3), tr(0:3)
 
75
      double precision nn2, ct, st, mct
 
76
 
 
77
      ct=p1(1)*p2(1)+p1(2)*p2(2)+p1(3)*p2(3)
 
78
      ct=ct/dsqrt(p1(1)**2+p1(2)**2+p1(3)**2)
 
79
      ct=ct/dsqrt(p2(1)**2+p2(2)**2+p2(3)**2)
 
80
      mct=ct
 
81
c     catch bad numerics
 
82
      if (mct-1d0>0d0) mct=0d0
 
83
      st=dsqrt(1d0-mct*mct)
 
84
 
 
85
      call crossp(p1,p2,n)
 
86
      nn2=n(1)**2+n(2)**2+n(3)**2
 
87
c     don't rotate if nothing to rotate
 
88
      if (nn2.le.1d-34) then
 
89
         nn2=0d0
 
90
         return
 
91
      endif
 
92
 
 
93
c     check rotation
 
94
c      call rotate(p1(0),tr(0),n(0),nn2,ct,st,1)
 
95
c      write(*,*)'p1 (',p1(0),',',p1(1),',',p1(2),',',p1(3),')'
 
96
c      write(*,*)'p2 (',p2(0),',',p2(1),',',p2(2),',',p2(3),')'
 
97
c      write(*,*)'nn (',n(0),',',n(1),',',n(2),',',n(3),')'
 
98
c      write(*,*)'nn (',n(0),',',n(1),',',n(2),',',n(3),')'
 
99
c      write(*,*)'nn2 = ',nn2,', ct = ',ct,', st = ',st
 
100
c      write(*,*)'tr (',tr(0),',',tr(1),',',tr(2),',',tr(3),')'
 
101
      
 
102
      return 
 
103
      end
 
104
 
 
105
 
 
106
      Subroutine mapids(ids,id)
 
107
c**************************************************************************
 
108
c     input:
 
109
c            ids       array of particle ids
 
110
c            id        compressed particle id
 
111
c**************************************************************************
 
112
      implicit none
 
113
      include 'nexternal.inc'
 
114
      integer i, id, ids(nexternal)
 
115
      
 
116
      id=0
 
117
      do i=1,nexternal
 
118
         if (ids(i).ne.0) then
 
119
            id=id+ishft(1,i)
 
120
         endif
 
121
      enddo
 
122
c      write(*,*) 'cluster.f: compressed code is ',id
 
123
 
 
124
      return
 
125
      end
 
126
 
 
127
 
 
128
      subroutine mapid(id,ids)
 
129
c**************************************************************************
 
130
c     input:
 
131
c            id        compressed particle id
 
132
c            ids       array of particle ids
 
133
c**************************************************************************
 
134
      implicit none
 
135
      include 'nexternal.inc'
 
136
      integer i, icd, id, ids(nexternal)
 
137
      
 
138
      icd=id
 
139
      do i=1,nexternal
 
140
         ids(i)=0
 
141
         if (btest(id,i)) then
 
142
            ids(i)=1
 
143
         endif
 
144
c         write(*,*) 'cluster.f: uncompressed code ',i,' is ',ids(i)
 
145
      enddo
 
146
 
 
147
      return
 
148
      end
 
149
 
 
150
 
 
151
      integer function combid(i,j)
 
152
c**************************************************************************
 
153
c     input:
 
154
c            i,j       legs to combine
 
155
c     output:
 
156
c            index of combined leg
 
157
c**************************************************************************
 
158
      implicit none
 
159
      include 'nexternal.inc'
 
160
      integer i, j
 
161
 
 
162
c      combid=min(i+j,ishft(1,nexternal+1)-2-i-j)
 
163
      combid = i+j
 
164
     
 
165
      return
 
166
      end
 
167
 
 
168
 
 
169
      subroutine filprp(iproc,ignum,idij)
 
170
c**************************************************************************
 
171
c     Include graph ignum in list for propagator idij
 
172
c**************************************************************************
 
173
      implicit none
 
174
      include 'nexternal.inc'
 
175
      include 'maxamps.inc'
 
176
      include 'cluster.inc'
 
177
      include 'message.inc'
 
178
      integer ignum, idij, iproc, i
 
179
 
 
180
      if(idij.gt.n_max_cl) return
 
181
      do i=1,id_cl(iproc,idij,0)
 
182
         if (id_cl(iproc,idij,i).eq.ignum) return
 
183
      enddo
 
184
      id_cl(iproc,idij,0)=id_cl(iproc,idij,0)+1
 
185
      id_cl(iproc,idij,id_cl(iproc,idij,0))=ignum
 
186
      if(btest(mlevel,5))
 
187
     $     write(*,*)'Adding graph ',ignum,' to prop ',idij,' for proc ',iproc
 
188
      return
 
189
      end
 
190
 
 
191
      logical function filgrp(ignum,ipnum,ipids)
 
192
c**************************************************************************
 
193
c     input:
 
194
c            ignum      number of graph to be analysed
 
195
c            ipnum      number of level to be analysed, 
 
196
c                       starting with nexternal
 
197
c            ipids      particle number, iforest number, 
 
198
c                       daughter1, daughter2
 
199
c     output:
 
200
c            true if no errors
 
201
c**************************************************************************
 
202
      implicit none
 
203
      include 'genps.inc'
 
204
      include 'maxconfigs.inc'
 
205
      include 'nexternal.inc'
 
206
      include 'maxamps.inc'
 
207
      include 'cluster.inc'
 
208
      include 'coupl.inc'
 
209
      include 'message.inc'
 
210
      integer ignum, ipnum, ipids(nexternal,4,2:nexternal)
 
211
C $B$ IFOREST $B$ !this is a tag for MadWeight
 
212
      integer i, iforest(2,-max_branch:-1,lmaxconfigs)
 
213
      common/to_forest/iforest
 
214
      integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
 
215
      integer tprid(-max_branch:-1,lmaxconfigs)
 
216
      common/to_sprop/sprop,tprid
 
217
C $E$ IFOREST $E$ !this is a tag for MadWeight
 
218
      INTEGER    n_max_cl_cg
 
219
      PARAMETER (n_max_cl_cg=n_max_cl*n_max_cg)
 
220
      data resmap/n_max_cl_cg*.false./
 
221
 
 
222
      Integer j, k, l, icmp(2), iproc
 
223
 
 
224
      double precision ZERO
 
225
      parameter (ZERO=0d0)
 
226
      double precision prmass(-nexternal:0,lmaxconfigs)
 
227
      double precision prwidth(-nexternal:0,lmaxconfigs)
 
228
      integer pow(-nexternal:0,lmaxconfigs)
 
229
      logical first_time
 
230
      save prmass,prwidth,pow
 
231
      INTEGER CONFSUB(MAXSPROC,LMAXCONFIGS)
 
232
      INCLUDE 'config_subproc_map.inc'
 
233
      data first_time /.true./
 
234
 
 
235
      integer combid
 
236
      logical isjet
 
237
      external combid,isjet
 
238
 
 
239
      if (first_time) then
 
240
         include 'props.inc'
 
241
         first_time=.false.
 
242
      endif
 
243
 
 
244
      if(btest(mlevel,4))
 
245
     $     write(*,*) 'graph,level: ',ignum,ipnum
 
246
 
 
247
      filgrp=.false.
 
248
C   Follow diagram tree down to last clustering
 
249
      do i=1,ipnum
 
250
         do j=i+1,ipnum
 
251
            if(btest(mlevel,4))
 
252
     $           write(*,*)'at ids   (',ipids(i,1,ipnum),',',ipids(i,2,ipnum),'), (',
 
253
     $           ipids(j,1,ipnum),',',ipids(j,2,ipnum),'), ',i,j
 
254
            do k=-nexternal+1,-1
 
255
               if ((iforest(1,k,ignum).eq.ipids(i,2,ipnum).and.
 
256
     &              iforest(2,k,ignum).eq.ipids(j,2,ipnum)).or.
 
257
     &             (iforest(2,k,ignum).eq.ipids(i,2,ipnum).and.
 
258
     &              iforest(1,k,ignum).eq.ipids(j,2,ipnum))) then
 
259
c                 Add the combined propagator
 
260
                  icmp(1)=combid(ipids(i,1,ipnum),ipids(j,1,ipnum))
 
261
c                 Add also the same propagator but from the other direction
 
262
                  icmp(2)=ishft(1,nexternal)-1-icmp(1)
 
263
c     Set pdg code for propagator
 
264
                  do l=1,2
 
265
                     do iproc=1,maxsproc
 
266
                     if(confsub(iproc,ignum).eq.0) cycle
 
267
                     if(sprop(iproc,k,ignum).ne.0)then
 
268
                        ipdgcl(icmp(l),ignum,iproc)=sprop(iproc,k,ignum)
 
269
c                       If this is radiation off heavy FS particle, set heavyrad to true
 
270
                        if(isjet(ipdgcl(ipids(i,1,ipnum),ignum,iproc)).and.
 
271
     $                       .not.isjet(ipdgcl(ipids(j,1,ipnum),ignum,iproc)).and.
 
272
     $                       ipdgcl(ipids(j,1,ipnum),ignum,iproc).eq.sprop(iproc,k,ignum).or.
 
273
     $                       isjet(ipdgcl(ipids(j,1,ipnum),ignum,iproc)).and.
 
274
     $                       .not.isjet(ipdgcl(ipids(i,1,ipnum),ignum,iproc)).and.
 
275
     $                       ipdgcl(ipids(i,1,ipnum),ignum,iproc).eq.sprop(iproc,k,ignum))then
 
276
                           heavyrad(ignum) = .true.
 
277
                        endif
 
278
                     else if(tprid(k,ignum).ne.0)then
 
279
                        ipdgcl(icmp(l),ignum,iproc)=tprid(k,ignum)
 
280
                     else if(ipnum.eq.3)then
 
281
                        ipdgcl(icmp(l),ignum,iproc)=ipdgcl(2,ignum,iproc)
 
282
                     else
 
283
                        ipdgcl(icmp(l),ignum,iproc)=0
 
284
                     endif
 
285
                     if(btest(mlevel,4))
 
286
     $                    write(*,*) 'add table entry for (',ipids(i,1,ipnum),
 
287
     &                    ',',ipids(j,1,ipnum),',',icmp(l),')',
 
288
     $                    'proc: ',iproc,
 
289
     $                    ', pdg: ',ipdgcl(icmp(l),ignum,iproc)
 
290
                     call filprp(iproc,ignum,icmp(l))
 
291
c               Insert graph in list of propagators
 
292
                     if(prwidth(k,ignum).gt.ZERO) then
 
293
                     if(btest(mlevel,4))
 
294
     $                       write(*,*)'Adding resonance ',ignum,icmp(l)
 
295
                        resmap(icmp(l),ignum)=.true.
 
296
                     endif
 
297
                  enddo
 
298
                  enddo
 
299
c     proceed w/ next table, since there is no possibility,
 
300
c     to combine the same particle in another way in this graph
 
301
                  ipids(i,1,ipnum-1)=icmp(1)
 
302
                  ipids(i,2,ipnum-1)=k
 
303
                  ipids(i,3,ipnum-1)=i
 
304
                  ipids(i,4,ipnum-1)=j
 
305
                  ipnum=ipnum-1
 
306
                  do l=1,j-1
 
307
                    if(l.eq.i) cycle
 
308
                    ipids(l,1,ipnum)=ipids(l,1,ipnum+1)
 
309
                    ipids(l,2,ipnum)=ipids(l,2,ipnum+1)
 
310
                    ipids(l,3,ipnum)=l
 
311
                    ipids(l,4,ipnum)=0
 
312
                  enddo
 
313
                  do l=j,ipnum
 
314
                     ipids(l,1,ipnum)=ipids(l+1,1,ipnum+1)
 
315
                     ipids(l,2,ipnum)=ipids(l+1,2,ipnum+1)
 
316
                     ipids(l,3,ipnum)=l+1
 
317
                     ipids(l,4,ipnum)=0
 
318
                  enddo
 
319
                  if(ipnum.eq.2)then
 
320
c                 Done with this diagram
 
321
                     return
 
322
                  else
 
323
                     filgrp=.true.
 
324
                     return
 
325
                  endif
 
326
               endif
 
327
            enddo
 
328
         enddo
 
329
      enddo
 
330
      return
 
331
      end
 
332
 
 
333
 
 
334
      logical function filmap()
 
335
c**************************************************************************
 
336
c     output:
 
337
c            true if no errors
 
338
c**************************************************************************
 
339
      implicit none
 
340
      include 'genps.inc'
 
341
      include 'maxconfigs.inc'
 
342
      include 'nexternal.inc'
 
343
      include 'maxamps.inc'
 
344
      include 'cluster.inc'
 
345
      include 'run.inc'
 
346
      include 'message.inc'
 
347
C $B$ IFOREST $B$ !this is a tag for MadWeight
 
348
      integer mapconfig(0:lmaxconfigs), this_config
 
349
      common/to_mconfigs/mapconfig, this_config
 
350
C $E$ IFOREST $E$ !this is a tag for MadWeight
 
351
      integer i, j, inpids, iproc, ipids(nexternal,4,2:nexternal)
 
352
      integer start_config,end_config
 
353
      integer idup(nexternal,maxproc,maxsproc)
 
354
      integer mothup(2,nexternal)
 
355
      integer icolup(2,nexternal,maxflow,maxsproc)
 
356
      include 'leshouche.inc'
 
357
      
 
358
      logical filgrp
 
359
      external filgrp
 
360
 
 
361
      if(chcluster) then
 
362
         start_config=this_config
 
363
         end_config=this_config
 
364
      else
 
365
         start_config=1
 
366
         end_config=mapconfig(0)
 
367
      endif
 
368
      do iproc=1,maxsproc
 
369
         do i=1,n_max_cl
 
370
            id_cl(iproc,i,0)=0
 
371
         enddo
 
372
      enddo
 
373
      do i=start_config,end_config
 
374
         heavyrad(i)=.false.
 
375
c         write (*,*) ' at graph ',i
 
376
         do j=1,nexternal
 
377
            ipids(j,1,nexternal)=ishft(1,j-1)
 
378
            ipids(j,2,nexternal)=j
 
379
            ipids(j,3,nexternal)=0
 
380
            ipids(j,4,nexternal)=0
 
381
            do iproc=1,maxsproc
 
382
               ipdgcl(ipids(j,1,nexternal),i,iproc)=idup(j,1,iproc)
 
383
            enddo
 
384
         enddo
 
385
         inpids=nexternal
 
386
c         print *,'Inserting graph ',i
 
387
 10      if (filgrp(i,inpids,ipids)) goto 10
 
388
         if(btest(mlevel,4).and.heavyrad(i)) then
 
389
            write(*,*)' set heavyrad of ',i,' to T'
 
390
         endif
 
391
      enddo
 
392
c     Ensure that there are some allowed clusterings
 
393
      do i=start_config,end_config
 
394
         if(.not.heavyrad(i)) goto 20
 
395
      enddo
 
396
      if(btest(mlevel,4)) then
 
397
         write(*,*)' Reset all heavyrad to .false.'
 
398
      endif      
 
399
      do i=start_config,end_config
 
400
         heavyrad(i)=.false.
 
401
      enddo
 
402
 20   continue
 
403
      filmap=.true.
 
404
      return
 
405
      end
 
406
 
 
407
 
 
408
      subroutine checkbw(nbw,ibwlist,isbw)
 
409
c**************************************************************************
 
410
c      Checks if any resonances are on the BW for this configuration
 
411
c**************************************************************************
 
412
      implicit none
 
413
      include 'genps.inc'
 
414
      include 'maxconfigs.inc'
 
415
      include 'nexternal.inc'
 
416
C $B$ NGRAPHS $E$ !this is a tag for MadWeight
 
417
 
 
418
      integer nbw,ibwlist(nexternal)
 
419
      logical isbw(*)
 
420
 
 
421
      logical             OnBW(-nexternal:0)     !Set if event is on B.W.
 
422
      common/to_BWEvents/ OnBW
 
423
      integer            mapconfig(0:lmaxconfigs), this_config
 
424
      common/to_mconfigs/mapconfig, this_config
 
425
C $B$ IFOREST $B$ !this is a tag for MadWeight
 
426
      integer i, iforest(2,-max_branch:-1,lmaxconfigs)
 
427
      common/to_forest/ iforest
 
428
C $E$ IFOREST $E$ !this is a tag for MadWeight
 
429
C $B$ DECAYBW $E$ !this is a tag for MadWeight
 
430
 
 
431
      integer icl(-(nexternal-3):nexternal)
 
432
      integer ibw
 
433
 
 
434
      nbw=0
 
435
      do i=-1,-(nexternal-3),-1
 
436
C $B$ ONBW $B$ !this is a tag for MadWeight
 
437
        if(OnBW(i)) then 
 
438
C $E$ ONBW $E$ !this is a tag for MadWeight
 
439
           nbw=nbw+1
 
440
        endif
 
441
      enddo
 
442
      if(nbw.eq.0)then
 
443
c        print *,'No BW found'
 
444
        return
 
445
      endif
 
446
 
 
447
      do i=1,nexternal
 
448
        icl(i)=ishft(1,i-1)
 
449
      enddo
 
450
      ibw=0
 
451
      do i=-1,-(nexternal-3),-1
 
452
        icl(i)=icl(iforest(1,i,this_config))+
 
453
     $     icl(iforest(2,i,this_config))
 
454
        isbw(icl(i))=.false.
 
455
C $B$ ONBW $B$ !this is a tag for MadWeight
 
456
        if(OnBW(i))then
 
457
C $E$ ONBW $E$ !this is a tag for MadWeight
 
458
          ibw=ibw+1
 
459
          ibwlist(ibw)=icl(i)
 
460
          isbw(icl(i))=.true.
 
461
c          print *,'Added BW for resonance ',i,icl(i),this_config
 
462
          if(ibw.eq.nbw) return
 
463
        endif
 
464
      enddo
 
465
      
 
466
      end
 
467
 
 
468
      logical function findmt(idij,icgs,nbw,ibwlist)
 
469
c**************************************************************************
 
470
c     input:
 
471
c            idij, icgs
 
472
c     output:
 
473
c            true if tree structure identified
 
474
c**************************************************************************
 
475
      implicit none
 
476
      include 'nexternal.inc'
 
477
      include 'maxamps.inc'
 
478
      include 'cluster.inc'
 
479
      include 'message.inc'
 
480
 
 
481
      integer idij,nbw,ibwlist(nexternal),icgs(0:n_max_cg)
 
482
      logical foundbw
 
483
      integer i, ii, j, jj, il, igsbk(0:n_max_cg)
 
484
 
 
485
c     IPROC has the present process number
 
486
      INTEGER IMIRROR,IPROC
 
487
      COMMON/TO_MIRROR/IMIRROR, IPROC
 
488
 
 
489
      findmt=.false.
 
490
c     if first clustering, set possible graphs
 
491
      if (icgs(0).eq.0) then
 
492
         ii=0
 
493
         do i=1,id_cl(iproc,idij,0)
 
494
c        check if we have constraint from onshell resonances
 
495
           foundbw=.true.
 
496
           do j=1,nbw
 
497
             if(resmap(ibwlist(j),id_cl(iproc,idij,i)))then
 
498
               cycle
 
499
             endif
 
500
             foundbw=.false.
 
501
 10        enddo
 
502
c        check if this diagram has radiation off heavy particle
 
503
c        For now, turn this off if there is a bw in the process,
 
504
c        since this created problems when the bw included a radiated part.
 
505
c        This is something that might need to be thought more about later
 
506
c        (in order for b matching to work...)
 
507
           if(nbw.eq.0.and.heavyrad(id_cl(iproc,idij,i))) cycle
 
508
           if((nbw.eq.0.or.foundbw))then
 
509
              ii=ii+1
 
510
              icgs(ii)=id_cl(iproc,idij,i)
 
511
           endif
 
512
         enddo
 
513
         icgs(0)=ii
 
514
         if (icgs(0).gt.0)then
 
515
           findmt=.true.
 
516
         endif
 
517
         if (btest(mlevel,5))
 
518
     $        write (*,*)'findmt: ',findmt,' IPROC=',IPROC,' icgs(0)=',icgs(0),
 
519
     $        ' icgs = ',(icgs(i),i=1,icgs(0))
 
520
         return
 
521
      else
 
522
c     Check for common graphs
 
523
         j=1
 
524
         ii=0
 
525
         do i=1,icgs(0)
 
526
            if(j.le.id_cl(iproc,idij,0).and.icgs(i).eq.id_cl(iproc,idij,j))then
 
527
               ii=ii+1
 
528
               icgs(ii)=id_cl(iproc,idij,j)
 
529
               j=j+1
 
530
            else if(j.le.id_cl(iproc,idij,0).and.icgs(i).gt.id_cl(iproc,idij,j)) then
 
531
               do while(icgs(i).gt.id_cl(iproc,idij,j).and.j.le.id_cl(iproc,idij,0))
 
532
                  j=j+1
 
533
               enddo
 
534
               if(j.le.id_cl(iproc,idij,0).and.icgs(i).eq.id_cl(iproc,idij,j))then
 
535
                  ii=ii+1
 
536
                  icgs(ii)=id_cl(iproc,idij,j)
 
537
               endif
 
538
            endif
 
539
         enddo
 
540
         icgs(0)=ii
 
541
         findmt=(icgs(0).gt.0)
 
542
         return
 
543
      endif
 
544
      end
 
545
 
 
546
 
 
547
      logical function cluster(p)
 
548
c**************************************************************************
 
549
c     input:
 
550
c            p(0:3,i)           momentum of i'th parton
 
551
c     output:
 
552
c            true if tree structure identified
 
553
c**************************************************************************
 
554
      implicit none
 
555
      include 'run.inc'
 
556
      include 'genps.inc'
 
557
      include 'nexternal.inc'
 
558
      include 'maxamps.inc'
 
559
      include 'cluster.inc'
 
560
      include 'message.inc'
 
561
      include 'maxconfigs.inc'
 
562
 
 
563
      real*8 p(0:3,nexternal), pcmsp(0:3), p1(0:3)
 
564
      real*8 pi(0:3), nr(0:3), pz(0:3)
 
565
      integer i, j, k, n, idi, idj, idij, icgs(0:n_max_cg)
 
566
      integer nleft, iwin, jwin, iwinp, imap(nexternal,2) 
 
567
      double precision nn2,ct,st
 
568
      double precision minpt2ij,pt2ij(n_max_cl),zij(n_max_cl)
 
569
 
 
570
      integer mapconfig(0:lmaxconfigs), this_config
 
571
      common/to_mconfigs/mapconfig, this_config
 
572
 
 
573
      integer nbw,ibwlist(nexternal)
 
574
      data isbw/n_max_cl*.false./
 
575
 
 
576
      data (pz(i),i=0,3)/1d0,0d0,0d0,1d0/
 
577
 
 
578
      integer combid
 
579
      logical findmt
 
580
      external findmt
 
581
      double precision dj, pydj, djb, pyjb, dot, SumDot, zclus
 
582
      external dj, pydj, djb, pyjb, dot, SumDot, zclus, combid
 
583
 
 
584
      if (btest(mlevel,1))
 
585
     $   write (*,*)'New event'
 
586
 
 
587
      cluster=.false.
 
588
      clustered=.false.
 
589
      do i=0,3
 
590
        pcmsp(i)=0
 
591
      enddo
 
592
c     Check if any resonances are on the BW, store results in to_checkbw
 
593
      call checkbw(nbw,ibwlist,isbw)
 
594
      if(btest(mlevel,4).and.nbw.gt.0)
 
595
     $     write(*,*) 'Found BWs: ',(ibwlist(i),i=1,nbw)
 
596
 
 
597
c     initialize index map
 
598
      do i=1,nexternal
 
599
         imap(i,1)=i
 
600
         imap(i,2)=ishft(1,i-1)
 
601
         mt2ij(i)=0
 
602
      enddo   
 
603
      mt2last=0
 
604
      minpt2ij=1.0d37
 
605
      do i=1,nexternal
 
606
c     initialize momenta
 
607
         idi=ishft(1,i-1)
 
608
         do j=0,3
 
609
            pcl(j,idi)=p(j,i)
 
610
         enddo
 
611
c     give mass to external particles
 
612
         pcl(4,idi)=dot(p(0,i),p(0,i))
 
613
c     never combine the two beams
 
614
         if (i.gt.2) then
 
615
c     fill combine table, first pass, determine all ptij
 
616
            do j=1,i-1
 
617
               idj=ishft(1,j-1)
 
618
               if (btest(mlevel,4))
 
619
     $              write (*,*)'i = ',i,'(',idi,'), j = ',j,'(',idj,')'
 
620
c     cluster only combinable legs (acc. to diagrams)
 
621
               icgs(0)=0
 
622
               idij=combid(idi,idj)
 
623
               pt2ij(idij)=1.0d37
 
624
               if (findmt(idij,icgs,nbw,ibwlist)) then
 
625
                  if (btest(mlevel,4)) then
 
626
                     write(*,*)'diagrams: ',(icgs(k),k=1,icgs(0))
 
627
                  endif
 
628
                  if (j.ne.1.and.j.ne.2) then
 
629
c     final state clustering                     
 
630
                     if(isbw(idij))then
 
631
                       pt2ij(idij)=SumDot(pcl(0,idi),pcl(0,idj),1d0)
 
632
                       if (btest(mlevel,4))
 
633
     $                    write(*,*)'Mother ',idij,' has ptij ',
 
634
     $                    sqrt(pt2ij(idij))
 
635
                     else
 
636
                       if(ktscheme.eq.2)then
 
637
                         pt2ij(idij)=pydj(pcl(0,idi),pcl(0,idj))
 
638
                       else
 
639
                         pt2ij(idij)=dj(pcl(0,idi),pcl(0,idj))
 
640
                       endif
 
641
                     endif
 
642
                     zij(idij)=0d0
 
643
                  else
 
644
c     initial state clustering, only if hadronic collision
 
645
c     check whether 2->(n-1) process w/ cms energy > 0 remains
 
646
                     iwinp=imap(3-j,2);
 
647
                     if(ickkw.eq.2.or.ktscheme.eq.2)then
 
648
                        pt2ij(idij)=pyjb(pcl(0,idi),
 
649
     $                    pcl(0,idj),pcl(0,iwinp),zij(idij))
 
650
                        zij(idij)=0d0
 
651
                     else
 
652
                        pt2ij(idij)=djb(pcl(0,idi))
 
653
                        zij(idij)=zclus(pcl(0,idi),pcl(0,idj),pcl(0,iwinp))
 
654
                     endif
 
655
c     prefer clustering when outgoing in direction of incoming
 
656
                     if(sign(1d0,pcl(3,idi)).ne.sign(1d0,pcl(3,idj)))
 
657
     $                    pt2ij(idij)=pt2ij(idij)*(1d0+1d-6)
 
658
                  endif
 
659
                  if (btest(mlevel,4)) then
 
660
                     write(*,*)'         ',idi,'&',idj,' part ',iwinp,
 
661
     &                              ' -> ',idij,' pt2ij = ',pt2ij(idij)
 
662
                     if(j.eq.1.or.j.eq.2)then
 
663
                       write(*,*)'     cf. djb: ',djb(pcl(0,idi))
 
664
                     endif
 
665
                  endif
 
666
c     Check if smallest pt2 ("winner")
 
667
                  if (pt2ij(idij).lt.minpt2ij) then
 
668
                     iwin=j
 
669
                     jwin=i
 
670
                     minpt2ij=pt2ij(idij)
 
671
                  endif                 
 
672
               endif
 
673
            enddo
 
674
         endif
 
675
      enddo
 
676
c     Take care of special 2 -> 1 case
 
677
      if (nexternal.eq.3.and.nincoming.eq.2) then
 
678
         n=1
 
679
c     Make sure that initial-state particles are daughters
 
680
         idacl(n,1)=imap(1,2)
 
681
         idacl(n,2)=imap(2,2)
 
682
         imocl(n)=imap(3,2)
 
683
         pt2ijcl(n)=pcl(4,imocl(n))
 
684
         zcl(n)=0.
 
685
         igraphs(0)=1
 
686
         igraphs(1)=this_config
 
687
         cluster=.true.
 
688
         clustered=.true.
 
689
         return
 
690
      endif
 
691
c     initialize graph storage
 
692
      igraphs(0)=0
 
693
      nleft=nexternal
 
694
c     cluster 
 
695
      do n=1,nexternal-2
 
696
c     combine winner
 
697
         imocl(n)=imap(iwin,2)+imap(jwin,2)
 
698
         idacl(n,1)=imap(iwin,2)
 
699
         idacl(n,2)=imap(jwin,2)
 
700
         pt2ijcl(n)=minpt2ij
 
701
         zcl(n)=zij(imocl(n))
 
702
         if (btest(mlevel,2)) then
 
703
            write(*,*)'winner ',n,': ',idacl(n,1),'&',idacl(n,2),
 
704
     &           ' -> ',minpt2ij,', z = ',zcl(n)
 
705
         endif
 
706
c     Reset igraphs with new mother
 
707
         if (.not.findmt(imocl(n),igraphs,nbw,ibwlist)) then
 
708
            write(*,*) 'cluster.f: Error. Invalid combination.' 
 
709
            return
 
710
         endif
 
711
         if (btest(mlevel,4)) then
 
712
            write(*,*)'graphs: ',(igraphs(k),k=1,igraphs(0))
 
713
         endif
 
714
         if (iwin.lt.3) then
 
715
c     is clustering
 
716
c     Set mt2ij to m^2+pt^2 
 
717
            mt2ij(n)=djb(pcl(0,idacl(n,2)))
 
718
            if (btest(mlevel,1)) then
 
719
               write(*,*)'mtij(',n,') for ',idacl(n,2),' is ',sqrt(mt2ij(n)),
 
720
     $              ' (cf ',sqrt(pt2ijcl(n)),')'
 
721
            endif
 
722
            iwinp=imap(3-iwin,2);
 
723
            do i=0,3
 
724
               pcl(i,imocl(n))=pcl(i,idacl(n,1))-pcl(i,idacl(n,2))
 
725
c            enddo
 
726
c     set incoming particle on-shell
 
727
c            pcl(0,imocl(n))=sqrt(pcl(1,imocl(n))**2+
 
728
c     $         pcl(2,imocl(n))**2+pcl(3,imocl(n))**2)
 
729
c            do i=0,3
 
730
               pcmsp(i)=-pcl(i,imocl(n))-pcl(i,iwinp)
 
731
            enddo
 
732
            pcmsp(0)=-pcmsp(0)
 
733
            pcl(4,imocl(n))=0
 
734
            if(pcl(4,idacl(n,1)).gt.0.or.pcl(4,idacl(n,2)).gt.0.and..not.
 
735
     $         (pcl(4,idacl(n,1)).gt.0.and.pcl(4,idacl(n,2)).gt.0))
 
736
     $         pcl(4,imocl(n))=max(pcl(4,idacl(n,1)),pcl(4,idacl(n,2)))
 
737
 
 
738
c       Don't boost if boost vector too lightlike or last vertex 
 
739
            if (pcmsp(0)**2-pcmsp(1)**2-pcmsp(2)**2-pcmsp(3)**2.gt.100d0.and.
 
740
     $           nleft.gt.4) then
 
741
               call boostx(pcl(0,imocl(n)),pcmsp(0),p1(0))
 
742
               call constr(p1(0),pz(0),nr(0),nn2,ct,st)
 
743
               do j=1,nleft
 
744
                  call boostx(pcl(0,imap(j,2)),pcmsp(0),p1(0))
 
745
                  call rotate(p1(0),pi(0),nr(0),nn2,ct,st,1)
 
746
                  do k=0,3
 
747
                     pcl(k,imap(j,2))=pi(k)
 
748
                  enddo
 
749
               enddo
 
750
               call boostx(pcl(0,imocl(n)),pcmsp(0),p1(0))
 
751
               call rotate(p1(0),pi(0),nr(0),nn2,ct,st,1)
 
752
               do k=0,3
 
753
                  pcl(k,imocl(n))=pi(k)
 
754
               enddo
 
755
            endif
 
756
         else
 
757
c     fs clustering
 
758
           do i=0,3
 
759
             pcl(i,imocl(n))=pcl(i,idacl(n,1))+pcl(i,idacl(n,2))
 
760
           enddo
 
761
           pcl(4,imocl(n))=0
 
762
           if(pcl(4,idacl(n,1)).gt.0.or.pcl(4,idacl(n,2)).gt.0.and..not.
 
763
     $        (pcl(4,idacl(n,1)).gt.0.and.pcl(4,idacl(n,2)).gt.0))
 
764
     $        pcl(4,imocl(n))=max(pcl(4,idacl(n,1)),pcl(4,idacl(n,2)))
 
765
           if(isbw(imocl(n)))then
 
766
             pcl(4,imocl(n))=pt2ijcl(n)
 
767
             if (btest(mlevel,4))
 
768
     $          write(*,*) 'Mother ',imocl(n),' has mass**2 ',
 
769
     $          pcl(4,imocl(n))
 
770
           endif
 
771
         endif
 
772
         
 
773
         nleft=nleft-1
 
774
c     create new imap
 
775
         imap(iwin,2)=imocl(n)
 
776
         do i=jwin,nleft
 
777
            imap(i,1)=imap(i+1,1)
 
778
            imap(i,2)=imap(i+1,2)
 
779
         enddo
 
780
         if (nleft.le.3) then
 
781
c           If last clustering is FS, store also average transverse mass
 
782
c           of the particles combined (for use if QCD vertex, e.g. tt~ or qq~)
 
783
            if(iwin.gt.2)then
 
784
               mt2last=sqrt(djb(pcl(0,idacl(n,1)))*djb(pcl(0,idacl(n,2))))
 
785
               if (btest(mlevel,3)) then
 
786
                  write(*,*)'Set mt2last to ',mt2last
 
787
               endif              
 
788
c         Boost and rotate back to get m_T for final particle
 
789
               if (pcmsp(0)**2-pcmsp(1)**2-pcmsp(2)**2-pcmsp(3)**2.gt.100d0) then
 
790
                  call rotate(pcl(0,imap(3,2)),p1(0),nr(0),nn2,ct,st,-1)
 
791
                  do k=1,3
 
792
                     pcmsp(k)=-pcmsp(k)
 
793
                  enddo
 
794
                  call boostx(p1(0),pcmsp(0),pcl(0,imap(3,2)))
 
795
               endif
 
796
            endif
 
797
c         Make sure that initial-state particle is always among daughters
 
798
            idacl(n+1,1)=imap(1,2)
 
799
            idacl(n+1,2)=imap(2,2)
 
800
            imocl(n+1)=imap(3,2)
 
801
 
 
802
c            if(pcl(0,imocl(n)).gt.0d0)then
 
803
            pt2ijcl(n+1)=djb(pcl(0,imap(3,2)))
 
804
            if (btest(mlevel,3)) then
 
805
              write(*,*) 'Last vertex is ',imap(1,2),imap(2,2),imap(3,2)
 
806
              write(*,*) '            -> ',pt2ijcl(n+1),sqrt(pt2ijcl(n+1))
 
807
            endif
 
808
c     If present channel among graphs, use only this channel
 
809
c     This is important when we have mixed QED-QCD
 
810
            do i=1,igraphs(0)
 
811
               if (igraphs(i).eq.this_config) then
 
812
                  igraphs(0)=1
 
813
                  igraphs(1)=this_config
 
814
                  exit
 
815
               endif
 
816
            enddo
 
817
c            if(pt2ijcl(n).gt. pt2ijcl(n+1))then
 
818
c              pt2ijcl(n+1)=pt2ijcl(n)
 
819
c              if (btest(mlevel,3)) then
 
820
c                write(*,*)'Reset scale for vertex ',n+1,' to ',pt2ijcl(n+1)
 
821
c              endif              
 
822
c            endif
 
823
            zcl(n+1)=1
 
824
c            else
 
825
c              pt2ijcl(n+1)=pt2ijcl(n)
 
826
c            endif
 
827
c           Pick out the found graphs
 
828
c            print *,'Clustering succeeded, found graph ',igscl(1)
 
829
            cluster=.true.
 
830
            clustered=.true.
 
831
            return
 
832
         endif
 
833
c     calculate new ptij
 
834
c            write(*,*)'is case'
 
835
c     recalculate all in is case due to rotation & boost
 
836
         minpt2ij=1.0d37
 
837
            do i=1,nleft
 
838
               idi=imap(i,2)
 
839
c     never combine the two beams
 
840
               if (i.gt.2) then
 
841
c     determine all ptij
 
842
                  do j=1,i-1
 
843
                     idj=imap(j,2)
 
844
                     if (btest(mlevel,4))
 
845
     $                    write (*,*)'i = ',i,'(',idi,'), j = ',j,'(',idj,')'
 
846
c     Reset diagram list icgs
 
847
                     do k=0,igraphs(0)
 
848
                        icgs(k)=igraphs(k)
 
849
                     enddo
 
850
                     if (btest(mlevel,4))
 
851
     $                    write (*,*)'Reset diagrams to: ',(icgs(k),k=1,icgs(0))
 
852
c     cluster only combinable legs (acc. to diagrams)
 
853
                     idij=combid(idi,idj)
 
854
c                     write (*,*) 'RECALC !!! ',idij
 
855
                     pt2ij(idij)=1.0d37
 
856
                     if (findmt(idij,icgs,nbw,ibwlist)) then
 
857
                        if (btest(mlevel,4)) then
 
858
                           write(*,*)'diagrams: ',(icgs(k),k=1,icgs(0))
 
859
                       endif
 
860
                        if (j.ne.1.and.j.ne.2) then
 
861
c     final state clustering                     
 
862
                           if(isbw(idij))then
 
863
                             pt2ij(idij)=SumDot(pcl(0,idi),pcl(0,idj),1d0)
 
864
                             if (btest(mlevel,4))
 
865
     $                          write(*,*) 'Mother ',idij,' has ptij ',
 
866
     $                          sqrt(pt2ij(idij))
 
867
                           else
 
868
                             if(ktscheme.eq.2)then
 
869
                               pt2ij(idij)=pydj(pcl(0,idi),pcl(0,idj))
 
870
                             else
 
871
                               pt2ij(idij)=dj(pcl(0,idi),pcl(0,idj))
 
872
                             endif
 
873
                           endif
 
874
                           zij(idij)=0d0
 
875
                        else
 
876
c     initial state clustering, only if hadronic collision
 
877
c     check whether 2->(n-1) process w/ cms energy > 0 remains
 
878
                          iwinp=imap(3-j,2);
 
879
                          do k=0,3
 
880
                             pcl(k,idij)=pcl(k,idj)-pcl(k,idi)
 
881
c                           pcmsp(k)=pcl(k,idij)+pcl(k,iwinp)
 
882
                          enddo
 
883
c                       ecms2=pcmsp(0)**2-pcmsp(1)**2-
 
884
c                       $                          pcmsp(2)**2-pcmsp(3)**2
 
885
c                       if (ecms2.gt.0.1d0.and.
 
886
c                       if ((nleft.eq.4.or.ecms2.gt.0.1d0).and.
 
887
c                         if((lpp(j).ne.0)) then
 
888
                            if(ickkw.eq.2.or.ktscheme.eq.2)then
 
889
                              pt2ij(idij)=pyjb(pcl(0,idi),
 
890
     $                           pcl(0,idj),pcl(0,iwinp),zij(idij))
 
891
                            else
 
892
                              pt2ij(idij)=djb(pcl(0,idi))
 
893
                              zij(idij)=zclus(pcl(0,idi),pcl(0,idj),pcl(0,iwinp))
 
894
                            endif
 
895
c                 prefer clustering when outgoing in direction of incoming
 
896
                            if(sign(1d0,pcl(3,idi)).ne.sign(1d0,pcl(3,idj)))
 
897
     $                         pt2ij(idij)=pt2ij(idij)*(1d0+1d-6)
 
898
c                          endif
 
899
                        endif
 
900
                        if (btest(mlevel,4)) then
 
901
                          write(*,*)'         ',idi,'&',idj,' part ',iwinp,' -> ',idij,
 
902
     &                       ' pt2ij = ',pt2ij(idij)
 
903
                           if(j.eq.1.or.j.eq.2)then
 
904
                             write(*,*)'     cf. djb: ',djb(pcl(0,idi))
 
905
                           endif
 
906
                        endif
 
907
                        if (pt2ij(idij).lt.minpt2ij) then
 
908
                           iwin=j
 
909
                           jwin=i
 
910
                           minpt2ij=pt2ij(idij)
 
911
                        endif                 
 
912
                     endif
 
913
                  enddo
 
914
               endif
 
915
            enddo
 
916
      enddo
 
917
 
 
918
      return
 
919
      end