← Back to branch summary

~alifson/chiralityflow/trunk

« back to all changes in this revision

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

  • Committer: andrew.lifson at lu
  • Date: 2021-09-01 15:34:39 UTC
  • Revision ID: andrew.lifson@thep.lu.se-20210901153439-7fasjhav4cp4m88r
testing a new repository of a madgraph folder

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/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
c
 
20
c     LOCAL
 
21
c
 
22
      integer i,j,k
 
23
C     
 
24
C     GLOBAL
 
25
C
 
26
c--masses and poles
 
27
      double precision pmass(nexternal)
 
28
      common/to_mass/  pmass
 
29
c
 
30
c     les houches accord stuff to identify particles
 
31
c
 
32
      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
 
33
     &     icolup(2,nexternal,maxflow),niprocs
 
34
      common /c_leshouche_inc/idup,mothup,icolup,niprocs
 
35
C
 
36
      LOGICAL  IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
 
37
      LOGICAL  IS_A_PH(NEXTERNAL)
 
38
      COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
 
39
c 
 
40
      double precision etmin(nincoming+1:nexternal-1)
 
41
      double precision etmax(nincoming+1:nexternal-1)
 
42
      double precision mxxmin(nincoming+1:nexternal-1,nincoming+1:nexternal-1)
 
43
      common /to_cuts/etmin,etmax, mxxmin
 
44
c
 
45
c     setup masses for the final-state particles (fills the /to_mass/ common block)
 
46
c
 
47
      include 'pmass.inc'
 
48
 
 
49
C-----
 
50
C  BEGIN CODE
 
51
C-----
 
52
c
 
53
c     No pdfs for decay processes - set here since here we know the nincoming
 
54
c     Also set stot here, and use mass of incoming particle for ren scale
 
55
c
 
56
         if(nincoming.eq.1)then
 
57
            lpp(1)=0
 
58
            lpp(2)=0
 
59
            ebeam(1)=pmass(1)/2d0
 
60
            ebeam(2)=pmass(1)/2d0
 
61
         endif
 
62
c-check consistency of maxjetflavor in the run_card and with Nf
 
63
c specified in coupl.inc
 
64
      if (maxjetflavor.lt.int(Nf)) then
 
65
         write(*,'(a,i3,a,i3)') "WARNING: the value of maxjetflavor"/
 
66
     $        /"specified in the run_card (",maxjetflavor,") is"/
 
67
     $        /" inconsistent with the number of light flavours in"/
 
68
     $        /"the model. Hence it will be set to:", int(Nf)
 
69
          maxjetflavor = int(Nf)
 
70
      endif
 
71
 
 
72
      do i=nincoming+1,nexternal
 
73
         is_a_j(i)=.false.
 
74
         is_a_lp(i)=.false.
 
75
         is_a_lm(i)=.false.
 
76
         is_a_ph(i)=.false.
 
77
         
 
78
c     -light-jets
 
79
         if (abs(idup(i,1)).le.maxjetflavor) then
 
80
              is_a_j(i)=.true.
 
81
         endif
 
82
         if (abs(idup(i,1)).eq.21)  is_a_j(i)=.true. ! gluon is a jet
 
83
 
 
84
c-charged-leptons
 
85
         if (idup(i,1).eq.11)  is_a_lm(i)=.true. !  e-
 
86
         if (idup(i,1).eq.13)  is_a_lm(i)=.true. !  mu-
 
87
         if (idup(i,1).eq.15)  is_a_lm(i)=.true. !  ta-
 
88
         if (idup(i,1).eq.-11) is_a_lp(i)=.true. !  e-
 
89
         if (idup(i,1).eq.-13) is_a_lp(i)=.true. !  mu-
 
90
         if (idup(i,1).eq.-15) is_a_lp(i)=.true. !  ta-
 
91
 
 
92
c-photons
 
93
         if (idup(i,1).eq.22.and..not.gamma_is_j)  is_a_ph(i)=.true. ! iso photon
 
94
         if (idup(i,1).eq.22.and.gamma_is_j)  is_a_j(i)=.true. !  photon in jets
 
95
      enddo
 
96
 
 
97
c
 
98
c     check for pdg specific cut (pt/eta)
 
99
c
 
100
      do i=nincoming+1, nexternal-1 ! never include last particle
 
