← Back to branch summary

~maddevelopers/mg5amcnlo/2.7.1.3

« back to all changes in this revision

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

  • Committer: olivier Mattelaer
  • Date: 2017-12-12 21:11:08 UTC
  • mfrom: (274.1.52 2.6.1)
  • Revision ID: olivier.mattelaer@uclouvain.be-20171212211108-xfh9rt0j9c9m6nxp
pass to 2.6.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/setcuts.f
19
19
c
20
20
c     LOCAL
21
21
c
22
 
      integer i,j
 
22
      integer i,j,k
23
23
C     
24
24
C     GLOBAL
25
25
C
36
36
      LOGICAL  IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
37
37
      LOGICAL  IS_A_PH(NEXTERNAL)
38
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
39
44
c
40
45
c     setup masses for the final-state particles (fills the /to_mass/ common block)
41
46
c
87
92
c-photons
88
93
         if (idup(i,1).eq.22)  is_a_ph(i)=.true. !  photon
89
94
      enddo
 
95
c
 
96
c     check for pdg specific cut (pt/eta)
 
97
c
 
98
      do i=nincoming+1, nexternal-1 ! never include last particle
 
99
         etmin(i) = 0d0
 
100
         etmax(i) = -1d0
 
101
         do j = i, nexternal-1
 
102
            mxxmin(i,j) = 0d0
 
103
         enddo
 
104
      enddo
 
105
 
 
106
      if (pdg_cut(1).ne.0)then
 
107
         do j=1,pdg_cut(0)
 
108
            do i=nincoming+1, nexternal-1 ! never include last particle
 
109
               if (abs(idup(i,1)).ne.pdg_cut(j)) cycle
 
110
c     fully ensure that only massive particles are allowed at NLO
 
111
               if(pmass(i).eq.0d0) then
 
112
                  write(*,*) 'Illegal use of pdg specific cut.'
 
113
                  write(*,*) 'For NLO process, '/
 
114
     $                 /'only massive particle can be included'
 
115
                  stop 1
 
116
               endif
 
117
c     fully ensure that this is not a jet/lepton/photon
 
118
               if(is_a_lp(i) .or. is_a_lm(i) .or. is_a_j(i) .or.
 
119
     $              is_a_ph(i)) then
 
120
                  write(*,*) 'Illegal use of pdg specific cut.'
 
121
                  write(*,*) 'This can not be used for '/
 
122
     $                 /'jet/lepton/photon/gluon'
 
123
                  stop 1
 
124
               endif
 
125
               etmin(i) = ptmin4pdg(j)
 
126
               etmax(i) = ptmax4pdg(j)
 
127
!     add the invariant mass cut
 
128
               if(mxxmin4pdg(j).ne.0d0)then
 
129
                  do k=i+1, nexternal-1
 
130
                     if (mxxpart_antipart(j))then
 
131
                        if (idup(k, 1).eq.-1*idup(i,1))then
 
132
                           mxxmin(i,k) = mxxmin4pdg(j)
 
133
                        endif
 
134
                     else
 
135
                        if (abs(idup(k, 1)).eq.pdg_cut(j))then
 
136
                           mxxmin(i,k) = mxxmin4pdg(j)
 
137
                        endif
 
138
                     endif
 
139
                  enddo
 
140
               endif
 
141
            enddo
 
142
         enddo
 
143
      endif
 
144
 
90
145
 
91
146
      RETURN
92
147
      END
157
212
      common /c_leshouche_inc/idup,mothup,icolup,niprocs
158
213
      real*8         emass(nexternal)
159
214
      common/to_mass/emass
 
215
c     block for the (simple) cut bsed on the pdg
 
216
      double precision etmin(nincoming+1:nexternal-1)
 
217
      double precision etmax(nincoming+1:nexternal-1)
 
218
      double precision mxxmin(nincoming+1:nexternal-1,nincoming+1:nexternal-1)
 
219
      common /to_cuts/etmin,etmax,mxxmin
 
220
      double precision smin_update , mxx
 
221
      integer nb_iden_pdg
 
222
c
160
223
      logical firsttime,firsttime_chans(maxchannels)
161
224
      data firsttime /.true./
162
225
      data firsttime_chans/maxchannels*.true./
