← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_427003/PROC_427003/SubProcesses/myamp.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_427003/PROC_427003/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 'maxconfigs.inc'
 
11
      include 'nexternal.inc'
 
12
      double precision   zero
 
13
      parameter (zero = 0d0)
 
14
c
 
15
c     Arguments
 
16
c
 
17
      double precision p(0:3,nexternal)
 
18
c      integer iconfig
 
19
c
 
20
c     Local
 
21
c
 
22
      double precision xp(0:3,-nexternal:nexternal)
 
23
      double precision mpole(-nexternal:0),shat,tsgn
 
24
      integer i,j,iconfig
 
25
 
 
26
      double precision prmass(-nexternal:0,lmaxconfigs)
 
27
      double precision prwidth(-nexternal:0,lmaxconfigs)
 
28
      integer pow(-nexternal:0,lmaxconfigs)
 
29
      logical first_time
 
30
c
 
31
c     Global
 
32
c
 
33
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
34
      common/to_forest/ iforest
 
35
      integer            mapconfig(0:lmaxconfigs), this_config
 
36
      common/to_mconfigs/mapconfig, this_config
 
37
      
 
38
      include 'coupl.inc'
 
39
c
 
40
c     External
 
41
c
 
42
      double precision dot
 
43
 
 
44
      save prmass,prwidth,pow
 
45
      data first_time /.true./
 
46
c-----
 
47
c  Begin Code
 
48
c-----      
 
49
      iconfig = this_config
 
50
      if (first_time) then
 
51
c         include 'props.inc'
 
52
         first_time=.false.
 
53
      endif
 
54
 
 
55
      do i=1,nexternal
 
56
         mpole(-i)=0d0
 
57
         do j=0,3
 
58
            xp(j,i)=p(j,i)
 
59
         enddo
 
60
      enddo
 
61
c      mpole(-3) = 174**2
 
62
c      shat = dot(p(0,1),p(0,2))/(1800)**2
 
63
      shat = dot(p(0,1),p(0,2))/(10)**2
 
64
c      shat = 1d0
 
65
      testamp = 1d0
 
66
      tsgn    = +1d0
 
67
      do i=-1,-(nexternal-3),-1              !Find all the propagotors
 
68
         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
69
         do j=0,3
 
70
            xp(j,i) = xp(j,iforest(1,i,iconfig))
 
71
     $           +tsgn*xp(j,iforest(2,i,iconfig))
 
72
         enddo
 
73
         if (prwidth(i,iconfig) .ne. 0d0 .and. .false.) then
 
74
            testamp=testamp/((dot(xp(0,i),xp(0,i))
 
75
     $                        -prmass(i,iconfig)**2)**2
 
76
     $         -(prmass(i,iconfig)*prwidth(i,iconfig))**2)
 
77
         else
 
78
            testamp = testamp/((dot(xp(0,i),xp(0,i)) -
 
79
     $                          prmass(i,iconfig)**2)
 
80
     $                          **(pow(i,iconfig)))
 
81
         endif
 
82
        testamp=testamp*shat**(pow(i,iconfig))
 
83
c        write(*,*) i,iconfig,pow(i,iconfig),prmass(i,iconfig)
 
84
      enddo
 
85
c      testamp = 1d0/dot(xp(0,-1),xp(0,-1))
 
86
      testamp=abs(testamp)
 
87
c      testamp = 1d0
 
88
      end
 
89
 
 
90
      logical function cut_bw(p)
 
91
c*****************************************************************************
 
92
c     Approximates matrix element by propagators
 
93
c*****************************************************************************
 
94
      implicit none
 
95
c
 
96
c     Constants
 
97
c     
 
98
      include 'genps.inc'
 
99
      include 'maxconfigs.inc'
 
100
      include 'nexternal.inc'
 
101
      double precision   zero
 
102
      parameter (zero = 0d0)
 
103
      include 'run.inc'
 
104
c
 
105
c     Arguments
 
106
c
 
107
      double precision p(0:3,nexternal)
 
108
c
 
109
c     Local
 
110
c
 
111
      double precision xp(0:3,-nexternal:nexternal)
 
112
      double precision mpole(-nexternal:0),shat,tsgn
 
113
      integer i,j,iconfig,iproc
 