101
         etmin(i) = 0d0
 
102
         etmax(i) = -1d0
 
103
         do j = i, nexternal-1
 
104
            mxxmin(i,j) = 0d0
 
105
         enddo
 
106
      enddo
 
107
 
 
108
      if (pdg_cut(1).ne.0)then
 
109
         do j=1,pdg_cut(0)
 
110
            do i=nincoming+1, nexternal-1 ! never include last particle
 
111
               if (abs(idup(i,1)).ne.pdg_cut(j)) cycle
 
112
c     fully ensure that only massive particles are allowed at NLO
 
113
               if(pmass(i).eq.0d0) then
 
114
                  write(*,*) 'Illegal use of pdg specific cut.'
 
115
                  write(*,*) 'For NLO process, '/
 
116
     $                 /'only massive particle can be included'
 
117
                  stop 1
 
118
               endif
 
119
c     fully ensure that this is not a jet/lepton/photon
 
120
               if(is_a_lp(i) .or. is_a_lm(i) .or. is_a_j(i) .or.
 
121
     $              is_a_ph(i)) then
 
122
                  write(*,*) 'Illegal use of pdg specific cut.'
 
123
                  write(*,*) 'This can not be used for '/
 
124
     $                 /'jet/lepton/photon/gluon'
 
125
                  stop 1
 
126
               endif
 
127
               etmin(i) = ptmin4pdg(j)
 
128
               etmax(i) = ptmax4pdg(j)
 
129
!     add the invariant mass cut
 
130
               if(mxxmin4pdg(j).ne.0d0)then
 
131
                  do k=i+1, nexternal-1
 
132
                     if (mxxpart_antipart(j))then
 
133
                        if (idup(k, 1).eq.-1*idup(i,1))then
 
134
                           mxxmin(i,k) = mxxmin4pdg(j)
 
135
                        endif
 
136
                     else
 
137
                        if (abs(idup(k, 1)).eq.pdg_cut(j))then
 
138
                           mxxmin(i,k) = mxxmin4pdg(j)
 
139
                        endif
 
140
                     endif
 
141
                  enddo
 
142
               endif
 
143
            enddo
 
144
         enddo
 
145
      endif
 
146
 
 
147
 
 
148
      RETURN
 
149
      END
 
150
 
 
151
 
 
152
      subroutine set_tau_min()
 
153
c Sets the lower bound for tau=x1*x2, using information on particle
 
154
c masses and on the jet minimum pt, as entered in run_card.dat, 
 
155
c variable ptj
 
156
      use mint_module
 
157
      implicit none
 
158
      double precision zero,vtiny
 
159
      parameter (zero=0.d0,vtiny=1d-8)
 
160
      include 'cuts.inc'
 
161
      include 'run.inc'
 
162
      include 'genps.inc'
 
163
      include 'nexternal.inc'
 
164
      include 'coupl.inc'
 
165
      include 'nFKSconfigs.inc'
 
166
      include "fks_info.inc"
 
167
      LOGICAL  IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
 
168
      LOGICAL  IS_A_PH(NEXTERNAL)
 
169
      COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
 
170
c
 
171
      integer pow(-nexternal:0,lmaxconfigs)
 
172
      double precision pmass(-nexternal:0,lmaxconfigs)
 
173
      double precision pwidth(-nexternal:0,lmaxconfigs)
 
174
      integer itree(2,-max_branch:-1),iconf
 
175
      common /to_itree/itree,iconf
 
176
      INTEGER NFKSPROCESS
 
177
      COMMON/C_NFKSPROCESS/NFKSPROCESS
 
178
      double precision taumin(fks_configs,maxchannels)
 
179
     $     ,taumin_s(fks_configs,maxchannels),taumin_j(fks_configs
 
180
     $     ,maxchannels),stot,xk(-nexternal:nexternal)
 
181
      save  taumin,taumin_s,taumin_j,stot
 
182
      integer i,j,k,d1,d2,iFKS,nt
 
183
      double precision xm(-nexternal:nexternal),xm1,xm2,xmi
 
184
      double precision xw(-nexternal:nexternal),xw1,xw2,xwi
 
185
      integer tsign,i_fks,j_fks
 
186
      double precision tau_Born_lower_bound,tau_lower_bound_resonance
 
187
     &     ,tau_lower_bound
 
