← 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/setcuts.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/setcuts.f
 
1
      SUBROUTINE SETCUTS
 
2
C**************************************************************************
 
3
C     SET THE CUTS 
 
4
C**************************************************************************
 
5
      IMPLICIT NONE
 
6
c
 
7
c     INCLUDE
 
8
c
 
9
      include 'genps.inc'
 
10
      include 'nexternal.inc'
 
11
      include 'coupl.inc'
 
12
      include 'run.inc'
 
13
      include 'cuts.inc'
 
14
c
 
15
c     Constants
 
16
c
 
17
      double precision zero
 
18
      parameter       (ZERO = 0d0)
 
19
      real*8 Pi
 
20
      parameter( Pi = 3.14159265358979323846d0 )
 
21
      integer    lun
 
22
      parameter (lun=22)
 
23
c
 
24
c     LOCAL
 
25
c
 
26
      integer i,j
 
27
      integer icollider,detail_level
 
28
      logical  do_cuts(nexternal)
 
29
      integer ncheck
 
30
      logical done,fopened
 
31
      logical from_decay(-(nexternal+3):nexternal)
 
32
C     
 
33
C     GLOBAL
 
34
C
 
35
c--masses and poles
 
36
      double precision pmass(nexternal)
 
37
      common/to_mass/  pmass
 
38
      DOUBLE PRECISION       qmass(2)
 
39
      COMMON/to_qmass/qmass
 
40
      double precision      spole(maxinvar),swidth(maxinvar),bwjac
 
41
      common/to_brietwigner/spole          ,swidth          ,bwjac
 
42
c--cuts
 
43
      double precision etmin(nincoming+1:nexternal)
 
44
      double precision etamax(nincoming+1:nexternal)
 
45
      double precision emin(nincoming+1:nexternal)
 
46
      double precision r2min(nincoming+1:nexternal,nincoming+1:nexternal)
 
47
      double precision s_min(nexternal,nexternal)
 
48
      double precision etmax(nincoming+1:nexternal)
 
49
      double precision etamin(nincoming+1:nexternal)
 
50
      double precision emax(nincoming+1:nexternal)
 
51
      double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
 
52
      double precision s_max(nexternal,nexternal)
 
53
      double precision ptll_min(nexternal,nexternal),ptll_max(nexternal,nexternal)
 
54
      double precision inclHtmin,inclHtmax
 
55
      common/to_cuts/  etmin, emin, etamax, r2min, s_min,
 
56
     $     etmax, emax, etamin, r2max, s_max, ptll_min, ptll_max, inclHtmin,inclHtmax
 
57
 
 
58
      double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
 
59
      logical jetor
 
60
      common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
 
61
 
 
62
      double precision ptlmin4(4),ptlmax4(4)
 
63
      common/to_lepton_cuts/ ptlmin4,ptlmax4
 
64
 
 
65
      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
 
66
      common/to_xqcuts/xqcutij,xqcuti
 
67
c
 
68
c     les houches accord stuff to identify neutrinos
 
69
c
 
70
      include 'maxamps.inc'
 
71
      integer idup(nexternal,maxproc,maxsproc)
 
72
      integer mothup(2,nexternal)
 
73
      integer icolup(2,nexternal,maxflow,maxsproc)
 
74
      include 'leshouche.inc'
 
75
C
 
76
      LOGICAL  IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
 
77
      LOGICAL  IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
 
78
      LOGICAL  IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
 
79
      COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
 
80
     . IS_A_ONIUM 
 
81
c
 
82
c
 
83
c     reading parameters
 
84
      include '../../Source/run_config.inc'
 
85
      character*20 param(maxpara),value(maxpara)
 
86
      integer npara
 
87
c
 
88
c     setup masses for the final-state particles
 
89
c
 
90
      include 'pmass.inc'
 
91
      include 'qmass.inc'
 
92
 
 
93
C-----
 
94
C  BEGIN CODE
 
95
C-----
 
96
c
 
97
c     read the cuts from the run_card.dat - this should already be done in main program
 
98
c
 
99
c      call setrun
 
100
 
 
101
c
 
102
c     No pdfs for decay processes - set here since here we know the nincoming
 