201
264
c Add the minimal jet pTs to tau
202
265
               if(IS_A_J(i) .and. i.ne.nexternal)then
203
266
                  if  (j_fks.gt.nincoming .and. j_fks.lt.nexternal) then
204
 
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+max(ptj,emass(i))
205
 
                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+max(ptj,emass(i))
206
 
                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+max(ptj,emass(i))
 
267
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
 
268
                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
 
269
                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
207
270
                  elseif (j_fks.ge.1 .and. j_fks.le.nincoming) then
208
271
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
209
 
                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+max(ptj,emass(i))
210
 
                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+max(ptj,emass(i))
 
272
                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
 
273
                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptj**2 +emass(i)**2)
211
274
                  elseif (j_fks.eq.nexternal) then
212
275
                     write (*,*)
213
276
     &                    'ERROR, j_fks cannot be the final parton'
234
297
                  xm(i)=emass(i)+ptgmin
235
298
               elseif (is_a_lp(i)) then
236
299
c Add the postively charged lepton pTs to tau
237
 
                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
238
 
                  if (j_fks.gt.nincoming)
239
 
     &                 taumin(iFKS,ichan)=taumin(iFKS,ichan)+ptl
240
 
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)+ptl
241
 
                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)+ptl
 
300
                  if (j_fks.gt.nincoming) then
 
301
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
 
302
                  else
 
303
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
 
304
                  endif
 
305
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(emass(i)**2+ptl**2)
 
306
                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(emass(i)**2+ptl**2)
242
307
                  xm(i)=emass(i)+ptl
243
308
c Add the lepton invariant mass to tau if there is at least another
244
309
c lepton of opposite charge. (Only add half of it, i.e. 'the part
248
313
                     if (is_a_lm(j) .and. idup(i,1).eq.-idup(j,1) .and.
249
314
     $                    (mll_sf.ne.0d0 .or. mll.ne.0d0) ) then
250
315
                        if (j_fks.gt.nincoming)
251
 
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
252
 
     &                              max(mll/2d0,mll_sf/2d0,ptl+emass(i))
253
 
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
254
 
     $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
255
 
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
256
 
     $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
 
316
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
 
317
     &                              max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
318
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
319
     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
320
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
321
     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
257
322
                        xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,mll_sf/2d0
258
323
     $                       ,ptl+emass(i))
259
324
                        exit
260
325
                     elseif (is_a_lm(j) .and. mll.ne.0d0) then
261
326
                        if (j_fks.gt.nincoming)
262
 
     &                       taumin(iFKS,ichan)= taumin(iFKS,ichan)-ptl-emass(i) +
263
 
     &                                     max(mll/2d0,ptl+emass(i))
264
 
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
265
 
     $                       + max(mll/2d0,ptl+emass(i))
266
 
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
267
 
     $                       + max(mll/2d0,ptl+emass(i))
 
327
     &                       taumin(iFKS,ichan)= taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
 
328
     &                                     max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
 
329
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
330
     $                       + max(mll/2d0, dsqrt(ptl**2+emass(i)**2))
 
331
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
332
     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
268
333
                        xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
269
334
     $                       +emass(i))
270
335
                        exit
272
337
                  enddo
273
338
               elseif (is_a_lm(i)) then
274
339
c Add the negatively charged lepton pTs to tau
275
 
                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
276
 
                  if (j_fks.gt.nincoming)
277
 
     &                 taumin(iFKS,ichan)=taumin(iFKS,ichan)+ptl
278
 
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)+ptl
279
 
                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)+ptl
 
340
                  if (j_fks.gt.nincoming) then
 
341
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
 
342
                  else
 
343
                     taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
 
344
                  endif
 
345
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
 
346
                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+dsqrt(ptl**2+emass(i)**2)
280
347
                  xm(i)=emass(i)+ptl