188
      common/ctau_lower_bound/tau_Born_lower_bound
 
189
     &     ,tau_lower_bound_resonance,tau_lower_bound
 
190
c BW stuff
 
191
      double precision mass_min(-nexternal:nexternal,maxchannels)
 
192
     $     ,masslow(-nexternal:-1),widthlow(-nexternal:-1),sum_all_s
 
193
      save mass_min
 
194
      integer t_channel
 
195
      integer cBW_FKS_level_max(fks_configs,maxchannels),
 
196
     &     cBW_FKS(fks_configs,-nexternal:-1,maxchannels),
 
197
     &     cBW_FKS_level(fks_configs,-nexternal:-1,maxchannels)
 
198
      double precision cBW_FKS_mass(fks_configs,-1:1,-nexternal:-1
 
199
     $     ,maxchannels),cBW_FKS_width(fks_configs,-1:1,-nexternal:-1
 
200
     $     ,maxchannels)
 
201
      save cBW_FKS_level_max,cBW_FKS,cBW_FKS_level,cBW_FKS_mass
 
202
     $     ,cBW_FKS_width
 
203
      integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
 
204
      double precision cBW_mass(-1:1,-nexternal:-1),
 
205
     &     cBW_width(-1:1,-nexternal:-1)
 
206
      common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
 
207
     $     ,cBW_level
 
208
      double precision s_mass(-nexternal:nexternal)
 
209
     $     ,s_mass_FKS(fks_configs,-nexternal:nexternal,maxchannels)
 
210
      save s_mass_FKS
 
211
      common/to_phase_space_s_channel/s_mass
 
212
c Les Houches common block
 
213
      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
 
214
     &     icolup(2,nexternal,maxflow),niprocs
 
215
      common /c_leshouche_inc/idup,mothup,icolup,niprocs
 
216
      real*8         emass(nexternal)
 
217
      common/to_mass/emass
 
218
c     block for the (simple) cut bsed on the pdg
 
219
      double precision etmin(nincoming+1:nexternal-1)
 
220
      double precision etmax(nincoming+1:nexternal-1)
 
221
      double precision mxxmin(nincoming+1:nexternal-1,nincoming+1:nexternal-1)
 
222
      common /to_cuts/etmin,etmax,mxxmin
 
223
      double precision smin_update , mxx
 
224
      integer nb_iden_pdg
 
225
c
 
226
      logical firsttime,firsttime_chans(maxchannels)
 
227
      data firsttime /.true./
 
228
      data firsttime_chans/maxchannels*.true./
 
229
      if (firsttime) then
 
230
         do i = 1,lmaxconfigs
 
231
            do j = -nexternal,0
 
232
               pmass(j,i) = 0d0
 
233
               pwidth(j,i) = 0d0
 
234
            enddo
 
235
         enddo
 
236
         firsttime=.false.
 
237
      endif
 
238
      include "born_props.inc"
 
239
 
 
240
c The following assumes that light QCD particles are at the end of the
 
241
c list. Exclude one of them (i_fks) to set tau bound at the Born level 
 
242
c This sets a hard cut in the minimal shat of the Born phase-space
 
243
c generation.
 
244
c
 
245
c The contribution from ptj should be treated only as a 'soft lower
 
246
c bound' if j_fks is initial state: the real-emission i_fks parton is
 
247
c not necessarily the softest.  Therefore, it could be that even though
 
248
c the Born does not have enough energy to pass the cuts set by ptj, the
 
249
c event could.
 
250
      if (firsttime_chans(ichan)) then
 
251
         do i=-nexternal,nexternal
 
252
            xm(i)=0d0
 
253
            xw(i)=0d0
 
254
            mass_min(i,ichan)=0d0
 
255
         end do
 
256
         firsttime_chans(ichan)=.false.
 
257
         do iFKS=1,fks_configs
 
258
            j_fks=FKS_J_D(iFKS)
 
259
            i_fks=FKS_I_D(iFKS)
 
260
            taumin(iFKS,ichan)=0.d0
 
261
            taumin_s(iFKS,ichan)=0.d0
 
262
            taumin_j(iFKS,ichan)=0.d0
 
263
            do i=nincoming+1,nexternal
 
264
C Skip i_fks
 
265
               if (i.eq.i_fks) cycle
 