103
c     Also set stot here, and use mass of incoming particle for ren scale
 
104
c
 
105
         if(nincoming.eq.1)then
 
106
            lpp(1)=0
 
107
            lpp(2)=0
 
108
            ebeam(1)=pmass(1)/2d0
 
109
            ebeam(2)=pmass(1)/2d0
 
110
            scale=pmass(1)
 
111
            fixed_ren_scale=.true.
 
112
            fixed_fac_scale=.true.
 
113
         endif
 
114
c
 
115
c     set ptj and s_min if xqcut and ktscheme = 1, to improve
 
116
c     integration speed, and set drjj and drjl to 0.
 
117
c
 
118
        if(xqcut.gt.0) then
 
119
           if(auto_ptj_mjj.and.ptj.ge.0d0.and.ktscheme.eq.1)then
 
120
            ptj=xqcut
 
121
            write(*,*) 'Warning! ptj set to xqcut=',xqcut,
 
122
     $            ' to improve integration efficiency'
 
123
            write(*,*) 'Note that this might affect non-radiated jets,'
 
124
            write(*,*) 'e.g. from decays. Use cut_decays=F in run_card.'
 
125
          else if(ptj.gt.xqcut)then
 
126
            ptj=0d0
 
127
            write(*,*) 'Warning! ptj set to 0 since xqcut > 0 and'
 
128
            write(*,*) '         auto_ptj_mjj = F or ktscheme > 1'
 
129
          endif
 
130
          if(auto_ptj_mjj.and.mmjj.ge.0d0)then
 
131
            mmjj=xqcut
 
132
            write(*,*) 'Warning! mmjj set to xqcut=',xqcut,
 
133
     $            ' to improve integration efficiency'
 
134
            write(*,*) 'Note that this might affect non-radiated jets,'
 
135
            write(*,*) 'e.g. from decays. Use cut_decays=F in run_card.'
 
136
          else if(mmjj.gt.xqcut)then
 
137
            mmjj=0d0
 
138
            write(*,*) 'Warning! mmjj set to 0 since xqcut > 0 and'
 
139
            write(*,*) '         auto_ptj_mjj = F'
 
140
          endif
 
141
          if(drjj.gt.0d0) then
 
142
            write(*,*) 'Warning! drjj > 0 with xqcut > 0, set to 0'
 
143
            drjj = 0d0
 
144
          endif
 
145
          if(drjl.gt.0d0) then
 
146
            write(*,*) 'Warning! drjl > 0 with xqcut > 0, set to 0'
 
147
            drjl = 0d0
 
148
          endif
 
149
        endif
 
150
 
 
151
c     Check which particles come from decays
 
152
      if(.not.cut_decays)
 
153
     $       call check_decay(from_decay)
 
154
 
 
155
c
 
156
c     check if I have to apply cuts on the particles
 
157
c
 
158
      do i=nincoming+1,nexternal
 
159
         do_cuts(i)=.true.
 
160
         if(nincoming.eq.1) do_cuts(i)=.false.
 
161
         if(.not.cut_decays.and.from_decay(i)) do_cuts(i)=.false.
 
162
         is_a_j(i)=.false.
 
163
         is_a_l(i)=.false.
 
164
         is_a_b(i)=.false.
 
165
         is_a_a(i)=.false.
 
166
         is_a_nu(i)=.false.
 
167
 
 
168
 
 
169
c-do not apply cuts to these
 
170
         if (pmass(i).gt.20d0)     do_cuts(i)=.false.  ! no cuts on top,W,Z,H
 
171
         if (abs(idup(i,1,1)).eq.12) do_cuts(i)=.false.  ! no cuts on ve ve~
 
172
         if (abs(idup(i,1,1)).eq.14) do_cuts(i)=.false.  ! no cuts on vm vm~
 
173
         if (abs(idup(i,1,1)).eq.16) do_cuts(i)=.false.  ! no cuts on vt vt~
 
174
c-flavor-jets
 
175
         if (abs(idup(i,1,1)).le.min(maxjetflavor,5)) then
 
176
              is_a_j(i)=.true.
 
177
c              write(*,*)'jet:ithe pdg is ',abs(idup(i,1,1)),' maxflavor=',maxjetflavor
 
