← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2-alpha0

« back to all changes in this revision

Viewing changes to Template/SubProcesses/myamp.f

Added Template and HELAS into bzr

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/SubProcesses/myamp.f
 
1
      double precision function testamp(p)
 
2
c*****************************************************************************
 
3
c     Approximates matrix element by propagators
 
4
c*****************************************************************************
 
5
      implicit none
 
6
c
 
7
c     Constants
 
8
c     
 
9
      include 'genps.inc'
 
10
      include 'nexternal.inc'
 
11
      double precision   zero
 
12
      parameter (zero = 0d0)
 
13
c
 
14
c     Arguments
 
15
c
 
16
      double precision p(0:3,nexternal)
 
17
c      integer iconfig
 
18
c
 
19
c     Local
 
20
c
 
21
      double precision xp(0:3,-nexternal:nexternal)
 
22
      double precision mpole(-nexternal:0),shat,tsgn
 
23
      integer i,j,iconfig
 
24
 
 
25
      double precision pmass(-nexternal:0,lmaxconfigs)
 
26
      double precision pwidth(-nexternal:0,lmaxconfigs)
 
27
      integer pow(-nexternal:0,lmaxconfigs)
 
28
      logical first_time
 
29
c
 
30
c     Global
 
31
c
 
32
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
33
      common/to_forest/ iforest
 
34
      integer            mapconfig(0:lmaxconfigs), this_config
 
35
      common/to_mconfigs/mapconfig, this_config
 
36
      
 
37
      include 'coupl.inc'
 
38
c
 
39
c     External
 
40
c
 
41
      double precision dot
 
42
 
 
43
      save pmass,pwidth,pow
 
44
      data first_time /.true./
 
45
c-----
 
46
c  Begin Code
 
47
c-----      
 
48
      iconfig = this_config
 
49
      if (first_time) then
 
50
c         include 'props.inc'
 
51
         first_time=.false.
 
52
      endif
 
53
 
 
54
      do i=1,nexternal
 
55
         mpole(-i)=0d0
 
56
         do j=0,3
 
57
            xp(j,i)=p(j,i)
 
58
         enddo
 
59
      enddo
 
60
c      mpole(-3) = 174**2
 
61
c      shat = dot(p(0,1),p(0,2))/(1800)**2
 
62
      shat = dot(p(0,1),p(0,2))/(10)**2
 
63
c      shat = 1d0
 
64
      testamp = 1d0
 
65
      tsgn    = +1d0
 
66
      do i=-1,-(nexternal-3),-1              !Find all the propagotors
 
67
         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
68
         do j=0,3
 
69
            xp(j,i) = xp(j,iforest(1,i,iconfig))
 
70
     $           +tsgn*xp(j,iforest(2,i,iconfig))
 
71
         enddo
 
72
         if (pwidth(i,iconfig) .ne. 0d0 .and. .false.) then
 
73
            testamp=testamp/((dot(xp(0,i),xp(0,i))
 
74
     $                        -pmass(i,iconfig)**2)**2
 
75
     $         -(pmass(i,iconfig)*pwidth(i,iconfig))**2)
 
76
         else
 
77
            testamp = testamp/((dot(xp(0,i),xp(0,i)) -
 
78
     $                          pmass(i,iconfig)**2)
 
79
     $                          **(pow(i,iconfig)))
 
80
         endif
 
81
        testamp=testamp*shat**(pow(i,iconfig))
 
82
c        write(*,*) i,iconfig,pow(i,iconfig),pmass(i,iconfig)
 
83
      enddo
 
84
c      testamp = 1d0/dot(xp(0,-1),xp(0,-1))
 
85
      testamp=abs(testamp)
 
86
c      testamp = 1d0
 
87
      end
 
88
 
 
89
      logical function cut_bw(p)
 
90
c*****************************************************************************
 
91
c     Approximates matrix element by propagators
 
92
c*****************************************************************************
 
93
      implicit none
 
94
c
 
95
c     Constants
 
96
c     
 
97
      include 'genps.inc'
 
98
      include 'nexternal.inc'
 
99
      double precision   zero
 
100
      parameter (zero = 0d0)
 
101
      include 'run.inc'
 
102
c
 
103
c     Arguments
 
104
c
 
105
      double precision p(0:3,nexternal)
 
106
c
 
107
c     Local
 
108
c
 
109
      double precision xp(0:3,-nexternal:nexternal)
 
110
      double precision mpole(-nexternal:0),shat,tsgn
 
111
      integer i,j,iconfig
 
112
 
 
113
      double precision pmass(-nexternal:0,lmaxconfigs)
 
114
      double precision pwidth(-nexternal:0,lmaxconfigs)
 
115
      integer pow(-nexternal:0,lmaxconfigs)
 
116
      logical first_time, onshell
 
117
      double precision xmass
 
118
      integer nbw
 
119
 
 
120
      integer ida(2),idenpart
 
121
c
 
122
c     Global
 
123
c
 
124
      include 'maxamps.inc'
 
125
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
126
      common/to_forest/ iforest
 
127
      integer sprop(-max_branch:-1,lmaxconfigs)
 
128
      integer tprid(-max_branch:-1,lmaxconfigs)
 
129
      common/to_sprop/sprop,tprid
 
130
      integer            mapconfig(0:lmaxconfigs), this_config
 
131
      common/to_mconfigs/mapconfig, this_config
 
132
 
 
133
      integer        lbw(0:nexternal)  !Use of B.W.
 
134
      common /to_BW/ lbw
 
135
 
 
136
      logical             OnBW(-nexternal:0)     !Set if event is on B.W.
 
137
      common/to_BWEvents/ OnBW
 
138
      
 
139
      include 'coupl.inc'
 
140
 
 
141
      integer idup(nexternal,maxproc)
 
142
      integer mothup(2,nexternal,maxproc)
 
143
      integer icolup(2,nexternal,maxflow)
 