266
c Add the minimal jet pTs to tau
 
267
               if(IS_A_J(i)) then
 
268
                  if  (j_fks.gt.nincoming .and. j_fks.lt.nexternal) then
 
269
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
 
270
                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
 
271
                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
 
272
                  elseif (j_fks.ge.1 .and. j_fks.le.nincoming) then
 
273
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
 
274
                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
 
275
                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
 
276
                  elseif (j_fks.eq.nexternal) then
 
277
                     write (*,*)
 
278
     &                    'ERROR, j_fks cannot be the final parton'
 
279
     &                    ,j_fks
 
280
                     stop
 
281
                  else
 
282
                     write (*,*) 'ERROR, j_fks not correctly defined'
 
283
     &                    ,j_fks
 
284
                     stop
 
285
                  endif
 
286
                  xm(i)=emass(i)+ptj
 
287
c Add the minimal photon pTs to tau
 
288
               elseif(IS_A_PH(i))then
 
289
                  if (abs(emass(i)).gt.vtiny) then
 
290
                     write (*,*) 'Error in set_tau_min in setcuts.f:'
 
291
                     write (*,*) 'mass of a photon should be zero',i
 
292
     &                    ,emass(i)
 
293
                     stop
 
294
                  endif
 
295
                  if  (j_fks.gt.nincoming)
 
296
     &                 taumin(iFKS,ichan)=taumin(iFKS,ichan)+ptgmin
 
297
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+ptgmin
 
298
                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+ptgmin
 
299
                  xm(i)=emass(i)+ptgmin
 
300
               elseif (is_a_lp(i)) then
 
301
c Add the postively charged lepton pTs to tau
 
302
                  if (j_fks.gt.nincoming) then
 
303
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
 
304
                  else
 
305
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
 
306
                  endif
 
307
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(emass(i)**2+ptl**2)
 
308
                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(emass(i)**2+ptl**2)
 
309
                  xm(i)=emass(i)+ptl
 
310
c Add the lepton invariant mass to tau if there is at least another
 
311
c lepton of opposite charge. (Only add half of it, i.e. 'the part
 
312
c contributing from this lepton'). Remove possible overcounting with the
 
313
c lepton pT
 
314
                  do j=nincoming+1,nexternal
 
315
                     if (is_a_lm(j) .and. idup(i,1).eq.-idup(j,1) .and.
 
316
     $                    (mll_sf.ne.0d0 .or. mll.ne.0d0) ) then
 
317
                        if (j_fks.gt.nincoming)
 
318
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
 
319
     &                              max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
320
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
321
     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
322
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
323
     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
324
                        xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,mll_sf/2d0
 
325
     $                       ,ptl+emass(i))
 
326
                        exit
 
327
                     elseif (is_a_lm(j) .and. mll.ne.0d0) then
 
328
                        if (j_fks.gt.nincoming)
 
329
     &                       taumin(iFKS,ichan)= taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
 
330
     &                                     max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
 
331
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
332
     $                       + max(mll/2d0, dsqrt(ptl**2+emass(i)**2))
 
333
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
334
     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
 
335
                        xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
 
336
     $                       +emass(i))
 
337
                        exit
 
338
                     endif
 
339
                  enddo
 
340
               elseif (is_a_lm(i)) then
 
341
c Add the negatively charged lepton pTs to tau
 
342
                  if (j_fks.gt.nincoming) then
 
343
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
 
344
                  else
 
345
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
 
346
                  endif
 
347
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
 
348
                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
 
349
                  xm(i)=emass(i)+ptl
 
350
c Add the lepton invariant mass to tau if there is at least another
 
351
c lepton of opposite charge. (Only add half of it, i.e. 'the part
 
352
c contributing from this lepton'). Remove possible overcounting with the
 
353
c lepton pT
 
354
                  do j=nincoming+1,nexternal
 
355
                     if (is_a_lp(j) .and. idup(i,1).eq.-idup(j,1) .and.
 
356
     $                    (mll_sf.ne.0d0 .or. mll.ne.0d0) ) then
 
357
                        if (j_fks.gt.nincoming)
 
358
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
 
359
     &                              max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
360
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
361
     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
362
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
363
     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
364
                        xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,mll_sf/2d0
 
365
     $                       ,ptl+emass(i))
 