178
         else if (abs(idup(i,1,1)).ge.maxjetflavor+1 .and. abs(idup(i,1,1)).le.5) then
 
179
              is_a_b(i)=.true.
 
180
c              write(*,*)'bjet:the pdg is ',abs(idup(i,1,1)),' maxflavor=',maxjetflavor
 
181
         endif
 
182
 
 
183
         if (abs(idup(i,1,1)).eq.21)  is_a_j(i)=.true. ! gluon is a jet
 
184
c-charged-leptons
 
185
         if (abs(idup(i,1,1)).eq.11)  is_a_l(i)=.true. ! e+  e-
 
186
         if (abs(idup(i,1,1)).eq.13)  is_a_l(i)=.true. ! mu+ mu-
 
187
         if (abs(idup(i,1,1)).eq.15)  is_a_l(i)=.true. ! ta+ ta-
 
188
c-b-quarks
 
189
c         if (abs(idup(i,1,1)).eq.5)   is_a_b(i)=.true. ! b b~
 
190
c-photon
 
191
         if (idup(i,1,1).eq.22)   is_a_a(i)=.true. ! photon
 
192
c-neutrino's for missing et
 
193
         if (abs(idup(i,1,1)).eq.12) is_a_nu(i)=.true.  ! no cuts on ve ve~
 
194
         if (abs(idup(i,1,1)).eq.14) is_a_nu(i)=.true.  ! no cuts on vm vm~
 
195
         if (abs(idup(i,1,1)).eq.16) is_a_nu(i)=.true.  ! no cuts on vt vt~
 
196
         if (pmass(i).gt.10d0)     is_heavy(i)=.true. ! heavy fs particle
 
197
c-onium
 
198
         if (idup(i,1,1).eq.441)   is_a_onium(i)=.true. ! charmonium
 
199
         if (idup(i,1,1).eq.10441)   is_a_onium(i)=.true. ! charmonium
 
200
         if (idup(i,1,1).eq.100441)   is_a_onium(i)=.true. ! charmonium
 
201
         if (idup(i,1,1).eq.443)   is_a_onium(i)=.true. ! charmonium
 
202
         if (idup(i,1,1).eq.10443)   is_a_onium(i)=.true. ! charmonium
 
203
         if (idup(i,1,1).eq.20443)   is_a_onium(i)=.true. ! charmonium
 
204
         if (idup(i,1,1).eq.100443)   is_a_onium(i)=.true. ! charmonium
 
205
         if (idup(i,1,1).eq.30443)   is_a_onium(i)=.true. ! charmonium
 
206
         if (idup(i,1,1).eq.9000443)   is_a_onium(i)=.true. ! charmonium
 
207
         if (idup(i,1,1).eq.9010443)   is_a_onium(i)=.true. ! charmonium
 
208
         if (idup(i,1,1).eq.9020443)   is_a_onium(i)=.true. ! charmonium
 
209
         if (idup(i,1,1).eq.445)   is_a_onium(i)=.true. ! charmonium
 
210
         if (idup(i,1,1).eq.9000445)   is_a_onium(i)=.true. ! charmonium
 
211
 
 
212
         if (idup(i,1,1).eq.551)   is_a_onium(i)=.true. ! bottomonium
 
213
         if (idup(i,1,1).eq.10551)   is_a_onium(i)=.true. ! bottomonium
 
214
         if (idup(i,1,1).eq.100551)   is_a_onium(i)=.true. ! bottomonium
 
215
         if (idup(i,1,1).eq.110551)   is_a_onium(i)=.true. ! bottomonium
 
216
         if (idup(i,1,1).eq.200551)   is_a_onium(i)=.true. ! bottomonium
 
217
         if (idup(i,1,1).eq.210551)   is_a_onium(i)=.true. ! bottomonium
 
218
         if (idup(i,1,1).eq.553)   is_a_onium(i)=.true. ! bottomonium
 
219
         if (idup(i,1,1).eq.10553)   is_a_onium(i)=.true. ! bottomonium
 
220
         if (idup(i,1,1).eq.20553)   is_a_onium(i)=.true. ! bottomonium
 