144
      include 'leshouche.inc'
 
145
 
 
146
      logical gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
 
147
      include 'decayBW.inc'
 
148
c
 
149
c     External
 
150
c
 
151
      double precision dot
 
152
 
 
153
      save pmass,pwidth,pow
 
154
      data first_time /.true./
 
155
c-----
 
156
c  Begin Code
 
157
c-----      
 
158
      cut_bw = .false.    !Default is we passed the cut
 
159
      iconfig = this_config
 
160
      if (first_time) then
 
161
         include 'props.inc'
 
162
         nbw = 0
 
163
         tsgn = 1d0
 
164
         do i=-1,-(nexternal-3),-1
 
165
            if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
166
            if (pwidth(i,iconfig) .gt. 0 .and. tsgn .eq. 1d0) then
 
167
               nbw=nbw+1
 
168
               if (lbw(nbw) .eq. 1) then
 
169
                  write(*,*) 'Requiring BW ',i,nbw
 
170
               elseif(lbw(nbw) .eq. 2) then
 
171
                  write(*,*) 'Excluding BW ',i,nbw
 
172
               else
 
173
                  write(*,*) 'No cut BW ',i,nbw
 
174
               endif
 
175
            endif
 
176
         enddo
 
177
         first_time=.false.
 
178
      endif
 
179
 
 
180
      do i=1,nexternal
 
181
         mpole(-i)=0d0
 
182
         do j=0,3
 
183
            xp(j,i)=p(j,i)
 
184
         enddo
 
185
      enddo
 
186
      nbw = 0
 
187
      tsgn    = +1d0
 
188
      do i=-1,-(nexternal-3),-1              !Loop over propagators
 
189
         onbw(i) = .false.
 
190
         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
191
         do j=0,3
 
192
            xp(j,i) = xp(j,iforest(1,i,iconfig))
 
193
     $           +tsgn*xp(j,iforest(2,i,iconfig))
 
194
         enddo
 
195
         if (tsgn .gt. 0d0 .and. pwidth(i,iconfig) .gt. 0d0 ) then !This is B.W.
 
196
            nbw = nbw+1
 
197
c            write(*,*) 'Checking BW',nbw
 
198
            xmass = sqrt(dot(xp(0,i),xp(0,i)))
 
199
c            write(*,*) 'xmass',xmass,pmass(i,iconfig)
 
200
            onshell = (abs(xmass - pmass(i,iconfig)) .lt.
 
201
     $           bwcutoff*pwidth(i,iconfig))
 
202
 
 
203
c
 
204
c           Here we set if the BW is "on-shell" for LesHouches
 
205
c
 
206
            if(onshell)then
 
207
c           Only allow onshell if no "decay" to identical particle
 
208
              OnBW(i) = .true.
 
209
              idenpart=0
 
210
              do j=1,2
 
211
                ida(j)=iforest(j,i,iconfig)
 
212
                if((ida(j).lt.0.and.sprop(i,iconfig).eq.sprop(ida(j),iconfig))
 
213
     $             .or.(ida(j).gt.0.and.sprop(i,iconfig).eq.IDUP(ida(j),1)))
 
214
     $             idenpart=ida(j)
 
215
              enddo
 
216
c           Always remove if daughter final-state
 
217
              if(idenpart.gt.0) then
 
218
                OnBW(i)=.false.
 
219
c           Else remove if daughter forced to be onshell
 
220
              elseif(idenpart.lt.0.and.gForceBW(idenpart, iconfig)) then
 
221
                OnBW(i)=.false.
 
222
c           Else remove daughter if forced to be onshell
 
223
              elseif(idenpart.lt.0.and.gForceBW(i, iconfig)) then
 
224
                OnBW(idenpart)=.false.
 
225
c           Else remove either this resonance or daughter, which is closer to mass shell
 
226
              elseif(idenpart.lt.0.and.abs(xmass-pmass(i,iconfig)).gt.
 
227
     $             abs(sqrt(dot(xp(0,idenpart),xp(0,idenpart)))-
 
228
     $             pmass(i,iconfig))) then
 
229
                OnBW(i)=.false.
 
230
c           Else remove OnBW for daughter
 
231
              else if(idenpart.lt.0) then
 
232
                OnBW(idenpart)=.false.
 
233
              endif
 
234
            endif
 
235
            if (onshell .and. (lbw(nbw).eq. 2) ) cut_bw=.true.
 
236
            if (.not. onshell .and. (lbw(nbw).eq. 1)) cut_bw=.true.
 
237
c            write(*,*) nbw,xmass,onshell,lbw(nbw),cut_bw
 
238
         endif
 
239
 
 
240
      enddo
 
241
      end
 
242
 
 
243
 
 
244
      subroutine set_peaks
 
245
c*****************************************************************************
 
246
c     Attempts to determine peaks for this configuration
 
247
c*****************************************************************************
 
248
      implicit none
 
249
c
 
250
c     Constants
 
251
c     
 
252
      include 'genps.inc'
 
253
      include 'nexternal.inc'
 
254
      double precision   zero
 
255
      parameter (zero = 0d0)
 
256
c
 
257
c     Arguments
 
258
c
 
259
c
 
260
c     Local
 
261
c
 
262
      double precision  xm(-nexternal:nexternal)
 
263
      double precision  xe(-nexternal:nexternal)
 
264
      double precision tsgn, xo, a
 
265
      double precision x1,x2,xk(nexternal)
 
266
      double precision dr,mtot,etot,stot,xqfact
 
267
      integer i, iconfig, l1, l2, j, nt, nbw
 
268
 
 
269
      double precision pmass(-nexternal:0,lmaxconfigs)
 
270
      double precision pwidth(-nexternal:0,lmaxconfigs)
 
271
      integer pow(-nexternal:0,lmaxconfigs)
 
272
 
 
273
c
 
274
c     Global
 
275
c
 