114
 
 
115
      double precision prmass(-nexternal:0,lmaxconfigs)
 
116
      double precision prwidth(-nexternal:0,lmaxconfigs)
 
117
      integer pow(-nexternal:0,lmaxconfigs)
 
118
      logical first_time, onshell
 
119
      double precision xmass
 
120
      integer nbw
 
121
 
 
122
      integer ida(2),idenpart
 
123
c
 
124
c     Global
 
125
c
 
126
      include 'maxamps.inc'
 
127
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
128
      common/to_forest/ iforest
 
129
      integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
 
130
      integer tprid(-max_branch:-1,lmaxconfigs)
 
131
      common/to_sprop/sprop,tprid
 
132
      integer            mapconfig(0:lmaxconfigs), this_config
 
133
      common/to_mconfigs/mapconfig, this_config
 
134
 
 
135
      integer        lbw(0:nexternal)  !Use of B.W.
 
136
      common /to_BW/ lbw
 
137
 
 
138
      logical             OnBW(-nexternal:0)     !Set if event is on B.W.
 
139
      common/to_BWEvents/ OnBW
 
140
      
 
141
      include 'coupl.inc'
 
142
 
 
143
      integer idup(nexternal,maxproc,maxsproc)
 
144
      integer mothup(2,nexternal)
 
145
      integer icolup(2,nexternal,maxflow,maxsproc)
 
146
      include 'leshouche.inc'
 
147
 
 
148
      integer gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
 
149
      include 'decayBW.inc'
 
150
c
 
151
c     External
 
152
c
 
153
      double precision dot
 
154
 
 
155
      save prmass,prwidth,pow
 
156
      data first_time /.true./
 
157
c-----
 
158
c  Begin Code
 
159
c-----      
 
160
      cut_bw = .false.    !Default is we passed the cut
 
161
      iconfig = this_config
 
162
      if (first_time) then
 
163
         include 'props.inc'
 
164
         nbw = 0
 
165
         do i=-1,-(nexternal-3),-1
 
166
            if (iforest(1,i,iconfig) .eq. 1 .or. prwidth(i,iconfig).le.0) then
 
167
              cycle
 
168
            endif
 
169
            nbw=nbw+1
 
170
            if (lbw(nbw) .eq. 1) then
 
171
               write(*,*) 'Requiring BW ',i,nbw
 
172
            elseif(lbw(nbw) .eq. 2) then
 
173
               write(*,*) 'Excluding BW ',i,nbw
 
174
            else
 
175
               write(*,*) 'No cut BW ',i,nbw
 
176
            endif
 
177
         enddo
 
178
         first_time=.false.
 
179
      endif
 
180
 
 
181
      do i=1,nexternal
 
182
         mpole(-i)=0d0
 
183
         do j=0,3
 
184
            xp(j,i)=p(j,i)
 
185
         enddo
 
186
      enddo
 
187
      nbw = 0
 
188
      tsgn    = +1d0
 
189
c     Find non-zero process number
 
190
      do iproc=1,maxsproc
 
191
         if(sprop(iproc,-1,iconfig).ne.0) goto 10
 
192
      enddo
 
193
 10   continue
 
194
c     If no non-zero sprop, set iproc to 1
 
195
      if(iproc.gt.maxsproc) iproc=1
 
196
c     Start loop over propagators
 
197
      do i=-1,-(nexternal-3),-1
 
198
         onbw(i) = .false.
 
199
         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
200
         do j=0,3
 
201
            xp(j,i) = xp(j,iforest(1,i,iconfig))
 
202
     $           +tsgn*xp(j,iforest(2,i,iconfig))
 
203
         enddo
 
204
         if (tsgn .lt. 0d0) cycle
 
205
         if (prwidth(i,iconfig) .gt. 0d0 ) then !This is B.W.
 
206
            nbw=nbw+1
 
207
c            write(*,*) 'Checking BW',nbw
 
208
            xmass = sqrt(dot(xp(0,i),xp(0,i)))
 
209
c            write(*,*) 'xmass',xmass,prmass(i,iconfig)
 
210
c
 
211
c           Here we set if the BW is "on-shell" for LesHouches
 
212
c
 
213
            onshell = (abs(xmass - prmass(i,iconfig)) .lt.
 
214
     $           bwcutoff*prwidth(i,iconfig))
 