221
         if (idup(i,1,1).eq.30553)   is_a_onium(i)=.true. ! bottomonium
 
222
         if (idup(i,1,1).eq.100553)   is_a_onium(i)=.true. ! bottomonium
 
223
         if (idup(i,1,1).eq.110553)   is_a_onium(i)=.true. ! bottomonium
 
224
         if (idup(i,1,1).eq.120553)   is_a_onium(i)=.true. ! bottomonium
 
225
         if (idup(i,1,1).eq.130553)   is_a_onium(i)=.true. ! bottomonium
 
226
         if (idup(i,1,1).eq.200553)   is_a_onium(i)=.true. ! bottomonium
 
227
         if (idup(i,1,1).eq.210553)   is_a_onium(i)=.true. ! bottomonium
 
228
         if (idup(i,1,1).eq.220553)   is_a_onium(i)=.true. ! bottomonium
 
229
         if (idup(i,1,1).eq.300553)   is_a_onium(i)=.true. ! bottomonium
 
230
         if (idup(i,1,1).eq.9000553)   is_a_onium(i)=.true. ! bottomonium
 
231
         if (idup(i,1,1).eq.9010553)   is_a_onium(i)=.true. ! bottomonium
 
232
         if (idup(i,1,1).eq.555)   is_a_onium(i)=.true. ! bottomonium
 
233
         if (idup(i,1,1).eq.10555)   is_a_onium(i)=.true. ! bottomonium
 
234
         if (idup(i,1,1).eq.20555)   is_a_onium(i)=.true. ! bottomonium
 
235
         if (idup(i,1,1).eq.100555)   is_a_onium(i)=.true. ! bottomonium
 
236
         if (idup(i,1,1).eq.110555)   is_a_onium(i)=.true. ! bottomonium
 
237
         if (idup(i,1,1).eq.200555)   is_a_onium(i)=.true. ! bottomonium
 
238
         if (idup(i,1,1).eq.557)   is_a_onium(i)=.true. ! bottomonium
 
239
         if (idup(i,1,1).eq.100557)   is_a_onium(i)=.true. ! bottomonium
 
240
 
 
241
         if (idup(i,1,1).eq.541)   is_a_onium(i)=.true. ! Bc
 
242
         if (idup(i,1,1).eq.10541)   is_a_onium(i)=.true. ! Bc
 
243
         if (idup(i,1,1).eq.543)   is_a_onium(i)=.true. ! Bc
 
244
         if (idup(i,1,1).eq.10543)   is_a_onium(i)=.true. ! Bc
 
245
         if (idup(i,1,1).eq.20543)   is_a_onium(i)=.true. ! Bc
 
246
         if (idup(i,1,1).eq.545)   is_a_onium(i)=.true. ! Bc
 
247
      enddo
 
248
 
 
249
c
 
250
c     et and eta cuts
 
251
c
 
252
      do i=nincoming+1,nexternal
 
253
 
 
254
         etmin(i)  = 0d0
 
255
         etmax(i)  = -1
 
256
 
 
257
         emin(i)   = 0d0
 
258
         emax(i)   = -1
 
259
 
 
260
         etamin(i) = 0d0
 
261
         etamax(i) = -1
 
262
 
 
263
         if(do_cuts(i)) then
 
264
 
 
265
            if(is_a_j(i)) etmin(i)=ptj
 
266
            if(is_a_l(i)) etmin(i)=ptl
 
267
            if(is_a_b(i)) etmin(i)=ptb
 
268
            if(is_a_a(i)) etmin(i)=pta
 
269
            if(is_a_onium(i)) etmin(i)=ptonium
 
270
 
 
271
            if(is_a_j(i)) etmax(i)=ptjmax
 
272
            if(is_a_l(i)) etmax(i)=ptlmax
 
273
            if(is_a_b(i)) etmax(i)=ptbmax
 
274
            if(is_a_a(i)) etmax(i)=ptamax
 
275
 
 
276
c
 
277
            if(is_a_j(i)) emin(i)=ej
 
278
            if(is_a_l(i)) emin(i)=el
 
279
            if(is_a_b(i)) emin(i)=eb
 