276
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
277
      common/to_forest/ iforest
 
278
 
 
279
      integer            mapconfig(0:lmaxconfigs), this_config
 
280
      common/to_mconfigs/mapconfig, this_config
 
281
 
 
282
      real*8         emass(nexternal)
 
283
      common/to_mass/emass
 
284
 
 
285
      include 'run.inc'
 
286
 
 
287
      double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
 
288
      double precision emin(nincoming+1:nexternal)
 
289
      double precision                    r2min(nincoming+1:nexternal,nincoming+1:nexternal)
 
290
      double precision s_min(nexternal,nexternal)
 
291
      common/to_cuts/  etmin, emin, etamax, r2min, s_min
 
292
 
 
293
      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
 
294
      common/to_xqcuts/xqcutij,xqcuti
 
295
 
 
296
      double precision      spole(maxinvar),swidth(maxinvar),bwjac
 
297
      common/to_brietwigner/spole          ,swidth          ,bwjac
 
298
 
 
299
      integer        lbw(0:nexternal)  !Use of B.W.
 
300
      common /to_BW/ lbw
 
301
      
 
302
      include 'coupl.inc'
 
303
 
 
304
c
 
305
c     External
 
306
c
 
307
 
 
308
c-----
 
309
c  Begin Code
 
310
c-----      
 
311
      include 'props.inc'
 
312
      stot = 4d0*ebeam(1)*ebeam(2)
 
313
c      etmin = 10
 
314
      nt = 0
 
315
      iconfig = this_config
 
316
      mtot = 0d0
 
317
      etot = 0d0   !Total energy needed
 
318
      xqfact=1d0
 
319
      if(ickkw.eq.2.or.ktscheme.eq.2) xqfact=0.3d0
 
320
      do i=nincoming+1,nexternal  !assumes 2 incoming
 
321
         xm(i)=emass(i)
 
322
c-fax
 
323
         xe(i)=max(emass(i),max(etmin(i),0d0))
 
324
         xe(i)=max(xe(i),max(emin(i),0d0))
 
325
c-JA 1/2009: Set grid also based on xqcut
 
326
         xe(i)=max(xe(i),xqfact*xqcuti(i))
 
327
         xk(i)= 0d0
 
328
         etot = etot+xe(i)
 
329
         mtot=mtot+xm(i)         
 
330
      enddo
 
331
      tsgn    = +1d0
 
332
      nbw = 0
 
333
      do i=-1,-(nexternal-3),-1              !Find all the propagotors
 
334
         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
335
         if (tsgn .eq. 1d0) then                         !s channel
 
336
            xm(i) = xm(iforest(1,i,iconfig))+xm(iforest(2,i,iconfig))
 
337
            xe(i) = xe(iforest(1,i,iconfig))+xe(iforest(2,i,iconfig))
 
338
            mtot = mtot - xm(i)
 
339
            etot = etot - xe(i)
 
340
            if (iforest(1,i,iconfig) .gt. 0
 
341
     &           .and. iforest(2,i,iconfig) .gt. 0) then
 
342
c-JA 1/2009: Set deltaR cuts here together with s_min cuts
 
343
              l1 = iforest(1,i,iconfig)
 
344
              l2 = iforest(2,i,iconfig)
 
345
              xm(i)=max(xm(i),sqrt(max(s_min(l1,l2),0d0)))
 
346
              dr = max(r2min(l1,l2)*dabs(r2min(l1,l2)),0d0)*0.8d0
 
347
              xm(i)=max(xm(i),
 
348
     &           sqrt(max(etmin(l2),0d0)*max(etmin(l1),0d0)*dr))
 
349
c-JA 1/2009: Set grid also based on xqcut
 
350
              xm(i)=max(xm(i),sqrt(max(xqcutij(l1,l2),0d0)))
 
351
              xe(i)=max(xe(i),xm(i))
 
352
            endif
 
353
c            write(*,*) 'iconfig,i',iconfig,i
 
354
c            write(*,*) pwidth(i,iconfig),pmass(i,iconfig)
 
355
            if (pwidth(i,iconfig) .gt. 0) nbw=nbw+1
 
356
            if (pwidth(i,iconfig) .gt. 0 .and. lbw(nbw) .le. 1) then         !B.W.
 
357
c               nbw = nbw +1
 
358
 
 
359
               if (i .eq. -(nexternal-(nincoming+1))) then  !This is s-hat
 
360
                  j = 3*(nexternal-2)-4+1    !set i to ndim+1
 
361
c-----
 
362
c tjs 11/2008 if require BW then force even if worried about energy
 
363
c----
 
364
                  if(pmass(i,iconfig).ge.xe(i) .or. lbw(nbw).eq.1) then
 
365
                     write(*,*) 'Setting PDF BW',j,nbw,pmass(i,iconfig)
 
366
                     spole(j)=pmass(i,iconfig)*pmass(i,iconfig)/stot
 
367
                     swidth(j) = pwidth(i,iconfig)*pmass(i,iconfig)/stot
 
368
                     xm(i) = pmass(i,iconfig)
 
369
                  else
 
370
                     spole(j)=0d0
 
371
                     swidth(j) = 0d0
 
372
                  endif
 
373
               else
 
374
                  write(*,*) 'Setting BW',i,nbw,pmass(i,iconfig)
 
375
                  spole(-i)=pmass(i,iconfig)*pmass(i,iconfig)/stot
 
376
                  swidth(-i) = pwidth(i,iconfig)*pmass(i,iconfig)/stot
 
377
                  xm(i) = pmass(i,iconfig)
 
378
c RF & TJS, should start from final state particle masses, not only at resonance.
 
379
c Therefore remove the next line.
 
380
c                  xe(i) = max(xe(i),xm(i))
 
381
               endif
 
382
            else                                  !1/x^pow
 
383
c-JA 1/2009: Comment out this whole section, since it only sets (wrong) xm(i)
 