366
                        exit
 
367
                     elseif (is_a_lp(j) .and. mll.ne.0d0) then
 
368
                        if (j_fks.gt.nincoming)
 
369
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
 
370
     &                                      max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
 
371
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
372
     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
 
373
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
374
     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
 
375
                        xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
 
376
     $                       +emass(i))
 
377
                        exit
 
378
                     endif
 
379
                  enddo
 
380
               else
 
381
                  if (i.eq.nexternal)then
 
382
                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + emass(i)
 
383
                        taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan) +  emass(i)
 
384
                        taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan) + emass(i)
 
385
                        xm(i) = emass(i)
 
386
                  else
 
387
                     smin_update = 0
 
388
                     nb_iden_pdg = 1
 
389
                     mxx = 0d0
 
390
c                    assume smin apply always on the same set of particle
 
391
                     do j=nincoming+1,nexternal-1
 
392
                        if (mxxmin(i,j).ne.0d0.or.mxxmin(j,i).ne.0d0) then
 
393
                           nb_iden_pdg = nb_iden_pdg +1
 
394
                           if (mxx.eq.0d0) mxx = max(mxxmin(i,j), mxxmin(j,i))
 
395
                        endif
 
396
                     enddo
 
397
                     ! S >= (2*N-N^2)*M1^2 + (N^2-N)/2 * Mxx^2
 
398
                     smin_update = nb_iden_pdg*((2-nb_iden_pdg)*emass(i)**2 + (nb_iden_pdg-1)/2.*mxx**2)
 
399
                     ! compare with the update from pt cut
 
400
                     if (smin_update.lt.nb_iden_pdg**2*(etmin(i)**2 + emass(i)**2))then
 
401
                        ! the pt is more restrictive
 
402
                        smin_update = dsqrt(etmin(i)**2 + emass(i)**2)
 
403
                     else
 
404
                        smin_update = dsqrt(smin_update)/nb_iden_pdg ! share over N particle, and change dimension
 
405
                     endif
 
406
                     ! update in sqrt(s) so take the 
 
407
                     if  (j_fks.gt.nincoming) then
 
408
                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + smin_update
 
409
                     else
 
410
                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + emass(i)
 
411
                     endif
 
412
                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan) + smin_update
 
413
                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan) + smin_update
 
414
                     xm(i) = smin_update
 
415
                  endif
 
416
               endif
 
417
               xw(i)=0d0
 
418
            enddo
 
419
            stot = 4d0*ebeam(1)*ebeam(2)
 
420
            tau_Born_lower_bound=taumin(iFKS,ichan)**2/stot
 
421
            tau_lower_bound=taumin_j(iFKS,ichan)**2/stot
 
422
c         
 
423
c Also find the minimum lower bound if all internal s-channel particles
 
424
c were on-shell
 
425
            tsign=-1
 
426
            nt=0
 
427
            do i=-1,-(nexternal-3),-1 ! All propagators
 
428
               if ( itree(1,i) .eq. 1 .or. itree(1,i) .eq. 2 ) tsign=1
 
429
               if (tsign.eq.-1) then ! s-channels
 
430
                  d1=itree(1,i)
 
431
                  d2=itree(2,i)
 
432
c If daughter is a jet, we should treat the ptj as a mass. Except if
 
433
c d1=nexternal, because we check the Born, so final parton should be
 
434
c skipped. [This is already done above; also for the leptons]
 
435
                  xm1=xm(d1)
 
436
                  xm2=xm(d2)
 
437
                  xw1=xw(d1)
 
438
                  xw2=xw(d2)
 
439
c On-shell mass of the intermediate resonance
 
440
                  xmi=pmass(i,iconf)
 
441
c Width of the intermediate resonance
 
442
                  xwi=pwidth(i,iconf)
 
443
c Set the intermediate mass equal to the max of its actual mass and
 
444
c the sum of the masses of the two daugters.
 
445
                  if (xmi.gt.xm1+xm2) then
 
446
                     xm(i)=xmi
 
447
                     xw(i)=xwi
 
448
                  else
 
449
                     xm(i)=xm1+xm2
 
450
                     xw(i)=xw1+xw2 ! just sum the widths
 
451
                  endif
 
452
c Add the new mass to the bound. To avoid double counting, we should
 