281
348
c Add the lepton invariant mass to tau if there is at least another
282
349
c lepton of opposite charge. (Only add half of it, i.e. 'the part
286
353
                     if (is_a_lp(j) .and. idup(i,1).eq.-idup(j,1) .and.
287
354
     $                    (mll_sf.ne.0d0 .or. mll.ne.0d0) ) then
288
355
                        if (j_fks.gt.nincoming)
289
 
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
290
 
     &                              max(mll/2d0,mll_sf/2d0,ptl+emass(i))
291
 
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
292
 
     $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
293
 
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
294
 
     $                       + max(mll/2d0,mll_sf/2d0,ptl+emass(i))
 
356
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
 
357
     &                              max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
358
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
359
     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
 
360
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
361
     $                       + max(mll/2d0,mll_sf/2d0,dsqrt(ptl**2+emass(i)**2))
295
362
                        xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,mll_sf/2d0
296
363
     $                       ,ptl+emass(i))
297
364
                        exit
298
365
                     elseif (is_a_lp(j) .and. mll.ne.0d0) then
299
366
                        if (j_fks.gt.nincoming)
300
 
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-ptl-emass(i) +
301
 
     &                                      max(mll/2d0,ptl+emass(i))
302
 
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-ptl-emass(i)
303
 
     $                       + max(mll/2d0,ptl+emass(i))
304
 
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-ptl-emass(i)
305
 
     $                       + max(mll/2d0,ptl+emass(i))
 
367
     &                       taumin(iFKS,ichan) = taumin(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2) +
 
368
     &                                      max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
 
369
                        taumin_s(iFKS,ichan) = taumin_s(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
370
     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
 
371
                        taumin_j(iFKS,ichan) = taumin_j(iFKS,ichan)-dsqrt(ptl**2+emass(i)**2)
 
372
     $                       + max(mll/2d0,dsqrt(ptl**2+emass(i)**2))
306
373
                        xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
307
374
     $                       +emass(i))
308
375
                        exit
309
376
                     endif
310
377
                  enddo
311
378
               else
312
 
                  taumin(iFKS,ichan)=taumin(iFKS,ichan)+emass(i)
313
 
                  taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan)+emass(i)
314
 
                  taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan)+emass(i)
315
 
                  xm(i)=emass(i)
 
379
                  if (i.eq.nexternal)then
 
380
                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + emass(i)
 
381
                        taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan) +  emass(i)
 
382
                        taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan) + emass(i)
 
383
                        xm(i) = emass(i)
 
384
                  else
 
385
                     smin_update = 0
 
386
                     nb_iden_pdg = 1
 
387
                     mxx = 0d0
 
388
c                    assume smin apply always on the same set of particle
 
389
                     do j=nincoming+1,nexternal-1
 
390
                        if (mxxmin(i,j).ne.0d0.or.mxxmin(j,i).ne.0d0) then
 
391
                           nb_iden_pdg = nb_iden_pdg +1
 
392
                           if (mxx.eq.0d0) mxx = max(mxxmin(i,j), mxxmin(j,i))
 
393
                        endif
 
394
                     enddo
 
395
                     ! S >= (2*N-N^2)*M1^2 + (N^2-N)/2 * Mxx^2
 
396
                     smin_update = nb_iden_pdg*((2-nb_iden_pdg)*emass(i)**2 + (nb_iden_pdg-1)/2.*mxx**2)
 
397
                     ! compare with the update from pt cut
 
398
                     if (smin_update.lt.nb_iden_pdg**2*(etmin(i)**2 + emass(i)**2))then
 
399
                        ! the pt is more restrictive
 
400
                        smin_update = dsqrt(etmin(i)**2 + emass(i)**2)
 
401
                     else
 
402
                        smin_update = dsqrt(smin_update)/nb_iden_pdg ! share over N particle, and change dimension
 
403
                     endif
 
404
                     smin_update = emass(i)
 
405
                     ! update in sqrt(s) so take the 
 
406
                     if  (j_fks.gt.nincoming) then
 
407
                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + smin_update
 
408
                     else
 
409
                        taumin(iFKS,ichan)=taumin(iFKS,ichan) + emass(i)
 
410
                     endif
 
411
                     taumin_s(iFKS,ichan)=taumin_s(iFKS,ichan) + smin_update
 
412
                     taumin_j(iFKS,ichan)=taumin_j(iFKS,ichan) + smin_update
 
413
                     xm(i) = smin_update
 
414
                  endif
316
415
               endif
317
416
               xw(i)=0d0
318
417
            enddo