280
            if(is_a_a(i)) emin(i)=ea
 
281
 
 
282
            if(is_a_j(i)) emax(i)=ejmax
 
283
            if(is_a_l(i)) emax(i)=elmax
 
284
            if(is_a_b(i)) emax(i)=ebmax
 
285
            if(is_a_a(i)) emax(i)=eamax
 
286
 
 
287
            if(is_a_j(i)) etamax(i)=etaj
 
288
            if(is_a_l(i)) etamax(i)=etal
 
289
            if(is_a_b(i)) etamax(i)=etab
 
290
            if(is_a_a(i)) etamax(i)=etaa
 
291
            if(is_a_onium(i)) etamax(i)=etaonium
 
292
 
 
293
            if(is_a_j(i)) etamin(i)=etajmin
 
294
            if(is_a_l(i)) etamin(i)=etalmin
 
295
            if(is_a_b(i)) etamin(i)=etabmin
 
296
            if(is_a_a(i)) etamin(i)=etaamin
 
297
         endif
 
298
      enddo
 
299
c
 
300
c     delta r cut
 
301
c
 
302
      do i=nincoming+1,nexternal-1
 
303
         do j=i+1,nexternal
 
304
            r2min(j,i)=0d0
 
305
            r2max(j,i)=-1
 
306
            if(do_cuts(i).and.do_cuts(j)) then
 
307
 
 
308
               if(is_a_j(i).and.is_a_j(j)) r2min(j,i)=drjj
 
309
               if(is_a_b(i).and.is_a_b(j)) r2min(j,i)=drbb
 
310
               if(is_a_l(i).and.is_a_l(j)) r2min(j,i)=drll
 
311
               if(is_a_a(i).and.is_a_a(j)) r2min(j,i)=draa
 
312
 
 
313
               if((is_a_b(i).and.is_a_j(j)).or.
 
314
     &           (is_a_j(i).and.is_a_b(j))) r2min(j,i)=drbj
 
315
               if((is_a_a(i).and.is_a_j(j)).or.
 
316
     &           (is_a_j(i).and.is_a_a(j))) r2min(j,i)=draj
 
317
               if((is_a_l(i).and.is_a_j(j)).or.
 
318
     &           (is_a_j(i).and.is_a_l(j))) r2min(j,i)=drjl
 
319
               if((is_a_b(i).and.is_a_a(j)).or.
 
320
     &           (is_a_a(i).and.is_a_b(j))) r2min(j,i)=drab
 
321
               if((is_a_b(i).and.is_a_l(j)).or.
 
322
     &           (is_a_l(i).and.is_a_b(j))) r2min(j,i)=drbl
 
323
               if((is_a_l(i).and.is_a_a(j)).or.
 
324
     &           (is_a_a(i).and.is_a_l(j))) r2min(j,i)=dral
 
325
 
 
326
               if(is_a_j(i).and.is_a_j(j)) r2max(j,i)=drjjmax
 
327
               if(is_a_b(i).and.is_a_b(j)) r2max(j,i)=drbbmax
 
328
               if(is_a_l(i).and.is_a_l(j)) r2max(j,i)=drllmax
 
329
               if(is_a_a(i).and.is_a_a(j)) r2max(j,i)=draamax
 
330
 
 
331
               if((is_a_b(i).and.is_a_j(j)).or.
 
332
     &           (is_a_j(i).and.is_a_b(j))) r2max(j,i)=drbjmax
 
333
               if((is_a_a(i).and.is_a_j(j)).or.
 
334
     &           (is_a_j(i).and.is_a_a(j))) r2max(j,i)=drajmax
 
335
               if((is_a_l(i).and.is_a_j(j)).or.
 
336
     &           (is_a_j(i).and.is_a_l(j))) r2max(j,i)=drjlmax
 
337
               if((is_a_b(i).and.is_a_a(j)).or.
 
338
     &           (is_a_a(i).and.is_a_b(j))) r2max(j,i)=drabmax
 
339
               if((is_a_b(i).and.is_a_l(j)).or.
 
340
     &           (is_a_l(i).and.is_a_b(j))) r2max(j,i)=drblmax
 