453
c subtract the daughters, because they are already included above or in
 
454
c the previous iteration of the loop
 
455
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+xm(i)-xm1-xm2
 
456
               else             ! t-channels
 
457
                  if (i.eq.-(nexternal-3)) then
 
458
                     xm(i)=0d0
 
459
                     cycle
 
460
                  endif
 
461
                  nt=nt+1
 
462
                  d1=itree(2,i) ! only use 2nd daughter (which is the outgoing one)
 
463
                  xm1=xm(d1)
 
464
                  if (nt.gt.1) xm1=max(xm1,xk(nt-1))
 
465
                  xk(nt)=xm1
 
466
                  j=i-1         ! this is the closest to p2
 
467
                  d2=itree(2,j)
 
468
                  xm2=xm(d2)
 
469
                  xm(i)=min(xm1,xm2)
 
470
               endif
 
471
            enddo
 
472
      
 
473
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
474
c Determine the "minimal" s-channel invariant masses
 
475
            do i=nincoming+1,nexternal-1
 
476
               s_mass_FKS(iFKS,i,ichan)=xm(i)**2
 
477
            enddo
 
478
            do i=-1,-(nexternal-3),-1 ! All propagators 
 
479
               s_mass_FKS(iFKS,i,ichan)=xm(i)**2
 
480
            enddo
 
481
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
482
c Determine the conflicting Breit-Wigner's. Note that xm(i) contains the
 
483
c mass of the BW
 
484
            do i=nincoming+1,nexternal-1
 
485
               mass_min(i,ichan)=xm(i) ! minimal allowed resonance mass (including masses set by cuts)
 
486
            enddo
 
487
            cBW_FKS_level_max(iFKS,ichan)=0
 
488
            t_channel=0
 
489
            do i=-1,-(nexternal-3),-1 ! All propagators
 
490
               cBW_FKS_mass(iFKS,1,i,ichan)=0d0
 
491
               cBW_FKS_width(iFKS,1,i,ichan)=0d0
 
492
               cBW_FKS_mass(iFKS,-1,i,ichan)=0d0
 
493
               cBW_FKS_width(iFKS,-1,i,ichan)=0d0
 
494
               masslow(i)=9d99
 
495
               widthlow(i)=0d0
 
496
               if ( itree(1,i).eq.1 .or. itree(1,i).eq.2 ) t_channel=i
 
497
               if (t_channel.ne.0) exit ! only s-channels
 
498
               mass_min(i,ichan)=mass_min(itree(1,i),ichan)
 
499
     $              +mass_min(itree(2,i),ichan)
 
500
               if (xm(i).lt.mass_min(i,ichan)-vtiny) then
 
501
                  write (*,*)
 
502
     $                 'ERROR in the determination of conflicting BW',i
 
503
     $                 ,xm(i),mass_min(i,ichan)
 
504
                  stop
 
505
               endif
 
506
               if (pmass(i,iconf).lt.xm(i) .and.
 
507
     $              pwidth(i,iconf).gt.0d0) then
 
508
c     Possible conflict in BW
 
509
                  if (pmass(i,iconf).lt.mass_min(i,ichan)) then
 
510
c     Resonance can never go on-shell due to the kinematics of the event
 
511
                     cBW_FKS(iFKS,i,ichan)=2
 
512
                     cBW_FKS_level(iFKS,i,ichan)=0
 
513
                  elseif(pmass(i,iconf).lt.xm(i)) then
 
514
c     Conflicting Breit-Wigner
 
515
                     cBW_FKS(iFKS,i,ichan)=1
 
516
                     cBW_FKS_level(iFKS,i,ichan)=1
 
517
                     cBW_FKS_level_max(iFKS,ichan)=max(cBW_FKS_level_max(iFKS,ichan)
 
518
     $                    ,cBW_FKS_level(iFKS,i,ichan))
 
519
c     Set here the mass (and width) of the alternative mass; it's the
 
520
c     sum of daughter masses. (2nd argument is '1', because this
 
521
c     alternative mass is LARGER than the resonance mass).
 
522
                     cBW_FKS_mass(iFKS,1,i,ichan)=xm(i)
 
523
                     cBW_FKS_width(iFKS,1,i,ichan)=xw(i)
 
524
                  endif
 
525
c     set the daughters also as conflicting (recursively)
 