384
c               if (xm(i) - pmass(i,iconfig) .le. 0d0) then !Can hit pole
 
385
cc                  write(*,*) 'Setting new min',i,xm(i),pmass(i,iconfig)
 
386
c                  l1 = iforest(1,i,iconfig)                  !need dr cut
 
387
c                  l2 = iforest(2,i,iconfig)
 
388
c                  if (l2 .lt. l1) then
 
389
c                     j = l1
 
390
c                     l1 = l2
 
391
c                     l2 = j
 
392
c                  endif
 
393
c                  dr = 0
 
394
cc-fax
 
395
c                  if (l1 .gt. 0) 
 
396
c     &  dr = max(r2min(l2,l1)*dabs(r2min(l2,l1)),0d0) !dr only for external
 
397
cc                  write(*,*) 'using r2min',l2,l1,sqrt(dr)
 
398
c                  dr = dr*.8d0                        !0.8 to hit peak hard
 
399
c                  xo = 0.5d0*pmass(i,iconfig)**2      !0.5 to hit peak hard
 
400
cc-fax
 
401
c                  if (dr .gt. 0d0) 
 
402
c     &            xo = max(xo,max(etmin(l2),0d0)*max(etmin(l1),0d0)*dr)
 
403
cc                  write(*,*) 'Got dr',i,l1,l2,dr
 
404
cc------
 
405
cctjs 11/2008  if explicitly missing pole, don't want to include mass in xm
 
406
cc-----
 
407
c                 if (pwidth(i,iconfig) .le. 0) then
 
408
cc                     write(*,*) "Including mass",i,pmass(i,iconfig)
 
409
c                    xo = xo+pmass(i,iconfig)**2
 
410
c                 else
 
411
cc                     write(*,*) "Skipping mass", i,pmass(i,iconfig),sqrt(xo)
 
412
c                 endif
 
413
cc                  write(*,*) 'Setting xm',i,xm(i),sqrt(xo)
 
414
c                  xm(i) = sqrt(xo)    !Reset xm to realistic minimum
 
415
c                  xo = xo/stot
 
416
cc                  xo = sqrt(pmass(i,iconfig)**2+ pmass(i,iconfig)**2)
 
417
cc-fax
 
418
cc                  xo = pmass(i,iconfig)+max(etmin,0d0)
 
419
c               else
 
420
c                  write(*,*) 'Using xm',i,xm(i)
 
421
c                  xo = xm(i)**2/stot
 
422
c               endif
 
423
              xo = xm(i)**2/stot
 
424
              a=pmass(i,iconfig)**2/stot
 
425
c               call setgrid(-i,xo,a,pow(i,iconfig))
 
426
c               write(*,*) 'Enter minimum for ',-i, xo
 
427
c               read(*,*) xo
 
428
 
 
429
 
 
430
               if (pwidth(i,iconfig) .eq. 0) call setgrid(-i,xo,a,1)
 
431
               if (pwidth(i,iconfig) .gt. 0) then
 
432
                  write(*,*) 'Using flat grid for BW',i,nbw,
 
433
     $                 pmass(i,iconfig)
 
434
               endif
 
435
            endif
 
436
c            xe(i) = max(xm(i),xe(i))               
 
437
            etot = etot+xe(i)
 
438
            mtot=mtot+max(xm(i),xm(i))
 
439
c            write(*,*) 'New mtot',i,mtot,xm(i)
 
440
         else                                        !t channel
 
441
c
 
442
c     Check closest to p1
 
443
c
 
444
            nt = nt+1
 
445
            l2 = iforest(2,i,iconfig) !need dr cut
 
446
            x1 = 0            
 
447
c-fax
 
448
c-JA 1/2009: Set grid also based on xqcut
 
449
            if (l2 .gt. 0) x1 = max(etmin(l2),max(xqfact*xqcuti(l2),0d0))
 
450
            x1 = max(x1, xm(l2)/1d0)
 
451
            if (nt .gt. 1) x1 = max(x1,xk(nt-1))
 
452
            xk(nt)=x1
 
453
c            write(*,*) 'Using 1',l2,x1
 
454
 
 
455
c
 
456
c     Check closest to p2
 
457
c
 
458
            j = i-1
 
459
            l2 = iforest(2,j,iconfig)
 
460
            x2 = 0
 
461
c-JA 1/2009: Set grid also based on xqcut
 
462
            if (l2 .gt. 0) x2 = max(etmin(l2),max(xqfact*xqcuti(l2),0d0))
 
463
c            if (l2 .gt. 0) x2 = max(etmin(l2),0d0)
 
464
            x2 = max(x2, xm(l2)/1d0)
 
465
c            if (nt .gt. 1) x2 = max(x2,xk(nt-1))
 
466
            
 
467
c            write(*,*) 'Using 2',l2,x2
 
468
 
 
469
            xo = min(x1,x2)
 
470
 
 
471
            xo = xo*xo/stot
 
472
            a=-pmass(i,iconfig)**2/stot
 
473
c            call setgrid(-i,xo,a,pow(i,iconfig))
 
474
 
 
475
c               write(*,*) 'Enter minimum for ',-i, xo
 
476
c               read(*,*) xo
 
477
 
 
478
             if (i .ne. -1 .or. .true.) call setgrid(-i,xo,a,1)
 
479
         endif
 
480
      enddo
 
481
      if (abs(lpp(1)) .eq. 1 .or. abs(lpp(2)) .eq. 1) then
 
482
         write(*,*) 'etot',etot,nexternal
 
483
         xo = etot**2/stot
 
484
         i = 3*(nexternal-2) - 4 + 1
 
485
c-----------------------
 
486
c     tjs  4/29/2008 use analytic transform for s-hat
 
487
c-----------------------
 
488
         if (swidth(i) .eq. 0d0) then
 
489
            swidth(i) = xo
 
490
            spole(i)= -2.0d0    ! 1/s pole
 