341
               if((is_a_l(i).and.is_a_a(j)).or.
 
342
     &           (is_a_a(i).and.is_a_l(j))) r2max(j,i)=dralmax
 
343
 
 
344
            endif
 
345
         enddo
 
346
      enddo
 
347
c     
 
348
c     smin cut
 
349
c
 
350
      do i=nincoming+1,nexternal-1
 
351
         do j=i+1,nexternal
 
352
            s_min(j,i)=0.0d0**2
 
353
            s_max(j,i)=-1
 
354
            if(do_cuts(i).and.do_cuts(j)) then
 
355
               if(is_a_j(i).and.is_a_j(j)) s_min(j,i)=mmjj*dabs(mmjj)   
 
356
               if(is_a_a(i).and.is_a_a(j)) s_min(j,i)=mmaa*dabs(mmaa)  
 
357
               if( is_a_b(i).and.is_a_b(j) ) s_min(j,i)=mmbb*dabs(mmbb)     
 
358
               if((is_a_l(i).and.is_a_l(j)).and.
 
359
     &            (abs(idup(i,1,1)).eq.abs(idup(j,1,1))).and.
 
360
     &            (idup(i,1,1)*idup(j,1,1).lt.0)) 
 
361
     &            s_min(j,i)=mmll*dabs(mmll)  !only on l+l- pairs (same flavour) 
 
362
 
 
363
               if(is_a_j(i).and.is_a_j(j)) s_max(j,i)=mmjjmax*dabs(mmjjmax)   
 
364
               if(is_a_a(i).and.is_a_a(j)) s_max(j,i)=mmaamax*dabs(mmaamax)  
 
365
               if( is_a_b(i).and.is_a_b(j) ) s_max(j,i)=mmbbmax*dabs(mmbbmax)     
 
366
               if((is_a_l(i).and.is_a_l(j)).and.
 
367
     &            (abs(idup(i,1,1)).eq.abs(idup(j,1,1))).and.
 
368
     &            (idup(i,1,1)*idup(j,1,1).lt.0)) 
 
369
     &            s_max(j,i)=mmllmax*dabs(mmllmax)  !only on l+l- pairs (same flavour)
 
370
            endif
 
371
         enddo
 
372
      enddo      
 
373
 
 
374
c     
 
375
c     ptll cut (min and max)
 
376
c
 
377
 
 
378
      do i=nincoming+1,nexternal-1
 
379
         do j=i+1,nexternal
 
380
            ptll_min(j,i)=0.0d0**2
 
381
            ptll_max(j,i)=-1
 
382
            if(((is_a_l(i).and.is_a_l(j)).and.
 
383
     &      (abs(idup(i,1,1)).eq.abs(idup(j,1,1))).and.
 
384
     &      (idup(i,1,1)*idup(j,1,1).lt.0))! Leptons from same flavor but different charge
 
385
     &       .or.(is_a_nu(i).and.is_a_l(j))
 
386
     &       .or.(is_a_l(i).and.is_a_nu(j))   !a lepton and a neutrino
 
387
     &       .or.(is_a_nu(i).and.is_a_nu(j))) then ! two neutrinos 
 
388
               ptll_min(j,i)=ptllmin*dabs(ptllmin)
 
389
               ptll_max(j,i)=ptllmax*dabs(ptllmax)
 
390
           endif
 
391
         enddo
 
392
      enddo
 
393
 
 
394
c
 
395
c   EXTRA JET CUTS
 
396
c
 
397
      ptjmin4(1)=ptj1min
 
398
      ptjmin4(2)=ptj2min
 
399
      ptjmin4(3)=ptj3min
 
400
      ptjmin4(4)=ptj4min
 
401
 
 
402
      ptjmax4(1)=ptj1max
 
403
      ptjmax4(2)=ptj2max
 
404
      ptjmax4(3)=ptj3max
 
405
      ptjmax4(4)=ptj4max
 
406
 
 
407
      Htjmin4(2)=ht2min
 
408
      htjmin4(3)=ht3min
 
409
      htjmin4(4)=ht4min
 
410
 
 
411
      htjmax4(2)=ht2max
 
412
      htjmax4(3)=ht3max
 
413
      htjmax4(4)=ht4max
 