526
                  masslow(i)=pmass(i,iconf)
 
527
                  widthlow(i)=pwidth(i,iconf)
 
528
                  do j=i,-1
 
529
                     if (cBW_FKS(iFKS,j,ichan).eq.0) cycle
 
530
                     do k=1,2   ! loop over the 2 daughters
 
531
                        if (itree(k,j).ge.0) cycle
 
532
                        if (cBW_FKS(iFKS,itree(k,j),ichan).eq.2) cycle
 
533
                        cBW_FKS(iFKS,itree(k,j),ichan)=1
 
534
                        cBW_FKS_level(iFKS,itree(k,j),ichan)=
 
535
     $                       cBW_FKS_level(iFKS,j,ichan)+1
 
536
                        cBW_FKS_level_max(iFKS,ichan)=
 
537
     $                       max(cBW_FKS_level_max(iFKS,ichan)
 
538
     $                       ,cBW_FKS_level(iFKS,itree(k,j),ichan))
 
539
c     Set here the mass (and width) of the alternative mass; it's the
 
540
c     difference between the mother and the sister masses. (3rd argument
 
541
c     is '-1', because this alternative mass is SMALLER than the
 
542
c     resonance mass).
 
543
                        masslow(itree(k,j))=min(masslow(itree(k,j)),
 
544
     &                       max(masslow(j)-xm(itree(3-k,j)),0d0)) ! mass difference
 
545
                        widthlow(itree(k,j))=max(widthlow(itree(k,j)),
 
546
     &                       widthlow(j)+xw(itree(3-k,j))) ! sum of widths
 
547
                        if (pwidth(itree(k,j),iconf).eq.0d0 .or.
 
548
     $                       masslow(itree(k,j)).ge.pmass(itree(k,j)
 
549
     $                       ,iconf)) cycle
 
550
                        cBW_FKS_mass(iFKS,-1,itree(k,j),ichan)=
 
551
     $                       masslow(itree(k,j))
 
552
                        cBW_FKS_width(iFKS,-1,itree(k,j),ichan)=
 
553
     $                       widthlow(itree(k,j))
 
554
                     enddo
 
555
                  enddo
 
556
               else
 
557
c     Normal Breit-Wigner
 
558
                  cBW_FKS(iFKS,i,ichan)=0
 
559
               endif
 
560
            enddo
 
561
c loop over t-channel to make sure that s-hat is consistent with sum of
 
562
c s-channel masses
 
563
            if (t_channel.ne.0) then
 
564
               sum_all_s=0d0
 
565
               do i=t_channel,-(nexternal-3),-1
 
566
c Breit-wigner can never go on-shell:
 
567
                  if (itree(2,i).gt.0) cycle
 
568
                  if ( pmass(itree(2,i),iconf).gt.sqrt(stot) .and.
 
569
     $                 pwidth(itree(2,i),iconf).gt.0d0) then
 
570
                     cBW_FKS(iFKS,itree(2,i),ichan)=2
 
571
                  endif
 
572
c     s-channel is always 2nd argument of itree, sum it to sum_all_s
 
573
                  sum_all_s=sum_all_s+xm(itree(2,i))
 
574
               enddo
 
575
               if (sum_all_s.gt.sqrt(stot)) then
 
576
c     conflicting BWs: set all s-channels as conflicting
 
577
                  do i=t_channel,-(nexternal-3),-1
 
578
                     if (itree(2,i).gt.0) cycle
 
579
                     if (cBW_FKS(iFKS,itree(2,i),ichan).ne.2) then
 
580
                        cBW_FKS(iFKS,itree(2,i),ichan)=1
 
581
                        cBW_FKS_mass(iFKS,-1,itree(2,i),ichan)=sqrt(stot)/2d0
 
582
                        cBW_FKS_width(iFKS,-1,itree(2,i),ichan)=xw(itree(2,i))
 
583
                     endif
 
584
                  enddo
 
585
               endif
 
586
            endif
 
587
 
 
588
 
 
589
c Conflicting BW's determined. They are saved in cBW_FKS
 
590
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
591
c
 
592
c If the lower bound found here is smaller than the hard bound,
 
593
c simply set the soft bound equal to the hard bound.
 
594
            taumin_s(iFKS,ichan)=
 