215
            if(onshell)then
 
216
c     Remove on-shell forbidden s-channels (gForceBW=2) (JA 2/10/11)
 
217
              if(gForceBW(i,iconfig).eq.2) then
 
218
                 cut_bw = .true.
 
219
                 return               
 
220
              endif
 
221
c           Only allow OnBW if no "decay" to identical particle
 
222
              OnBW(i) = .true.
 
223
              idenpart=0
 
224
              do j=1,2
 
225
                ida(j)=iforest(j,i,iconfig)
 
226
                if(ida(j).lt.0) then
 
227
                   if(sprop(iproc,i,iconfig).eq.sprop(iproc,ida(j),iconfig))
 
228
     $                  idenpart=ida(j)
 
229
                elseif (ida(j).gt.0) then
 
230
                   if(sprop(iproc,i,iconfig).eq.IDUP(ida(j),1,iproc))
 
231
     $                  idenpart=ida(j)
 
232
                endif
 
233
              enddo
 
234
c           Always remove if daughter final-state
 
235
              if(idenpart.gt.0) then
 
236
                 OnBW(i)=.false.
 
237
c           Else remove if daughter forced to be onshell
 
238
              elseif(idenpart.lt.0)then
 
239
                 if(gForceBW(idenpart, iconfig).eq.1) then
 
240
                    OnBW(i)=.false.
 
241
c           Else remove daughter if forced to be onshell
 
242
                 elseif(gForceBW(i, iconfig).eq.1) then
 
243
                    OnBW(idenpart)=.false.
 
244
c           Else remove either this resonance or daughter, which is closer to mass shell
 
245
                 elseif(abs(xmass-prmass(i,iconfig)).gt.
 
246
     $                   abs(sqrt(dot(xp(0,idenpart),xp(0,idenpart)))-
 
247
     $                   prmass(i,iconfig))) then
 
248
                    OnBW(i)=.false.
 
249
c           Else remove OnBW for daughter
 
250
                 else
 
251
                    OnBW(idenpart)=.false.
 
252
                 endif
 
253
              endif
 
254
            else if (gForceBW(i, iconfig).eq.1) then ! .not. onshell
 
255
c             Check if we are supposed to cut forced bw (JA 4/8/11)
 
256
              cut_bw = .true.
 
257
c              write(*,*) 'cut_bw: ',i,gForceBW(i,iconfig),OnBW(i),cut_bw
 
258
              return
 
259
            endif
 
260
c
 
261
c     Here we set onshell for phase space integration (JA 4/8/11)
 
262
c
 
263
            onshell = (abs(xmass - prmass(i,iconfig)) .lt.
 
264
     $           5d0*prwidth(i,iconfig))
 
265
 
 
266
            if (onshell .and. (lbw(nbw).eq. 2) .or.
 
267
     $          .not. onshell .and. (lbw(nbw).eq. 1)) then
 
268
               cut_bw=.true.
 
269
c               write(*,*) 'cut_bw: ',nbw,xmass,onshell,lbw(nbw),cut_bw
 
270
               return
 
271
            endif
 
272
         endif
 
273
c         write(*,*) 'final cut_bw: ',nbw,lbw(nbw),xmass,onshell,OnBW(i),cut_bw
 
274
      enddo
 
275
      end
 
276
 
 
277
 
 
278
      subroutine set_peaks
 
279
c*****************************************************************************
 
280
c     Attempts to determine peaks for this configuration
 
281
c*****************************************************************************
 
282
      implicit none
 
283
c
 
284
c     Constants
 
285
c     
 
286
      include 'genps.inc'
 
287
      include 'maxconfigs.inc'
 
288
      include 'nexternal.inc'
 
289
      include 'maxamps.inc'
 
290
      double precision   zero
 
291
      parameter (zero = 0d0)
 
292
c
 
293
c     Arguments
 
294
c
 
295
c
 
296
c     Local
 
297
c
 
298
      double precision  xm(-nexternal:nexternal)
 
299
      double precision  xe(-nexternal:nexternal)
 
300
      double precision tsgn, xo, a
 
301
      double precision x1,x2,xk(nexternal)
 
302
      double precision dr,mtot,etot,xqfact
 
303
      double precision spmass
 