491
            write(*,*) "Transforming s_hat 1/s ",i,xo
 
492
         else
 
493
            write(*,*) "Transforming s_hat BW ",spole(i),swidth(i)
 
494
         endif
 
495
c-----------------------
 
496
         if (swidth(i) .eq. 0d0) then
 
497
            call setgrid(i,xo,0d0,1)
 
498
         endif
 
499
      endif
 
500
 
 
501
      i=-8
 
502
c      write(*,*) 'Enter minimum for ',-i, xo
 
503
c      read(*,*) xo      
 
504
c      if (xo .gt. 0)      call setgrid(-i,xo,a,1)
 
505
 
 
506
      i=-10
 
507
c      write(*,*) 'Enter minimum for ',-i, xo
 
508
c      read(*,*) xo      
 
509
c      if (xo .gt. 0) call setgrid(-i,xo,a,1)
 
510
 
 
511
      end
 
512
 
 
513
 
 
514
 
 
515
c      subroutine find_matches(iconfig,isym)
 
516
cc*****************************************************************************
 
517
cc     Finds all of the matches to see what gains may be possible
 
518
cc*****************************************************************************
 
519
c      implicit none
 
520
cc
 
521
cc     Constants
 
522
cc     
 
523
c      include 'genps.inc'
 
524
c      double precision   zero
 
525
c      parameter (zero = 0d0)
 
526
cc
 
527
cc     Arguments
 
528
cc
 
529
c      integer iconfig,isym(0:10)
 
530
cc
 
531
cc     Local
 
532
cc
 
533
c      integer i,j, nc, imatch
 
534
c      double precision tprop(3,nexternal),xprop(3,nexternal)
 
535
cc
 
536
cc     Global
 
537
cc
 
538
c      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
539
c      common/to_forest/ iforest
 
540
c
 
541
c      integer            mapconfig(0:lmaxconfigs), this_config
 
542
c      common/to_mconfigs/mapconfig, this_config
 
543
c
 
544
cc      include 'coupl.inc'
 
545
cc      include 'configs.inc'
 
546
cc
 
547
cc     External
 
548
cc
 
549
c
 
550
cc-----
 
551
cc  Begin Code
 
552
cc-----      
 
553
cc
 
554
cc     First get peaks for configuration interested in
 
555
cc
 
556
c      if (.false.) then
 
557
c         isym(0)=1
 
558
c         isym(1)=iconfig
 
559
c      else
 
560
c         call get_peaks(iconfig,tprop)
 
561
c         isym(0)=0
 
562
c         do i=1,mapconfig(0)
 
563
c            call get_peaks(i,xprop)
 
564
cc        call sort_prop(xprop)
 
565
c            call match_peak(tprop,xprop,imatch)
 
566
c            if (imatch .eq. 1) then
 
567
c               write(*,*) 'Found Match',iconfig,i
 
568
c               isym(0)=isym(0)+1
 
569
c               isym(isym(0)) = i
 
570
c            endif
 
571
c         enddo
 
572
c      endif
 
573
c      write(*,'(20i4)') (isym(i),i=0,isym(0))
 
574
c      end
 
575
c
 
576
c      subroutine find_all_matches()
 
577
cc*****************************************************************************
 
578
cc     Finds all of the matches to see what gains may be possible
 
579
cc*****************************************************************************
 
580
c      implicit none
 
581
cc
 
582
cc     Constants
 
583
cc     
 
584
c      include 'genps.inc'
 
585
c      double precision   zero
 
586
c      parameter (zero = 0d0)
 
587
cc
 
588
cc     Arguments
 
589
cc
 
590
c      double precision tprop(3,nexternal),xprop(3,nexternal)
 
591
c      integer iconfig
 
592
cc
 
593
cc     Local
 
594
cc
 
595
c      integer i,j, nc, imatch
 
596
c      logical gm(lmaxconfigs)
 
597
cc
 
598
cc     Global
 
599
cc
 
600
c      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
601
c      common/to_forest/ iforest
 
602
c
 
603
c      integer            mapconfig(0:lmaxconfigs), this_config
 
604
c      common/to_mconfigs/mapconfig, this_config
 
605
c
 
606
cc      include 'coupl.inc'
 
607
cc      include 'configs.inc'
 
608
cc
 
609
cc     External
 
610
cc
 
611
c
 
612
cc-----
 
613
cc  Begin Code
 
614
cc-----      
 
615
c      nc = 0
 
616
c      do i=1,mapconfig(0)
 
617
c         gm(i)=.false.
 
618
c      enddo
 
619
c      do j=1,mapconfig(0)
 
620
c         if (.not. gm(j)) then
 
621
c            nc=nc+1
 
622
cc            write(*,*) 'Need config ',j
 
623
c            call get_peaks(j,tprop)
 
624
cc            call sort_prop(tprop)
 
625
c            write(*,'(i4,4e12.4)') j,(tprop(1,i), i=1,4)
 
626
c            do i=j+1,mapconfig(0)
 
627
c               call get_peaks(i,xprop)
 
628
cc               call sort_prop(xprop)
 
629
c               call match_peak(tprop,xprop,imatch)
 
630
c               if (imatch .eq. 1) then
 
631
c                  write(*,*) 'Found Match',j,i
 
632
c                  gm(i)=.true.
 
633
c               endif
 
634
c            enddo
 
635
c         endif
 
636
c      enddo
 
637
c      write(*,*) 'Found matches',mapconfig(0),nc
 
638
c      stop
 
639
c      end
 
640
 
 
641
c      subroutine sort_prop(xprop)
 
642
cc*****************************************************************************
 
643
cc     Sort props in order from min to max based on 1st component only
 
644
cc*****************************************************************************
 
645
c      implicit none
 
646
cc
 
647
cc     Constants
 
648
cc     
 
649
c      include 'genps.inc'
 
650
c      double precision   zero
 