595
     &           max(taumin_j(iFKS,ichan),taumin_s(iFKS,ichan))
 
596
c
 
597
c For the bound, we have to square and divide by stot.
 
598
            tau_lower_bound_resonance=taumin_s(iFKS,ichan)**2/stot
 
599
c
 
600
            if (j_fks.gt.nincoming) then
 
601
               write (*,'(a7,x,i3,x,i5,x,a1,3(e12.5,x)))') 'tau_min'
 
602
     $              ,iFKS,ichan,':',taumin(iFKS,ichan),taumin_j(iFKS
 
603
     $              ,ichan),taumin_s(iFKS,ichan)
 
604
            else
 
605
               write (*,'(a7,x,i3,x,i5,x,a1,e12.5,x,a13,e12.5,x))')
 
606
     $              'tau_min',iFKS,ichan,':',taumin(iFKS,ichan)
 
607
     $              ,'     --      ',taumin_s(iFKS,ichan)
 
608
            endif
 
609
         enddo
 
610
      endif
 
611
      tau_Born_lower_bound=taumin(nFKSprocess,ichan)**2/stot
 
612
      tau_lower_bound=taumin_j(nFKSprocess,ichan)**2/stot
 
613
      tau_lower_bound_resonance=taumin_s(nFKSprocess,ichan)**2/stot
 
614
      do i=-nexternal,-1
 
615
         cBW(i)=cBW_FKS(nFKSprocess,i,ichan)
 
616
         cBW_level(i)=cBW_FKS_level(nFKSprocess,i,ichan)
 
617
         do j=-1,1,2
 
618
            cBW_mass(j,i)=cBW_FKS_mass(nFKSprocess,j,i,ichan)
 
619
            cBW_width(j,i)=cBW_FKS_width(nFKSprocess,j,i,ichan)
 
620
         enddo
 
621
      enddo
 
622
      do i=-nexternal,nexternal
 
623
         s_mass(i)=s_mass_FKS(nFKSprocess,i,ichan)
 
624
      enddo
 
625
      cBW_level_max=cBW_FKS_level_max(nFKSprocess,ichan)
 
626
      call set_granny(nFKSprocess,iconf,mass_min(-nexternal,ichan))
 
627
      return
 
628
      end
 
629
 
 
630
 
 
631
      subroutine sChan_order(ns_channel,order)
 
632
      use mint_module
 
633
      implicit none
 
634
      include 'nexternal.inc'
 
635
      include 'maxparticles.inc'
 
636
      include 'maxconfigs.inc'
 
637
      integer itree(2,-max_branch:-1),iconf
 
638
      common /to_itree/itree,iconf
 
639
      double precision ran2,rnd
 
640
      integer i,j,order(-nexternal:0),ipos,ns_channel,npos
 
641
     $     ,pos(nexternal),ord(-nexternal:0)
 
642
      logical done(-nexternal:nexternal)
 
643
      external ran2
 
644
      save ord
 
645
      if (.not. new_point) then
 
646
         do j=-1,-ns_channel,-1
 
647
            order(j)=ord(j)
 
648
         enddo
 
649
         return
 
650
      endif
 
651
      do i=-ns_channel,0
 
652
         done(i)=.false.
 
653
      enddo
 
654
      do i=1,nexternal
 
655
         done(i)=.true.
 
656
      enddo
 
657
      do j=-1,-ns_channel,-1
 
658
         npos=0
 
659
         do i=-1,-ns_channel,-1
 
660
            if((.not. done(i)) .and.
 
661
     &           done(itree(1,i))  .and. done(itree(2,i))) then
 
662
               npos=npos+1
 
663
               pos(npos)=i
 
664
            endif
 
665
         enddo
 
666
         if (npos.gt.1) then
 
667
            rnd=ran2()
 
668
            ipos=min(int(rnd*npos)+1,npos)
 
669
            ord(j)=pos(ipos)
 
670
            done(pos(ipos))=.true.
 
671
         elseif (npos.eq.1) then
 
672
            ord(j)=pos(npos)
 
673
            done(pos(npos))=.true.
 
674
         else
 
675
            write (*,*) 'ERROR in sChan_order',npos
 
676
         endif
 
677
         order(j)=ord(j)
 
678
      enddo
 
679
      new_point=.false.
 
680
      return
 
681
      end
 
682
 
 
683