304
      integer i, iconfig, l1, l2, j, nt, nbw, iproc
 
305
      integer iden_part(-nexternal+1:nexternal)
 
306
 
 
307
      double precision prmass(-nexternal:0,lmaxconfigs)
 
308
      double precision prwidth(-nexternal:0,lmaxconfigs)
 
309
      integer pow(-nexternal:0,lmaxconfigs)
 
310
 
 
311
      integer idup(nexternal,maxproc,maxsproc)
 
312
      integer mothup(2,nexternal)
 
313
      integer icolup(2,nexternal,maxflow,maxsproc)
 
314
      include 'leshouche.inc'
 
315
 
 
316
      integer gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
 
317
      include 'decayBW.inc'
 
318
 
 
319
c
 
320
c     Global
 
321
c
 
322
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
323
      common/to_forest/ iforest
 
324
 
 
325
      integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
 
326
      integer tprid(-max_branch:-1,lmaxconfigs)
 
327
      common/to_sprop/sprop,tprid
 
328
 
 
329
      integer            mapconfig(0:lmaxconfigs), this_config
 
330
      common/to_mconfigs/mapconfig, this_config
 
331
 
 
332
      real*8         emass(nexternal)
 
333
      common/to_mass/emass
 
334
 
 
335
      include 'run.inc'
 
336
 
 
337
      double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
 
338
      double precision emin(nincoming+1:nexternal)
 
339
      double precision                    r2min(nincoming+1:nexternal,nincoming+1:nexternal)
 
340
      double precision s_min(nexternal,nexternal)
 
341
      common/to_cuts/  etmin, emin, etamax, r2min, s_min
 
342
 
 
343
      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
 
344
      common/to_xqcuts/xqcutij,xqcuti
 
345
 
 
346
      double precision      spole(maxinvar),swidth(maxinvar),bwjac
 
347
      common/to_brietwigner/spole          ,swidth          ,bwjac
 
348
 
 
349
      integer        lbw(0:nexternal)  !Use of B.W.
 
350
      common /to_BW/ lbw
 
351
 
 
352
      double precision stot,m1,m2
 
353
      common/to_stot/stot,m1,m2
 
354
 
 
355
      include 'coupl.inc'
 
356
      include 'cuts.inc'
 
357
c
 
358
c     External
 
359
c
 
360
 
 
361
c-----
 
362
c  Begin Code
 
363
c-----      
 
364
      include 'props.inc'
 
365
c      etmin = 10
 
366
      nt = 0
 
367
      iconfig = this_config
 
368
      mtot = 0d0
 
369
      etot = 0d0   !Total energy needed
 
370
      spmass = 0d0 !Keep track of BW masses for shat
 
371
      xqfact=1d0
 
372
      if(ickkw.eq.2.or.ktscheme.eq.2) xqfact=0.3d0
 
373
      do i=nincoming+1,nexternal  !assumes 2 incoming
 
374
         xm(i)=emass(i)
 
375
c-fax
 
376
         xe(i)=max(emass(i),max(etmin(i),0d0))
 
377
         xe(i)=max(xe(i),max(emin(i),0d0))
 
378
c-JA 1/2009: Set grid also based on xqcut
 
379
         xe(i)=max(xe(i),xqfact*xqcuti(i))
 
380
         xk(i)= 0d0
 
381
         etot = etot+xe(i)
 
382
         mtot=mtot+xm(i)         
 
383
      enddo
 
384
      spmass=mtot
 
385
      tsgn    = +1d0
 
386
c     Reset variables
 
387
      nbw = 0
 
388
      do i=1,nexternal-2
 
389
         spole(i)=0
 
390
         swidth(i)=0
 
391
      enddo
 
392
c     Find non-zero process number
 
393
      do iproc=1,maxsproc
 
394
         if(sprop(iproc,-1,iconfig).ne.0) goto 10
 
395
      enddo
 
396
 10   continue
 
397
c     If no non-zero sprop, set iproc to 1
 
398
      if(iproc.ge.maxsproc.and.sprop(maxsproc,-1,iconfig).eq.0)
 
399
     $     iproc=1
 
400
 
 
401
c     Look for identical particles to map radiation processes
 