414
   
 
415
      inclHtmin=ihtmin
 
416
      inclHtmax=ihtmax
 
417
 
 
418
      jetor = cutuse.eq.0d0
 
419
 
 
420
c
 
421
c   EXTRA LEPTON CUTS
 
422
c
 
423
      ptlmin4(1)=ptl1min
 
424
      ptlmin4(2)=ptl2min
 
425
      ptlmin4(3)=ptl3min
 
426
      ptlmin4(4)=ptl4min
 
427
 
 
428
      ptlmax4(1)=ptl1max
 
429
      ptlmax4(2)=ptl2max
 
430
      ptlmax4(3)=ptl3max
 
431
      ptlmax4(4)=ptl4max
 
432
 
 
433
c
 
434
c    ERROR TRAPS 
 
435
c
 
436
        do i=nincoming+1,nexternal
 
437
           if(is_a_j(i).and.etmin(i).eq.0.and.emin(i).eq.0) then
 
438
              write (*,*) "Warning: pt or E min of a jet should in general be >0"
 
439
           endif
 
440
           if(is_a_a(i).and.etmin(i).eq.0.and.emin(i).eq.0) then
 
441
              write (*,*) "Warning: pt or E min of a gamma should in general be >0"
 
442
           endif
 
443
        enddo
 
444
 
 
445
c    count number of jets to see if special cuts are applicable or not
 
446
 
 
447
        ncheck=0
 
448
        do i=nincoming+1,nexternal
 
449
           if(is_a_j(i)) ncheck=ncheck+1
 
450
        enddo
 
451
 
 
452
        if(ncheck.eq.0.and. xptj .gt. 0d0) then
 
453
           write (*,*) "Warning: cuts on the jet will be ignored"
 
454
           xptj = 0d0
 
455
        endif
 
456
 
 
457
        if(ncheck.lt.2.and. xetamin .gt. 0 .and. deltaeta .gt.0) then
 
458
           write (*,*) "Warning: WBF cuts not will be ignored"
 
459
           xetamin = 0d0
 
460
           deltaeta =0d0
 
461
        endif
 
462
 
 
463
c    count number of photons to see if special cuts are applicable or not
 
464
 
 
465
        ncheck=0
 
466
        do i=nincoming+1,nexternal
 
467
           if(is_a_a(i)) ncheck=ncheck+1
 
468
        enddo
 
469
 
 
470
        if(ncheck.eq.0.and. xpta .gt. 0d0) then
 
471
           write (*,*) "Warning: cuts on the photon will be ignored"
 
472
           xpta =0d0
 
473
        endif
 
474
 
 
475
c    count number of b-quarks to see if special cuts are applicable or not
 
476
 
 
477
        Ncheck=0
 
478
        do i=nincoming+1,nexternal
 
479
           if(is_a_b(i)) ncheck=ncheck+1
 
480
        enddo
 
481
 
 
482
        if(ncheck.eq.0.and. xptb .gt. 0d0) then
 
483
           write (*,*) "Warning: cuts on the b-quarks will be ignored"
 
484
           xptb = 0d0
 
485
        endif
 
486
 
 
487
c    count number of leptons to see if special cuts are applicable or not
 
488
 
 
489
        ncheck=0
 
490
        do i=nincoming+1,nexternal
 
491
           if(is_a_l(i)) ncheck=ncheck+1
 
492
        enddo
 
493
 
 
494
        if(ncheck.eq.0.and. xptl .gt. 0d0) then
 
495
           write (*,*) "Warning: cuts on the lepton will be ignored"
 
496
           xptl = 0d0
 
497
        endif
 
498
 
 
499
c     set possible xqcut combinations (for better grid preparation)
 
500
        if(xqcut.gt.0)
 
501
     $       call setxqcuts()
 
502
 
 
503
c      call write_cuts()
 
504
      RETURN
 
505
 
 
506
      END
 
507
 
 
508
 
 
509
      subroutine setxqcuts()
 
510
c**************************************************
 
511
c     Set xqcuti and xqcutij between all particles
 
512
c     to allow for grid preparation based on xqcut
 
513
c**************************************************
 
514
      implicit none
 
