← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.1

« back to all changes in this revision

Viewing changes to Template/NLO/SubProcesses/cluster.f

  • Committer: Marco Zaro
  • Date: 2014-01-27 16:54:10 UTC
  • mfrom: (78.124.55 MG5_aMC_2.1)
  • Revision ID: marco.zaro@gmail.com-20140127165410-5lma8c2hzbzm426j
merged with lp:~maddevelopers/madgraph5/MG5_aMC_2.1 r 267

Show diffs side-by-side

added added

removed removed

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