651
c      parameter (zero = 0d0)
 
652
cc
 
653
cc     Arguments
 
654
cc
 
655
c      double precision xprop(3,nexternal)
 
656
cc
 
657
cc     Local
 
658
cc
 
659
c      integer i,j,imin
 
660
c      double precision temp(3,nexternal),xmin
 
661
c      logical used(nexternal)
 
662
cc
 
663
cc     Global
 
664
cc
 
665
cc
 
666
cc     External
 
667
cc
 
668
cc-----
 
669
cc  Begin Code
 
670
cc-----      
 
671
c      do i=1,nexternal-3
 
672
c         used(i)=.false.
 
673
c      enddo
 
674
c      do j=1,nexternal-3
 
675
c         do i=1,nexternal-3
 
676
c            xmin = 2d0
 
677
c            if (.not. used(i) .and. xprop(1,i) .lt. xmin) then
 
678
c               xmin = xprop(1,i)
 
679
c               imin = i
 
680
c            endif
 
681
c         enddo
 
682
c         do i=1,3
 
683
c            temp(i,j)=xprop(i,imin)
 
684
c         enddo
 
685
c         used(i)=.true.
 
686
c      enddo
 
687
c      do i=1,nexternal-3
 
688
c         do j=1,3
 
689
c            xprop(j,i)=temp(j,i)
 
690
c         enddo
 
691
c      enddo
 
692
c      end
 
693
c
 
694
c
 
695
c      subroutine match_peak(tprop,xprop,imatch)
 
696
cc*****************************************************************************
 
697
cc     Determines if two sets of peaks are equivalent
 
698
cc*****************************************************************************
 
699
c      implicit none
 
700
cc
 
701
cc     Constants
 
702
cc     
 
703
c      include 'genps.inc'
 
704
c      double precision   zero
 
705
c      parameter (zero = 0d0)
 
706
cc
 
707
cc     Arguments
 
708
cc
 
709
c      double precision xprop(3,nexternal),tprop(3,nexternal)
 
710
c      integer imatch
 
711
cc
 
712
cc     Local
 
713
cc
 
714
c      integer i,j
 
715
cc
 
716
cc     Global
 
717
cc
 
718
cc
 
719
cc     External
 
720
cc
 
721
cc-----
 
722
cc  Begin Code
 
723
cc-----      
 
724
c      imatch = 1                     !By default assume match
 
725
c      do i=1,nexternal-3
 
726
c         do j=1,3
 
727
c            if (tprop(j,i) .ne. xprop(j,i)) imatch=0
 
728
c         enddo
 
729
c      enddo
 
730
c      end
 
731
 
 
732
c      subroutine get_peaks(iconfig,xt)
 
733
cc*****************************************************************************
 
734
cc     Attempts to determine peaks for this configuration
 
735
cc*****************************************************************************
 
736
c      implicit none
 
737
cc
 
738
cc     Constants
 
739
cc     
 
740
c      include 'genps.inc'
 
741
c      double precision   zero
 
742
c      parameter (zero = 0d0)
 
743
cc
 
744
cc     Arguments
 
745
cc
 
746
c      double precision xt(3,nexternal)
 
747
c      integer iconfig
 
748
c
 
749
cc
 
750
cc     Local
 
751
cc
 
752
c      double precision  xm(-nexternal:nexternal)
 
753
c      double precision tsgn, xo, a
 
754
c      double precision x1,x2
 
755
c      double precision dr,mtot, stot
 
756
c      integer i, l1, l2, j
 
757
c
 
758
c      double precision pmass(-nexternal:0,lmaxconfigs)
 
759
c      double precision pwidth(-nexternal:0,lmaxconfigs)
 
760
c      integer pow(-nexternal:0,lmaxconfigs)
 
761
c
 
762
c      integer imatch(0:lmaxconfigs)
 
763
c
 
764
cc
 
765
cc     Global
 
766
cc
 
767
c      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
768
c      common/to_forest/ iforest
 
769
c
 
770
c      integer            mapconfig(0:lmaxconfigs), this_config
 
771
c      common/to_mconfigs/mapconfig, this_config
 
772
c
 
773
c      real*8         emass(nexternal)
 
774
c      common/to_mass/emass
 
775
c
 
776
c      include 'run.inc'
 
777
c
 
778
c      double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
 
779
c      double precision emin(nincoming+1:nexternal)
 
780
c      double precision                    r2min(nincoming+1:nexternal,nincoming+1:nexternal)
 
781
c      double precision s_min(nexternal,nexternal)
 
782
c      common/to_cuts/  etmin, emin, etamax, r2min, s_min
 
783
c
 
784
c      double precision      spole(maxinvar),swidth(maxinvar),bwjac
 
785
c      common/to_brietwigner/spole          ,swidth          ,bwjac
 
786
c      
 
787
c      include 'coupl.inc'
 
788
cc      include 'props.inc'
 
789
c
 
790
cc
 
791
cc     External
 
792
cc
 
793
c
 
794
cc-----
 
795
cc  Begin Code
 
796
cc-----      
 
797
c      stot = 4d0*ebeam(1)*ebeam(2)
 
798
cc      iconfig = this_config
 
799
c      mtot = 0d0
 
800
c      do i=1,nexternal
 
801
c         xm(i)=emass(i)
 
802
c         mtot=mtot+xm(i)
 
803
c      enddo
 
804
c      tsgn    = +1d0
 
805
c      do i=-1,-(nexternal-3),-1              !Find all the propagotors
 
806
c         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
807
c         if (tsgn .eq. 1d0) then                         !s channel
 
808
c            xm(i) = xm(iforest(1,i,iconfig))+xm(iforest(2,i,iconfig))
 
809
c            mtot = mtot - xm(i)
 
810
c            if (pwidth(i,iconfig) .gt. 0) then    !B.W.
 
811
c               write(*,*) 'Setting BW',i,pmass(i,iconfig)
 