515
      include 'genps.inc'
 
516
      include 'maxconfigs.inc'
 
517
      include 'nexternal.inc'
 
518
      include 'cuts.inc'
 
519
 
 
520
      double precision pmass(nexternal)
 
521
      common/to_mass/  pmass
 
522
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
523
      common/to_forest/ iforest
 
524
      integer mapconfig(0:lmaxconfigs), this_config
 
525
      common/to_mconfigs/mapconfig, this_config
 
526
 
 
527
      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
 
528
      common/to_xqcuts/xqcutij,xqcuti
 
529
 
 
530
      
 
531
      integer i,j,k
 
532
      logical foundpartner
 
533
      include 'maxamps.inc'
 
534
      integer idup(nexternal,maxproc,maxsproc)
 
535
      integer mothup(2,nexternal)
 
536
      integer icolup(2,nexternal,maxflow,maxsproc)
 
537
      include 'leshouche.inc'
 
538
 
 
539
      do i=3,nexternal
 
540
         xqcuti(i)=0d0
 
541
         do j=3,nexternal
 
542
            xqcutij(i,j)=0d0
 
543
         enddo
 
544
      enddo
 
545
 
 
546
      do i=3,nexternal
 
547
         do j=1,mapconfig(0)
 
548
            foundpartner=.false.
 
549
            do k=-1,-(nexternal-3),-1
 
550
               if(iforest(1,k,j).eq.i.and.iforest(2,k,j).gt.2.or.
 
551
     $              iforest(2,k,j).eq.i.and.iforest(1,k,j).gt.2)then
 
552
                  foundpartner=.true.
 
553
                  if(iabs(idup(iforest(2,k,j),1,1)).le.maxjetflavor.or.
 
554
     $                 idup(iforest(2,k,j),1,1).eq.21.or.
 
555
     $                 iabs(idup(iforest(2,k,j),1,1)).le.maxjetflavor.or.
 
556
     $                 idup(iforest(2,k,j),1,1).eq.21)then
 
557
                     xqcutij(iforest(2,k,j),iforest(1,k,j))=xqcut
 
558
                     xqcutij(iforest(1,k,j),iforest(2,k,j))=xqcut
 
559
                  endif
 
560
               endif
 
561
            enddo
 
562
            if(.not.foundpartner.and.(iabs(idup(i,1,1)).le.maxjetflavor.or.
 
563
     $           idup(i,1,1).eq.21)) xqcuti(i)=max(0d0,sqrt(xqcut**2-pmass(i)**2))
 
564
         enddo
 
565
      enddo
 
566
 
 
567
      end
 
568
 
 
569
 
 
570
c************************************************************************
 
571
c     Subroutine to check which external particles are from decays
 
572
c************************************************************************
 
573
      subroutine check_decay(from_decay)
 
574
      implicit none
 
575
 
 
576
      include 'nexternal.inc'
 
577
      include 'maxconfigs.inc'
 
578
      include 'genps.inc'
 
579
      include 'maxamps.inc'
 
580
c
 
581
c    Arguments
 
582
c
 
583
      logical from_decay(-(nexternal+3):nexternal)
 
584
c
 
585
c     Global
 
586
c
 
587
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
588
      common/to_forest/ iforest
 
589
      integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
 
590
      integer tprid(-max_branch:-1,lmaxconfigs)
 
591
      common/to_sprop/sprop,tprid
 
592
      integer gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
 
593
      include 'decayBW.inc'
 
594
 
 
595
c
 
596
c     Local
 
597
c
 
598
      integer i,j
 
599
 
 
600
      do i=-(nexternal-3),nexternal
 
601
         from_decay(i)=.false.
 
602
      enddo
 
603
 
 
604
c     Set who comes from decay based on forced BW
 
605
      do i=-(nexternal-3),-1
 
606
         if(tprid(i,1).eq.0.and.gForceBW(i,1).eq.1.or.
 
607
     $        from_decay(i)) then
 
608
            from_decay(i)=.true.
 
609
            from_decay(iforest(1,i,1))=.true.
 
610
            from_decay(iforest(2,i,1))=.true.
 
611
         endif
 
612
      enddo
 
613
 
 
614
      end