402
      call idenparts(iden_part, iforest(1,-max_branch,iconfig),
 
403
     $     sprop(1,-max_branch,iconfig), gForceBW(-max_branch,iconfig),
 
404
     $     prwidth(-nexternal,iconfig))
 
405
 
 
406
c     Start loop over propagators
 
407
      do i=-1,-(nexternal-3),-1
 
408
         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
 
409
         if (tsgn .eq. 1d0) then                         !s channel
 
410
            xm(i) = xm(iforest(1,i,iconfig))+xm(iforest(2,i,iconfig))
 
411
            xe(i) = xe(iforest(1,i,iconfig))+xe(iforest(2,i,iconfig))
 
412
            mtot = mtot - xm(i)
 
413
            etot = etot - xe(i)
 
414
            if (iforest(1,i,iconfig) .gt. 0
 
415
     &           .and. iforest(2,i,iconfig) .gt. 0) then
 
416
c-JA 1/2009: Set deltaR cuts here together with s_min cuts
 
417
              l1 = iforest(1,i,iconfig)
 
418
              l2 = iforest(2,i,iconfig)
 
419
              xm(i)=max(xm(i),sqrt(max(s_min(l1,l2),0d0)))
 
420
              dr = max(r2min(l1,l2)*dabs(r2min(l1,l2)),0d0)*0.8d0
 
421
              xm(i)=max(xm(i),
 
422
     &           sqrt(max(etmin(l2),0d0)*max(etmin(l1),0d0)*dr))
 
423
c-JA 1/2009: Set grid also based on xqcut
 
424
              xm(i)=max(xm(i),max(xqcutij(l1,l2),0d0))
 
425
            endif
 
426
c            write(*,*) 'iconfig,i',iconfig,i
 
427
c            write(*,*) prwidth(i,iconfig),prmass(i,iconfig)
 
428
            if (prwidth(i,iconfig) .gt. 0 ) then
 
429
               nbw=nbw+1
 
430
c              JA 6/8/2011 Set xe(i) for resonances
 
431
               if (lbw(nbw).eq.1) then
 
432
                  xm(i) = max(xm(i), prmass(i,iconfig)-5d0*prwidth(i,iconfig))
 
433
               else if (gforcebw(i,iconfig).eq.1) then
 
434
                  xm(i) = max(xm(i), prmass(i,iconfig)-bwcutoff*prwidth(i,iconfig))
 
435
               endif
 
436
            endif
 
437
            xe(i)=max(xe(i),xm(i))
 
438
c     Check for impossible onshell configurations
 
439
c     Either: required onshell and daughter masses too large
 
440
c     Or: forced and daughter masses too large
 
441
c     Or: required offshell and forced, with bwcutoff.le.5
 
442
            if(prwidth(i,iconfig) .gt. 0.and.
 
443
     $         (lbw(nbw).eq.1.and.
 
444
     $          (prmass(i,iconfig)+5d0*prwidth(i,iconfig).lt.xm(i)
 
445
     $           .or.prmass(i,iconfig)-5d0*prwidth(i,iconfig).gt.stot)
 
446
     $          .or.gforcebw(i,iconfig).eq.1.and.
 
447
     $              prmass(i,iconfig)+bwcutoff*prwidth(i,iconfig).lt.xm(i)
 
448
     $          .or.lbw(nbw).eq.2.and.gforcebw(i,iconfig).eq.1 .and.
 
449
     $              bwcutoff.le.5d0))
 
450
     $        then
 
451
c     Write results.dat and quit
 
452
               call write_null_results()
 
453
               stop
 
454
            endif
 
455
            if (prwidth(i,iconfig) .gt. 0 .and. lbw(nbw) .le. 1) then         !B.W.
 
456
               if (i .eq. -(nexternal-(nincoming+1))) then  !This is s-hat
 
457
                  j = 3*(nexternal-2)-4+1    !set i to ndim+1
 
458
c-----
 
459
c tjs 11/2008 if require BW then force even if worried about energy
 
460
c JA 8/2011 don't use BW if mass is > CM energy
 
461
c----
 
462
                  if(prmass(i,iconfig).ge.xm(i).and.iden_part(i).eq.0.and.
 
463
     $                 prmass(i,iconfig).lt.sqrt(stot)
 
464
     $                 .or. lbw(nbw).eq.1) then
 
465
                     write(*,*) 'Setting PDF BW',j,nbw,prmass(i,iconfig)
 