812
c               spole(-i)=pmass(i,iconfig)*pmass(i,iconfig)/stot
 
813
c               swidth(-i) = pwidth(i,iconfig)*pmass(i,iconfig)/stot
 
814
c               xm(i) = pmass(i,iconfig)
 
815
c               xt(1,-i) = spole(-i)
 
816
c               xt(2,-i) = swidth(-i)
 
817
c               xt(3,-i) = 2
 
818
c            else                                  !1/x^pow
 
819
c               if (xm(i) - pmass(i,iconfig) .le. 0d0) then !Can hit pole
 
820
cc                  write(*,*) 'Setting new min',i,xm(i),pmass(i,iconfig)
 
821
c                  l1 = iforest(1,i,iconfig)                  !need dr cut
 
822
c                  l2 = iforest(2,i,iconfig)
 
823
c                  if (l2 .lt. l1) then
 
824
c                     j = l1
 
825
c                     l1 = l2
 
826
c                     l2 = j
 
827
c                  endif
 
828
c                  dr = 0
 
829
cc-fax
 
830
c                  if (l1 .gt. 0) 
 
831
c     &  dr = max(r2min(l2,l1)*dabs(r2min(l2,l1)),0d0) !dr only for external
 
832
c                  dr = dr*.8d0                        !0.8 to hit peak hard
 
833
c                  xo = 0.5d0*pmass(i,iconfig)**2      !0.5 to hit peak hard
 
834
cc-fax
 
835
c                  if (dr .gt. 0d0) 
 
836
c     &            xo = max(xo,max(etmin(l2),0d0)*max(etmin(l1),0d0)*dr)
 
837
cc                  write(*,*) 'Got dr',i,l1,l2,dr
 
838
c                  xo = xo+pmass(i,iconfig)**2
 
839
cc                  write(*,*) 'Setting xm',i,xm(i),sqrt(xo)
 
840
c                  xm(i) = sqrt(xo)    !Reset xm to realistic minimum
 
841
c                  xo = xo/stot
 
842
cc                  xo = sqrt(pmass(i,iconfig)**2+ pmass(i,iconfig)**2)
 
843
cc-fax
 
844
cc                  xo = pmass(i,iconfig)+max(etmin,0d0)
 
845
c               else
 
846
cc                  write(*,*) 'Using xm',i,xm(i)
 
847
c                  xo = xm(i)**2/stot
 
848
c               endif
 
849
c               a=pmass(i,iconfig)**2/stot
 
850
cc               call setgrid(-i,xo,a,pow(i,iconfig))
 
851
cc               call setgrid(-i,xo,a,1)
 
852
c               xt(1,-i) = xo
 
853
c               xt(2,-i) = a
 
854
c               xt(3,-i) = 1
 
855
c            endif
 
856
c            mtot=mtot+xm(i)
 
857
cc            write(*,*) 'New mtot',i,mtot,xm(i)
 
858
c         else                                        !t channel
 
859
cc
 
860
cc     Check closest to p1
 
861
cc
 
862
c            l2 = iforest(2,i,iconfig) !need dr cut
 
863
c            x1 = 0
 
864
cc-fax
 
865
c            if (l2 .gt. 0) x1 = max(etmin(l2),0d0)
 
866
c            x1 = max(x1, xm(l2)/2d0)
 
867
cc            write(*,*) 'Using 1',l2,x1
 
868
c
 
869
cc
 
870
cc     Check closest to p2
 
871
cc
 
872
c            j = i-1
 
873
c            l2 = iforest(2,j,iconfig)
 
874
c            x2 = 0
 
875
cc-fax
 
876
c            if (l2 .gt. 0) x2 = max(etmin(l2),0d0)
 
877
c            x2 = max(x2, xm(l2)/2d0)
 
878
c            
 
879
cc            write(*,*) 'Using 2',l2,x2
 
880
c
 
881
c            xo = min(x1,x2)
 
882
c
 
883
c            xo = xo*xo/stot
 
884
c            a=-pmass(i,iconfig)**2/stot
 
885
cc            call setgrid(-i,xo,a,pow(i,iconfig))
 
886
cc            call setgrid(-i,xo,a,1)
 
887
c               xt(1,-i) = xo
 
888
c               xt(2,-i) = a
 
889
c               xt(3,-i) = 1
 
890
c         endif
 
891
c      enddo
 
892
cc---------------------
 
893
cc     tjs routine for x-hat
 
894
cc------------------------
 
895
c      if (abs(lpp(1)) .eq. 1 .or. abs(lpp(2)) .eq. 1) then
 
896
c         write(*,*) "setting s_hat",mtot,sqrt(stot)
 
897
c         xo = mtot**2/stot
 
898
c         i = 3*(nexternal-2) - 4 + 1
 
899
cc         call setgrid(i,xo,0d0,1)
 
900
c      endif
 
901
c      end
 
902
 
 
903
 
 
904
c       subroutine writeamp(p)
 
905
cc
 
906
cc     Constants
 
907
cc     
 
908
c      include 'genps.inc'
 
909
cc
 
910
cc     Arguments
 
911
cc
 
912
c      double precision p(0:3,nexternal)
 
913
cc
 
914
cc     Local
 
915
cc
 
916
c      double precision xp(0:3,-nexternal:nexternal)
 
917
c      integer i,j,iconfig
 
918
cc
 
919
cc     Global
 
920
cc
 
921
c      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
922
c      common/to_forest/ iforest
 
923
c
 
924
c      integer            mapconfig(0:lmaxconfigs), this_config
 
925
c      common/to_mconfigs/mapconfig, this_config
 
926
c
 
927
cc
 
928
cc     External
 
929
cc
 
930
c      double precision dot
 
931
cc-----
 
932
cc  Begin Code
 
933
cc-----
 
934
c      iconfig = this_config
 
935
c      do i=1,nexternal
 