466
                     spole(j)=prmass(i,iconfig)*prmass(i,iconfig)/stot
 
467
                     swidth(j) = prwidth(i,iconfig)*prmass(i,iconfig)/stot
 
468
                  endif
 
469
               else if((prmass(i,iconfig)+5d0*prwidth(i,iconfig)).ge.xm(i)
 
470
     $                  .and. iden_part(i).eq.0 .or. lbw(nbw).eq.1) then
 
471
c              JA 02/13 Only allow BW if xm below M+5*Gamma
 
472
                  write(*,*) 'Setting BW',i,nbw,prmass(i,iconfig)
 
473
                  spole(-i)=prmass(i,iconfig)*prmass(i,iconfig)/stot
 
474
                  swidth(-i) = prwidth(i,iconfig)*prmass(i,iconfig)/stot
 
475
               endif
 
476
c     JA 4/1/2011 Set grid in case there is no BW (radiation process)
 
477
               if (swidth(-i) .eq. 0d0 .and.
 
478
     $              i.ne.-(nexternal-(nincoming+1)))then
 
479
                  a=prmass(i,iconfig)**2/stot
 
480
                  xo = min(xm(i)**2/stot, 1-1d-8)
 
481
                  if (xo.eq.0d0) xo=1d0/stot
 
482
                  call setgrid(-i,xo,a,1)
 
483
               endif
 
484
c     Set spmass for BWs
 
485
               if (swidth(-i) .ne. 0d0)
 
486
     $              spmass=spmass-xm(i) +
 
487
     $              max(xm(i),prmass(i,iconfig)-5d0*prwidth(i,iconfig))
 
488
            else                                  !1/x^pow
 
489
              a=prmass(i,iconfig)**2/stot
 
490
c     JA 4/1/2011 always set grid
 
491
              xo = min(xm(i)**2/stot, 1-1d-8)
 
492
              if (xo.eq.0d0) xo=1d0/stot
 
493
c              if (prwidth(i, iconfig) .eq. 0d0.or.iden_part(i).gt.0) then 
 
494
              call setgrid(-i,xo,a,1)
 
495
c              else 
 
496
c                 write(*,*) 'Using flat grid for BW',i,nbw,
 
497
c     $                prmass(i,iconfig)
 
498
c              endif
 
499
            endif
 
500
            etot = etot+xe(i)
 
501
            mtot=mtot+xm(i)
 
502
c            write(*,*) 'New mtot',i,mtot,xm(i)
 
503
         else                                        !t channel
 
504
c
 
505
c     Check closest to p1
 
506
c
 
507
            nt = nt+1
 
508
            l2 = iforest(2,i,iconfig) !need dr cut
 
509
            x1 = 0            
 
510
c-fax
 
511
c-JA 1/2009: Set grid also based on xqcut
 
512
            if (l2 .gt. 0) x1 = max(etmin(l2),max(xqfact*xqcuti(l2),0d0))
 
513
            x1 = max(x1, xe(l2)/1d0)
 
514
            if (nt .gt. 1) x1 = max(x1,xk(nt-1))
 
515
            xk(nt)=x1
 
516
c            write(*,*) 'Using 1',l2,x1
 
517
 
 
518
c
 
519
c     Check closest to p2
 
520
c
 
521
            j = i-1
 
522
            l2 = iforest(2,j,iconfig)
 
523
            x2 = 0
 
524
c-JA 1/2009: Set grid also based on xqcut
 
525
            if (l2 .gt. 0) x2 = max(etmin(l2),max(xqfact*xqcuti(l2),0d0))
 
526
c            if (l2 .gt. 0) x2 = max(etmin(l2),0d0)
 
527
            x2 = max(x2, xe(l2)/1d0)
 
528
c            if (nt .gt. 1) x2 = max(x2,xk(nt-1))
 
529
            
 
530
c            write(*,*) 'Using 2',l2,x2
 
531
 
 
532
            xo = min(x1,x2)
 
533
 
 
534
c           Use 1/10000 of sqrt(s) as minimum, to always get integration
 
535
            xo = xo*xo/stot
 
536
            if (xo.eq.0d0)then
 
537
               xo=1d0/stot
 
538
               write(*,*) 'Warning: No cutoff for shat integral found'
 