936
c         do j=0,3
 
937
c            xp(j,i)=p(j,i)
 
938
c         enddo
 
939
c      enddo
 
940
c      shat = dot(p(0,1),p(0,2))
 
941
c      testamp = 1d0
 
942
c      tsgn    = +1d0
 
943
c      do i=-1,-(nexternal-3),-1              !Find all the propagotors
 
944
c         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
945
c         do j=0,3
 
946
c            xp(j,i) = xp(j,iforest(1,i,iconfig))
 
947
c     $           +tsgn*xp(j,iforest(2,i,iconfig))
 
948
c         enddo
 
949
cc         testamp=testamp*shat/(dot(xp(0,i),xp(0,i)))**2
 
950
cc         testamp=testamp**2
 
951
c      enddo
 
952
c      if(nexternal.gt.3)
 
953
c     $   write(*,'(A,4e15.5)') 'V',(dot(xp(0,-i),xp(0,-i)),i=1,nexternal-3)
 
954
cc      testamp=abs(testamp)
 
955
c      end
 
956
c
 
957
c
 
958
c      subroutine histamp(p,dwgt)
 
959
cc
 
960
cc     Constants
 
961
cc     
 
962
c      include 'genps.inc'
 
963
cc
 
964
cc     Arguments
 
965
cc
 
966
c      double precision p(0:3,nexternal),dwgt
 
967
cc
 
968
cc     Local
 
969
cc
 
970
c      double precision xp(0:3,-nexternal:nexternal)
 
971
c      integer i,j,iconfig
 
972
c      real wgt
 
973
cc
 
974
cc     Global
 
975
cc
 
976
c      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
977
c      common/to_forest/ iforest
 
978
c
 
979
c      integer            mapconfig(0:lmaxconfigs), this_config
 
980
c      common/to_mconfigs/mapconfig, this_config
 
981
c
 
982
cc
 
983
cc     External
 
984
cc
 
985
c      double precision dot
 
986
cc-----
 
987
cc  Begin Code
 
988
cc-----
 
989
c      wgt = dwgt
 
990
c      iconfig = this_config
 
991
c      do i=1,nexternal
 
992
c         do j=0,3
 
993
c            xp(j,i)=p(j,i)
 
994
c         enddo
 
995
c      enddo
 
996
c      shat = dot(p(0,1),p(0,2))
 
997
c      testamp = 1d0
 
998
c      tsgn    = +1d0
 
999
c      do i=-1,-(nexternal-3),-1              !Find all the propagotors
 
1000
c         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
1001
c         do j=0,3
 
1002
c            xp(j,i) = xp(j,iforest(1,i,iconfig))
 
1003
c     $           +tsgn*xp(j,iforest(2,i,iconfig))
 
1004
c         enddo
 
1005
c      enddo
 
1006
c      do i=1,nexternal-3
 
1007
cc         write(*,*) sqrt(abs(dot(xp(0,-i),xp(0,-i)))),wgt
 
1008
c         call hfill(10+i,real(sqrt(abs(dot(xp(0,-i),xp(0,-i))))),0.,wgt)
 
1009
c      enddo
 
1010
c      end
 
1011
c
 
1012
c
 
1013
c      subroutine check_limits(p,xlimit, iconfig)
 
1014
cc*************************************************************************
 
1015
cc     Checks limits on all of the functions being integrated
 
1016
cc*************************************************************************
 
1017
c      implicit none
 
1018
cc
 
1019
cc     Constants
 
1020
cc
 
1021
c      include 'genps.inc'
 
1022
cc
 
1023
cc     Arguments
 
1024
cc
 
1025
c      double precision p(0:3,nexternal), xlimit(2,nexternal)
 
1026
c      integer iconfig
 
1027
cc
 
1028
cc     Local
 
1029
cc
 
1030
c      double precision xp(0:3,-nexternal:nexternal), tsgn
 
1031
c      double precision sm
 
1032
c      integer ik(4)
 
1033
c      integer i,j, k1,k2
 
1034
cc
 
1035
cc     Global
 
1036
cc
 
1037
c      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
1038
c      common/to_forest/ iforest
 
1039
c      integer            mapconfig(0:lmaxconfigs), this_config
 
1040
c      common/to_mconfigs/mapconfig, this_config
 
1041
cc
 
1042
cc     External
 
1043
cc
 
1044
c      double precision dot
 
1045
c      external dot
 
1046
c
 
1047
cc      data ik/1,2,3,0/
 
1048
cc-----
 
1049
cc  Begin Code
 
1050
cc-----
 
1051
cc
 
1052
cc     Transform from rambo(1:4) format to helas (0:3)
 
1053
cc     
 
1054
c      do i=1,nexternal
 
1055
c         do j=0,3
 
1056
c            xp(j,i)=p(j,i)
 
1057
c         enddo
 
1058
c      enddo
 
1059
cc
 
1060
cc     Now build propagators
 
1061
cc      
 
1062
c      tsgn=+1d0
 
1063
c      do i=-1,-(nexternal-3),-1
 
1064
c         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
1065
c         k1=iforest(1,i,iconfig)
 
1066
c         k2=iforest(2,i,iconfig)
 
1067
c         do j=0,3
 
1068
c            xp(j,i)=xp(j,k1)+tsgn*xp(j,k2)
 
1069
c         enddo
 
1070
c         sm = tsgn*dot(xp(0,i),xp(0,i))
 
1071
c         if (sm .lt. xlimit(1,-i)) then
 
1072
c            xlimit(1,-i)=sm
 
1073
cc            write(*,*) 'New limit',-i,sm
 
1074
c         endif
 
1075
c         if (sm .gt. xlimit(2,-i)) then
 
1076
c            xlimit(2,-i)=sm
 
1077
cc            write(*,*) 'New limit',-i,sm
 
1078
c         endif
 
1079
c      enddo
 
1080
c      end
 
1081
c