539
               write(*,*) '         Minimum set to ', xo
 
540
            endif
 
541
            a=-prmass(i,iconfig)**2/stot
 
542
c            call setgrid(-i,xo,a,pow(i,iconfig))
 
543
 
 
544
c               write(*,*) 'Enter minimum for ',-i, xo
 
545
c               read(*,*) xo
 
546
             if (i .ne. -1 .or. .true.) call setgrid(-i,xo,a,1)
 
547
         endif
 
548
      enddo
 
549
c     Perform setting for shat (PDF BW or 1/s)
 
550
      if (abs(lpp(1)) .eq. 1 .or. abs(lpp(2)) .eq. 1) then
 
551
c     Set minimum based on: 1) required energy 2) resonances 3) 1/10000 of sqrt(s)
 
552
         i = 3*(nexternal-2) - 4 + 1
 
553
         xo = max(min(etot**2/stot, 1d0-1d-8),1d0/stot)
 
554
c        Take into account special cuts
 
555
         xo = max(xo, xptj*dabs(xptj)/stot)
 
556
         xo = max(xo, xptb*dabs(xptb)/stot)
 
557
         xo = max(xo, xpta*dabs(xpta)/stot)
 
558
         xo = max(xo, xptl*dabs(xptl)/stot)
 
559
         xo = max(xo, xmtc*dabs(xmtc)/stot)
 
560
         xo = max(xo, htjmin**2/stot)
 
561
         xo = max(xo, ptj1min**2/stot)
 
562
         xo = max(xo, (2*ptj2min)**2/stot)
 
563
         xo = max(xo, (3*ptj3min)**2/stot)
 
564
         xo = max(xo, (4*ptj4min)**2/stot)
 
565
         xo = max(xo, ht2min**2/stot)
 
566
         xo = max(xo, ht3min**2/stot)
 
567
         xo = max(xo, ht4min**2/stot)
 
568
         xo = max(xo, misset**2/stot)
 
569
         xo = max(xo, ptllmin**2/stot)
 
570
         xo = max(xo, ptl1min**2/stot)
 
571
         xo = max(xo, (2*ptl2min)**2/stot)
 
572
         xo = max(xo, (3*ptl3min)**2/stot)
 
573
         xo = max(xo, (4*ptl4min)**2/stot)
 
574
         xo = max(xo, mmnl**2/stot)
 
575
c     Include mass scale from BWs
 
576
         xo = max(xo, spmass**2/stot)
 
577
         if (swidth(i).eq.0.and.xo.eq.1d0/stot) then
 
578
            write(*,*) 'Warning: No minimum found for integration'
 
579
            write(*,*) '         Setting minimum to ',1d0/stot
 
580
         endif
 
581
c-----------------------
 
582
c     tjs  4/29/2008 use analytic transform for s-hat
 
583
c-----------------------
 
584
         if (swidth(i) .eq. 0d0) then
 
585
            swidth(i) = xo
 
586
            spole(i)= -2.0d0    ! 1/s pole
 
587
            write(*,*) "Transforming s_hat 1/s ",i,xo
 
588
         else
 
589
            write(*,*) "Transforming s_hat BW ",spole(i),swidth(i)
 
590
         endif
 
591
      endif
 
592
 
 
593
      i=-8
 
594
c      write(*,*) 'Enter minimum for ',-i, xo
 
595
c      read(*,*) xo      
 
596
c      if (xo .gt. 0)      call setgrid(-i,xo,a,1)
 
597
 
 
598
      i=-10
 
599
c      write(*,*) 'Enter minimum for ',-i, xo
 
600
c      read(*,*) xo      
 
601
c      if (xo .gt. 0) call setgrid(-i,xo,a,1)
 
602
 
 
603
      end
 
604
 
 
605
      subroutine write_null_results()
 
606
      implicit none
 
607
 
 
608
      write(*,*),'Impossible BW configuration'
 
609
      open(unit=66,file='results.dat',status='unknown')
 
610
      write(66,'(3e12.5,2i9,i5,i9,3e10.3)')0.,0.,0.,0,0,1,0,0.,0.,0.
 
611
      write(66,'(i4,5e15.5)') 1,0.,0.,0.,0.,0.
 
612
      close(66)
 
613
      end