1
=== modified file 'Template/NLO/Cards/run_card.dat'
2
--- Template/NLO/Cards/run_card.dat	2017-12-08 09:49:16 +0000
3
+++ Template/NLO/Cards/run_card.dat	2018-09-06 19:46:41 +0000
4
@@ -110,6 +110,22 @@
5
110
#***********************************************************************
110
#***********************************************************************
6
111
 %(store_rwgt_info)s = store_rwgt_info ! Store info for reweighting in LHE file
111
 %(store_rwgt_info)s = store_rwgt_info ! Store info for reweighting in LHE file
7
112
#***********************************************************************
112
#***********************************************************************
8
113
# iOSsubtr parameter: used if the process is generated with            *
9
114
# remove_os = True                                                     *
10
115
#  iossubtr = 1 -> DR without interferece                              *
11
116
#  iossubtr = 2 -> DR with interferece                                 *
12
117
#  iossubtr = 3 -> DS with reshuffling on initial state, standard BW   *
13
118
#  iossubtr = 4 -> DS with reshuffling on initial state, running BW    *
14
119
#  iossubtr = 5 -> DS with reshuffling on all FS particles, standard BW*
15
120
#  iossubtr = 6 -> DS with reshuffling on all FS particles, running BW *
16
121
#  iossubtr = 7 -> DS with option C, standard BW                       *
17
122
#  iossubtr = 8 -> DS with option C, running BW                        *
18
123
#  iossubtr = 9 -> DS with reshuffling on spectator, standard BW       *
19
124
#  iossubtr = 10-> DS with reshuffling on spectator, running BW        *
20
125
#***********************************************************************
21
126
 %(iossubtr)s = iossubtr ! strategy to be used to remove resonances 
22
127
                         ! appearing in real emissions
23
128
#***********************************************************************
24
113
# ickkw parameter:                                                     *
129
# ickkw parameter:                                                     *
25
114
#   0: No merging                                                      *
130
#   0: No merging                                                      *
26
115
#   3: FxFx Merging - WARNING! Applies merging only at the hard-event  *
131
#   3: FxFx Merging - WARNING! Applies merging only at the hard-event  *
27
116
132
28
=== modified file 'Template/NLO/Source/run.inc'
29
--- Template/NLO/Source/run.inc	2017-10-05 13:10:51 +0000
30
+++ Template/NLO/Source/run.inc	2018-09-06 19:46:41 +0000
31
@@ -74,10 +74,15 @@
32
74
      common/to_rwgt/ do_rwgt_scale, rw_Fscale_down, rw_Fscale_up, rw_Rscale_down, rw_Rscale_up,
74
      common/to_rwgt/ do_rwgt_scale, rw_Fscale_down, rw_Fscale_up, rw_Rscale_down, rw_Rscale_up,
33
75
     # do_rwgt_pdf, pdf_set_min, pdf_set_max,
75
     # do_rwgt_pdf, pdf_set_min, pdf_set_max,
34
76
     # store_rwgt_info 
76
     # store_rwgt_info 
35
77
36
78
C for the OS subtraction
37
79
      integer iossubtr
38
80
      common /to_os_reshuf/iossubtr
39
77
c
81
c
40
78
c    For FO run (with lhe type of analysis
82
c    For FO run (with lhe type of analysis
41
79
c
83
c
42
80
      double precision FO_LHE_weight_ratio
84
      double precision FO_LHE_weight_ratio
43
85
<<<<<<< TREE
44
81
      common /FO_ANALYSIS_LHW/FO_LHE_weight_ratio
86
      common /FO_ANALYSIS_LHW/FO_LHE_weight_ratio
45
82
c
87
c
46
83
c
88
c
47
@@ -87,4 +92,6 @@
48
87
      double precision ptmax4pdg(0:25) 
92
      double precision ptmax4pdg(0:25) 
49
88
      double precision mxxmin4pdg(0:25) 
93
      double precision mxxmin4pdg(0:25) 
50
89
      logical mxxpart_antipart(1:25) 
94
      logical mxxpart_antipart(1:25) 
51
90
      common/TO_PDG_SPECIFIC_CUT/pdg_cut, ptmin4pdg,ptmax4pdg, mxxmin4pdg, mxxpart_antipart
52
91
\ No newline at end of file
95
\ No newline at end of file
53
96
      common/TO_PDG_SPECIFIC_CUT/pdg_cut, ptmin4pdg,ptmax4pdg, mxxmin4pdg, mxxpart_antipart=======
54
97
      common /FO_ANALYSIS_LHW/FO_LHE_weight_ratio
55
98
>>>>>>> MERGE-SOURCE
56
92
99
57
=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
58
--- Template/NLO/SubProcesses/driver_mintFO.f	2018-06-28 15:09:53 +0000
59
+++ Template/NLO/SubProcesses/driver_mintFO.f	2018-09-06 19:46:41 +0000
60
@@ -90,10 +90,17 @@
61
90
      logical              fixed_order,nlo_ps
90
      logical              fixed_order,nlo_ps
62
91
      common /c_fnlo_nlops/fixed_order,nlo_ps
91
      common /c_fnlo_nlops/fixed_order,nlo_ps
63
92
92
64
93
C stuff for the on-shell subtraction
65
94
      integer os_countall, os_countbad
66
95
      common /to_os_count/ os_countall, os_countbad
67
96
68
93
97
69
94
C-----
98
C-----
70
95
C  BEGIN CODE
99
C  BEGIN CODE
71
96
C-----  
100
C-----  
72
101
      os_countall = 0
73
102
      os_countbad = 0
74
103
75
97
c
104
c
76
98
c     Setup the timing variable
105
c     Setup the timing variable
77
99
c
106
c
78
@@ -301,6 +308,8 @@
79
301
      write(*,*) 'Time spent in Write_events : ',t_write
308
      write(*,*) 'Time spent in Write_events : ',t_write
80
302
      write(*,*) 'Time spent in Other_tasks : ',tOther
309
      write(*,*) 'Time spent in Other_tasks : ',tOther
81
303
      write(*,*) 'Time spent in Total : ',tTot
310
      write(*,*) 'Time spent in Total : ',tTot
82
311
      write(*,*) 'Fraction of points in OS subtraction dead zone : ',
83
312
     $ dble(os_countbad)/max(dble(os_countall),1d0)
84
304
313
85
305
      open (unit=12, file='res.dat',status='unknown')
314
      open (unit=12, file='res.dat',status='unknown')
86
306
      do kchan=0,nchans
315
      do kchan=0,nchans
87
@@ -380,6 +389,8 @@
88
380
      common /for_applgrid/ iappl
389
      common /for_applgrid/ iappl
89
381
      double precision       wgt_ME_born,wgt_ME_real
390
      double precision       wgt_ME_born,wgt_ME_real
90
382
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
391
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
91
392
      ! stuff for the on-shell reshuffling
92
393
      common /to_reshuffle/x
93
383
      if (ifl.ne.0) then
394
      if (ifl.ne.0) then
94
384
         write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
395
         write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
95
385
         stop 1
396
         stop 1
96
386
397
97
=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
98
--- Template/NLO/SubProcesses/driver_mintMC.f	2018-06-28 15:09:53 +0000
99
+++ Template/NLO/SubProcesses/driver_mintMC.f	2018-09-06 19:46:41 +0000
100
@@ -95,9 +95,16 @@
101
95
      logical              fixed_order,nlo_ps
95
      logical              fixed_order,nlo_ps
102
96
      common /c_fnlo_nlops/fixed_order,nlo_ps
96
      common /c_fnlo_nlops/fixed_order,nlo_ps
103
97
97
104
98
C stuff for the on-shell subtraction
105
99
      integer os_countall, os_countbad
106
100
      common /to_os_count/ os_countall, os_countbad
107
101
108
98
C-----
102
C-----
109
99
C  BEGIN CODE
103
C  BEGIN CODE
110
100
C-----  
104
C-----  
111
105
      os_countall = 0
112
106
      os_countbad = 0
113
107
114
101
      call cpu_time(tBefore)
108
      call cpu_time(tBefore)
115
102
      fixed_order=.false.
109
      fixed_order=.false.
116
103
      nlo_ps=.true.
110
      nlo_ps=.true.
117
@@ -481,6 +488,8 @@
118
481
      write(*,*) 'Time spent in Write_events : ',t_write
488
      write(*,*) 'Time spent in Write_events : ',t_write
119
482
      write(*,*) 'Time spent in Other_tasks : ',tOther
489
      write(*,*) 'Time spent in Other_tasks : ',tOther
120
483
      write(*,*) 'Time spent in Total : ',tTot
490
      write(*,*) 'Time spent in Total : ',tTot
121
491
      write(*,*) 'Fraction of points in OS subtraction dead zone : ',
122
492
     $ dble(os_countbad)/max(dble(os_countall),1d0)
123
484
493
124
485
      open (unit=12, file='res.dat',status='unknown')
494
      open (unit=12, file='res.dat',status='unknown')
125
486
      if (imode.eq.0) then
495
      if (imode.eq.0) then
126
@@ -789,6 +798,8 @@
127
789
      common /c_imode/imode,only_virt
798
      common /c_imode/imode,only_virt
128
790
      double precision       wgt_ME_born,wgt_ME_real
799
      double precision       wgt_ME_born,wgt_ME_real
129
791
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
800
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
130
801
      ! stuff for the on-shell reshuffling
131
802
      common /to_reshuffle/x
132
792
      sigintF=0d0
803
      sigintF=0d0
133
793
c Find the nFKSprocess for which we compute the Born-like contributions
804
c Find the nFKSprocess for which we compute the Born-like contributions
134
794
      if (firsttime) then
805
      if (firsttime) then
135
795
806
136
=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
137
--- Template/NLO/SubProcesses/fks_singular.f	2018-03-12 07:20:30 +0000
138
+++ Template/NLO/SubProcesses/fks_singular.f	2018-09-06 19:46:41 +0000
139
@@ -3011,22 +3011,23 @@
140
3011
         wgt=wgt*xi_i_fks**2*(1d0-y_ij_fks)
3011
         wgt=wgt*xi_i_fks**2*(1d0-y_ij_fks)
141
3012
      endif
3012
      endif
142
3013
3013
159
3014
      if(wgt.lt.0.d0)then
3014
CMZ : real ME can be <0 with OS removal
160
3015
         icount=icount+1
3015
CC      if(wgt.lt.0.d0)then
161
3016
         if (icount.le.10) then
3016
CC         icount=icount+1
162
3017
            write(*,*) 'Warning, numerical problem found in sreal. '/
3017
CC         if (icount.le.10) then
163
3018
     $           /'Setting weight to zero',wgt,xi_i_fks,y_ij_fks
3018
CC            write(*,*) 'Warning, numerical problem found in sreal. '/
164
3019
            do i=1,nexternal
3019
CC     $           /'Setting weight to zero',wgt,xi_i_fks,y_ij_fks
165
3020
               write(*,*) 'particle ',i,', ',(pp(j,i),j=0,3)
3020
CC            do i=1,nexternal
166
3021
            enddo
3021
CC               write(*,*) 'particle ',i,', ',(pp(j,i),j=0,3)
167
3022
            if (icount.eq.25) then
3022
CC            enddo
168
3023
               write (*,*) 'ERROR 25 problems found... '/
3023
CC            if (icount.eq.25) then
169
3024
     $              /'stopping the code'
3024
CC               write (*,*) 'ERROR 25 problems found... '/
170
3025
               stop
3025
CC     $              /'stopping the code'
171
3026
            endif
3026
CC               stop
172
3027
         endif
3027
CC            endif
173
3028
         wgt=0d0
3028
CC         endif
174
3029
      endif
3029
CC         wgt=0d0
175
3030
CC      endif
176
3030
3031
177
3031
      return
3032
      return
178
3032
      end
3033
      end
179
3033
3034
180
=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
181
--- Template/NLO/SubProcesses/genps_fks.f	2017-06-16 08:20:04 +0000
182
+++ Template/NLO/SubProcesses/genps_fks.f	2018-09-06 19:46:41 +0000
183
@@ -166,6 +166,10 @@
184
166
     &     cBW_width(-1:1,-nexternal:-1)
166
     &     cBW_width(-1:1,-nexternal:-1)
185
167
      common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
167
      common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
186
168
     $     ,cBW_level
168
     $     ,cBW_level
187
169
C force Shat, for some reshuffling techniques with OS subtraction
188
170
      logical force_shat
189
171
      double precision shat_forced
190
172
      common/to_shat_genps/force_shat, shat_forced
191
169
173
192
170
      pass=.true.
174
      pass=.true.
193
171
      do i=1,nexternal-1
175
      do i=1,nexternal-1
194
@@ -236,7 +240,7 @@
195
236
c Generate Bjorken x's if need be and update jacobian c
240
c Generate Bjorken x's if need be and update jacobian c
196
237
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
241
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
197
238
      if (abs(lpp(1)).ge.1 .and. abs(lpp(2)).ge.1 .and.
242
      if (abs(lpp(1)).ge.1 .and. abs(lpp(2)).ge.1 .and.
199
239
     &     .not.(softtest.or.colltest)) then
243
     &     (.not.(softtest.or.colltest).or.force_shat)) then
200
240
c x(ndim-1) -> tau_cnt(0); x(ndim) -> ycm_cnt(0)
244
c x(ndim-1) -> tau_cnt(0); x(ndim) -> ycm_cnt(0)
201
241
         if (one_body) then
245
         if (one_body) then
202
242
c tau is fixed by the mass of the final state particle
246
c tau is fixed by the mass of the final state particle
203
@@ -254,6 +258,10 @@
204
254
               call generate_tau(stot,ndim-4,x(ndim-4),tau_born,xjac0)
258
               call generate_tau(stot,ndim-4,x(ndim-4),tau_born,xjac0)
205
255
            endif
259
            endif
206
256
         endif
260
         endif
207
261
C for On-shell subtraction, if shat is forced from outside, set tau here
208
262
         if (force_shat) then
209
263
           tau_born=shat_forced/stot
210
264
         endif
211
257
c Generate the rapditity of the Born system
265
c Generate the rapditity of the Born system
212
258
         call generate_y(tau_born,x(ndim-3),ycm_born,ycmhat,xjac0)
266
         call generate_y(tau_born,x(ndim-3),ycm_born,ycmhat,xjac0)
213
259
      elseif (abs(lpp(1)).ge.1 .and.
267
      elseif (abs(lpp(1)).ge.1 .and.
214
260
268
215
=== modified file 'Template/NLO/SubProcesses/makefile_fks_dir'
216
--- Template/NLO/SubProcesses/makefile_fks_dir	2017-07-29 07:52:02 +0000
217
+++ Template/NLO/SubProcesses/makefile_fks_dir	2018-09-06 19:46:41 +0000
218
@@ -25,17 +25,29 @@
219
25
$(LIBDIR)libpdf.a $(LIBDIR)libcernlib.a
25
$(LIBDIR)libpdf.a $(LIBDIR)libcernlib.a
220
26
26
221
27
# Files for all executables
27
# Files for all executables
222
28
<<<<<<< TREE
223
28
FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) $(patsubst	\
29
FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) $(patsubst	\
224
29
	 %.f,%.o,$(wildcard matrix_*.f)) real_me_chooser.o		\
30
	 %.f,%.o,$(wildcard matrix_*.f)) real_me_chooser.o		\
225
30
	 leshouche_inc_chooser.o fks_inc_chooser.o			\
31
	 leshouche_inc_chooser.o fks_inc_chooser.o			\
226
32
=======
227
33
FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) 			\
228
34
	 $(patsubst %.f,%.o,$(wildcard wrapper_matrix_*.f)) 		\
229
35
	 $(patsubst %.f,%.o,$(wildcard matrix_*.f)) real_me_chooser.o	\
230
36
	 fks_inc_chooser.o leshouche_inc_chooser.o genps_fks.o		\
231
37
>>>>>>> MERGE-SOURCE
232
31
	 setcuts.o setscales.o born.o sborn_sf.o $(patsubst		\
38
	 setcuts.o setscales.o born.o sborn_sf.o $(patsubst		\
233
32
	 %.f,%.o,$(wildcard b_sf_???.f)) fks_Sij.o			\
39
	 %.f,%.o,$(wildcard b_sf_???.f)) fks_Sij.o			\
234
33
	 $(fastjetfortran_madfks) fks_singular.o veto_xsec.o		\
40
	 $(fastjetfortran_madfks) fks_singular.o veto_xsec.o		\
235
34
	 montecarlocounter.o reweight_xsec.o boostwdir2.o       	\
41
	 montecarlocounter.o reweight_xsec.o boostwdir2.o       	\
236
35
	 configs_and_props_inc_chooser.o initcluster.o cluster.o 	\
42
	 configs_and_props_inc_chooser.o initcluster.o cluster.o 	\
237
36
	 reweight.o get_color.o FKSParamReader.o iproc_map.o		\
43
	 reweight.o get_color.o FKSParamReader.o iproc_map.o		\
238
44
<<<<<<< TREE
239
37
	 MC_integer.o $(reweight_xsec_events_pdf_dummy)			\
45
	 MC_integer.o $(reweight_xsec_events_pdf_dummy)			\
240
38
	 $(applgrid_interface) genps_fks.o weight_lines.o
46
	 $(applgrid_interface) genps_fks.o weight_lines.o
241
47
=======
242
48
	 MC_integer.o transform_os.o $(reweight_xsec_events_pdf_dummy)	\
243
49
	 $(applgrid_interface)
244
50
>>>>>>> MERGE-SOURCE
245
39
51
246
40
# Files needed for mintFO & mintMC
52
# Files needed for mintFO & mintMC
247
41
RUN= $(FO_ANALYSE) $(FILES) cuts.o pythia_unlops.o recluster.o 		\
53
RUN= $(FO_ANALYSE) $(FILES) cuts.o pythia_unlops.o recluster.o 		\
248
42
54
249
=== added file 'Template/NLO/SubProcesses/symmetry_fks_test_MC.f.OTHER'
250
--- Template/NLO/SubProcesses/symmetry_fks_test_MC.f.OTHER	1970-01-01 00:00:00 +0000
251
+++ Template/NLO/SubProcesses/symmetry_fks_test_MC.f.OTHER	2018-09-06 19:46:41 +0000
252
@@ -0,0 +1,675 @@
253
1
      program symmetry
254
2
c*****************************************************************************
255
3
c     Given identical particles, and the configurations. This program identifies
256
4
c     identical configurations and specifies which ones can be skipped
257
5
c*****************************************************************************
258
6
      implicit none
259
7
c
260
8
c     Constants
261
9
c
262
10
      include 'genps.inc'      
263
11
      include 'nexternal.inc'
264
12
      include '../../Source/run_config.inc'
265
13
      include 'nFKSconfigs.inc'
266
14
      include 'fks_info.inc'
267
15
      include 'run.inc'
268
16
      include 'cuts.inc'
269
17
      include 'mint.inc'
270
18
      
271
19
      double precision ZERO,one
272
20
      parameter       (ZERO = 0d0)
273
21
      parameter       (one = 1d0)
274
22
      integer   maxswitch
275
23
      parameter(maxswitch=99)
276
24
c
277
25
c     Local
278
26
c
279
27
      integer iforest(2,-max_branch:-1,lmaxconfigs)
280
28
      integer mapconfig(0:lmaxconfigs)
281
29
      integer sprop(-max_branch:-1,lmaxconfigs)
282
30
      integer itree(2,-max_branch:-1)
283
31
      integer imatch
284
32
      integer i,j, k, n, nsym,l,ii,jj
285
33
      double precision diff,xi_i_fks
286
34
c$$$      double precision pmass(-max_branch:-1,lmaxconfigs)   !Propagotor mass
287
35
      double precision pmass(nexternal)
288
36
      double precision pwidth(-max_branch:-1,lmaxconfigs)  !Propagotor width
289
37
      integer pow(-max_branch:-1,lmaxconfigs)
290
38
291
39
      integer biforest(2,-max_branch:-1,lmaxconfigs)
292
40
      integer fksmother,fksgrandmother,fksaunt,compare
293
41
      integer fksconfiguration,mapbconf(0:lmaxconfigs)
294
42
      integer r2b(lmaxconfigs),b2r(lmaxconfigs)
295
43
      logical searchforgranny,is_beta_cms,is_granny_sch,topdown,non_prop
296
44
      integer nbranch,ns_channel,nt_channel
297
45
c      include "fks.inc"
298
46
      integer fks_j_from_i(nexternal,0:nexternal)
299
47
     &     ,particle_type(nexternal),pdg_type(nexternal)
300
48
      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
301
49
      double precision fxl,limit(15),wlimit(15)
302
50
      double precision lxp(0:3,nexternal+1),xp(15,0:3,nexternal+1)
303
51
      double precision fks_Sij
304
52
      double precision check,tolerance,zh,h_damp
305
53
      double precision gfactsf
306
54
      parameter (tolerance=1.d-4)
307
55
      integer kk,ll,bs,bs_min,bs_max,iconfig_in
308
56
309
57
      integer nsofttests,ncolltests,nerr,imax,iflag,iret,ilim
310
58
311
59
c alsf and besf are the parameters that control gfunsoft
312
60
      double precision alsf,besf
313
61
      common/cgfunsfp/alsf,besf
314
62
c alazi and beazi are the parameters that control gfunazi
315
63
      double precision alazi,beazi
316
64
      common/cgfunazi/alazi,beazi
317
65
318
66
c
319
67
c     Local for generating amps
320
68
c
321
69
      double precision p(0:3,99), wgt, x(99), fx
322
70
      double complex wgt1(2)
323
71
      double precision p1(0:3,99)
324
72
      integer ninvar, ndim,  minconfig, maxconfig
325
73
      common/tosigint/ndim
326
74
      integer ncall,itmax,nconfigs,ntry
327
75
      integer ic(nexternal,maxswitch), jc(12),nswitch
328
76
      integer nmatch, ibase
329
77
      logical mtc, even
330
78
331
79
      double precision totmass
332
80
333
81
      double precision xi_i_fks_fix_save,y_ij_fks_fix_save
334
82
      double precision xi_i_fks_fix,y_ij_fks_fix
335
83
      common/cxiyfix/xi_i_fks_fix,y_ij_fks_fix
336
84
c
337
85
c     Global
338
86
c
339
87
      include 'coupl.inc'
340
88
341
89
      logical calculatedBorn
342
90
      common/ccalculatedBorn/calculatedBorn
343
91
344
92
      integer i_fks,j_fks
345
93
      common/fks_indices/i_fks,j_fks
346
94
347
95
      double precision p1_cnt(0:3,nexternal,-2:2)
348
96
      double precision wgt_cnt(-2:2)
349
97
      double precision pswgt_cnt(-2:2)
350
98
      double precision jac_cnt(-2:2)
351
99
      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
352
100
353
101
      double precision p_born(0:3,nexternal-1)
354
102
      common/pborn/p_born
355
103
356
104
      double precision xi_i_fks_ev,y_ij_fks_ev
357
105
      double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
358
106
      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
359
107
360
108
      double precision xi_i_fks_cnt(-2:2)
361
109
      common /cxiifkscnt/xi_i_fks_cnt
362
110
363
111
      logical rotategranny
364
112
      common/crotategranny/rotategranny
365
113
366
114
      logical softtest,colltest
367
115
      common/sctests/softtest,colltest
368
116
      
369
117
      logical xexternal
370
118
      common /toxexternal/ xexternal
371
119
372
120
c Particle types (=color) of i_fks, j_fks and fks_mother
373
121
      integer i_type,j_type,m_type
374
122
      common/cparticle_types/i_type,j_type,m_type
375
123
376
124
c
377
125
c     External
378
126
c
379
127
      logical check_swap
380
128
      double precision dsig,ran2
381
129
      external dsig,ran2
382
130
      external check_swap, fks_Sij
383
131
384
132
c define here the maximum fraction of failures to consider the test
385
133
c   passed
386
134
      double precision max_fail, fail_frac
387
135
      parameter (max_fail=0.3d0)
388
136
389
137
c helicity stuff
390
138
      integer          isum_hel
391
139
      logical                    multi_channel
392
140
      common/to_matrix/isum_hel, multi_channel
393
141
      logical Hevents
394
142
      common/SHevents/Hevents
395
143
396
144
      integer fks_conf_number,fks_loop_min,fks_loop_max,fks_loop
397
145
      INTEGER NFKSPROCESS
398
146
      COMMON/C_NFKSPROCESS/NFKSPROCESS
399
147
      LOGICAL  IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
400
148
      LOGICAL  IS_A_PH(NEXTERNAL)
401
149
      COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
402
150
      
403
151
      character*10 MonteCarlo
404
152
      common/cMonteCarloType/MonteCarlo
405
153
406
154
      double precision shower_S_scale(fks_configs*2)
407
155
     &     ,shower_H_scale(fks_configs*2),ref_H_scale(fks_configs*2)
408
156
     &     ,pt_hardness
409
157
      common /cshowerscale2/shower_S_scale,shower_H_scale,ref_H_scale
410
158
     &     ,pt_hardness
411
159
412
160
      logical new_point
413
161
      common /c_new_point/new_point
414
162
      ! stuff for the on-shell reshuffling
415
163
      common /to_reshuffle/x
416
164
417
165
c     integer icomp
418
166
c
419
167
c     DATA
420
168
c
421
169
      integer tprid(-max_branch:-1,lmaxconfigs)
422
170
      include 'born_conf.inc'
423
171
c-----
424
172
c  Begin Code
425
173
c-----
426
174
      if (fks_configs.eq.1) then
427
175
         if (pdg_type_d(1,fks_i_d(1)).eq.-21) then
428
176
            write (*,*) 'Process generated with [LOonly=QCD]. '/
429
177
     $           /'No tests to do.'
430
178
            return
431
179
         endif
432
180
      endif
433
181
      write(*,*)'Enter the Monte Carlo name: possible choices are'
434
182
      write(*,*)'HERWIG6, HERWIGPP, PYTHIA6Q, PYTHIA6PT, PYTHIA8'
435
183
      read(*,*)MonteCarlo
436
184
      if(MonteCarlo.ne.'HERWIG6'.and.MonteCarlo.ne.'HERWIGPP'.and.
437
185
     &   MonteCarlo.ne.'PYTHIA6Q'.and.MonteCarlo.ne.'PYTHIA6PT'.and.
438
186
     &   MonteCarlo.ne.'PYTHIA8')then
439
187
         write(*,*)'Wrong name ',MonteCarlo,' during the tests'
440
188
         stop
441
189
      endif
442
190
443
191
      write(*,*)'Enter 0 to compute MC/MC(limit)'
444
192
      write(*,*)'      1 to compute MC/ME(limit)'
445
193
      read(*,*)ilim
446
194
447
195
      write(6,*)'Enter alpha, beta for G_soft'
448
196
      write(6,*)'  Enter alpha<0 to set G_soft=1 (no ME soft)'
449
197
      read(5,*)alsf,besf
450
198
451
199
      write(6,*)'Enter alpha, beta for G_azi'
452
200
      write(6,*)'  Enter alpha>0 to set G_azi=0 (no azi corr)'
453
201
      read(5,*)alazi,beazi
454
202
455
203
      write(*,*)'Enter xi_i, y_ij to be used in coll/soft tests'
456
204
      write(*,*)' Enter -2 to generate them randomly'
457
205
      read(*,*)xi_i_fks_fix_save,y_ij_fks_fix_save
458
206
459
207
      write(*,*)'Enter number of tests for soft and collinear limits'
460
208
      read(*,*)nsofttests,ncolltests
461
209
462
210
      write(*,*)'Sum over helicity (0), or random helicity (1)'
463
211
      read(*,*) isum_hel
464
212
465
213
      write(*,*)'  '
466
214
      write(*,*)'  '
467
215
      write(*,*)'****************************************************'
468
216
      write(*,*)'****************************************************'
469
217
      write(*,*)'            Testing limits for ',MonteCarlo
470
218
      write(*,*)'****************************************************'
471
219
      write(*,*)'****************************************************'
472
220
      write(*,*)'  '
473
221
      write(*,*)'  '
474
222
475
223
476
224
      call setrun                !Sets up run parameters
477
225
      call setpara('param_card.dat')   !Sets up couplings and masses
478
226
      call setcuts               !Sets up cuts 
479
227
480
228
c When doing hadron-hadron collision reduce the effect collision energy.
481
229
c Note that tests are always performed at fixed energy with Bjorken x=1.
482
230
      totmass = 0.0d0
483
231
      include 'pmass.inc' ! make sure to set the masses after the model has been included
484
232
      do i=nincoming+1,nexternal
485
233
         if (is_a_j(i) .and. i.ne.nexternal) then
486
234
            totmass = totmass + max(ptj,pmass(i))
487
235
         elseif ((is_a_lp(i).or.is_a_lm(i)) .and. i.ne.nexternal) then
488
236
            totmass = totmass + max(mll/2d0,mll_sf/2d0,ptl,pmass(i))
489
237
         else
490
238
            totmass = totmass + pmass(i)
491
239
         endif
492
240
      enddo
493
241
      if (lpp(1).ne.0) ebeam(1)=max(ebeam(1)/20d0,totmass)
494
242
      if (lpp(2).ne.0) ebeam(2)=max(ebeam(2)/20d0,totmass)
495
243
c
496
244
497
245
      write (*,*) 'Give FKS configuration number ("0" loops over all)'
498
246
      read (*,*) fks_conf_number
499
247
500
248
      if (fks_conf_number.eq.0) then
501
249
         fks_loop_min=1
502
250
         fks_loop_max=fks_configs
503
251
      else
504
252
         fks_loop_min=fks_conf_number
505
253
         fks_loop_max=fks_conf_number
506
254
      endif
507
255
508
256
      do fks_loop=fks_loop_min,fks_loop_max
509
257
         nFKSprocess=fks_loop
510
258
         write (*,*) ''
511
259
         write (*,*) '================================================='
512
260
         write (*,*) ''
513
261
         write (*,*) 'NEW FKS CONFIGURATION:'
514
262
515
263
         call fks_inc_chooser()
516
264
         call leshouche_inc_chooser()
517
265
         call setcuts
518
266
         write (*,*) 'FKS configuration number is ',nFKSprocess
519
267
         write (*,*) 'FKS partons are: i=',i_fks,'  j=',j_fks
520
268
         write (*,*) 'with PDGs:       i=',PDG_type(i_fks),'  j='
521
269
     $        ,PDG_type(j_fks)
522
270
c
523
271
      ndim = 55
524
272
      ncall = 10000
525
273
      itmax = 10
526
274
      ninvar = 35
527
275
      nconfigs = 1
528
276
529
277
c Set color types of i_fks, j_fks and fks_mother.
530
278
      i_type=particle_type(i_fks)
531
279
      j_type=particle_type(j_fks)
532
280
      if (abs(i_type).eq.abs(j_type)) then
533
281
         m_type=8
534
282
         if ( (j_fks.le.nincoming .and.
535
283
     &        abs(i_type).eq.3 .and. j_type.ne.i_type) .or.
536
284
     &        (j_fks.gt.nincoming .and.
537
285
     &        abs(i_type).eq.3 .and. j_type.ne.-i_type)) then
538
286
            write(*,*)'Flavour mismatch #1 in setfksfactor',
539
287
     &           i_fks,j_fks,i_type,j_type
540
288
            stop
541
289
         endif
542
290
      elseif(abs(i_type).eq.3 .and. j_type.eq.8)then
543
291
         if(j_fks.le.nincoming)then
544
292
            m_type=-i_type
545
293
         else
546
294
            write (*,*) 'Error in setfksfactor: (i,j)=(q,g)'
547
295
            stop
548
296
         endif
549
297
      elseif(i_type.eq.8 .and. abs(j_type).eq.3)then
550
298
         if (j_fks.le.nincoming) then
551
299
            m_type=j_type
552
300
         else
553
301
            m_type=j_type
554
302
         endif
555
303
      else
556
304
         write(*,*)'Flavour mismatch #2 in setfksfactor',
557
305
     &        i_type,j_type,m_type
558
306
         stop
559
307
      endif
560
308
561
309
562
310
c     
563
311
c     Get momentum configuration
564
312
c
565
313
566
314
c Set xexternal to true to use the x's from external vegas in the
567
315
c x_to_f_arg subroutine
568
316
      xexternal=.true.
569
317
      
570
318
      write(*,*)'  '
571
319
      write(*,*)'  '
572
320
      write(*,*)'Enter graph number (iconfig), '
573
321
     &     //"'0' loops over all graphs"
574
322
      read(*,*)iconfig_in
575
323
      
576
324
      if (iconfig_in.eq.0) then
577
325
         bs_min=1
578
326
         bs_max=mapconfig(0)
579
327
      elseif (iconfig_in.eq.-1) then
580
328
         bs_min=1
581
329
         bs_max=1
582
330
      else
583
331
         bs_min=iconfig_in
584
332
         bs_max=iconfig_in
585
333
      endif
586
334
587
335
      do iconfig=bs_min,bs_max       ! Born configurations
588
336
         ichan=1
589
337
         iconfigs(1)=iconfig
590
338
591
339
      call setfksfactor(.true.)
592
340
c      call set_mc_matrices
593
341
594
342
      wgt=1d0
595
343
      ntry=1
596
344
597
345
      softtest=.false.
598
346
      colltest=.false.
599
347
600
348
      do jj=1,ndim
601
349
         x(jj)=ran2()
602
350
      enddo
603
351
      new_point=.true.
604
352
      call generate_momenta(ndim,iconfig,wgt,x,p)
605
353
      calculatedBorn=.false.
606
354
      do while (( wgt.lt.0 .or. p(0,1).le.0d0 .or. p_born(0,1).le.0d0
607
355
     &           ) .and. ntry .lt. 1000)
608
356
         do jj=1,ndim
609
357
            x(jj)=ran2()
610
358
         enddo
611
359
         new_point=.true.
612
360
         wgt=1d0
613
361
         call generate_momenta(ndim,iconfig,wgt,x,p)
614
362
         calculatedBorn=.false.
615
363
         ntry=ntry+1
616
364
      enddo
617
365
618
366
      if (ntry.ge.1000) then
619
367
         write (*,*) 'No points passed cuts...'
620
368
         write (12,*) 'ERROR: no points passed cuts...'
621
369
     &        //' Cannot perform ME tests properly for config',iconfig
622
370
         exit
623
371
      endif
624
372
625
373
626
374
      call sborn(p_born,wgt1)
627
375
628
376
      write (*,*) ''
629
377
      write (*,*) ''
630
378
      write (*,*) ''
631
379
632
380
      Hevents=.true.
633
381
      softtest=.true.
634
382
      colltest=.false.
635
383
      nerr=0
636
384
      imax=10
637
385
      do j=1,nsofttests
638
386
         call get_helicity(i_fks,j_fks)
639
387
         if(nsofttests.le.10)then
640
388
           write (*,*) ' '
641
389
           write (*,*) ' '
642
390
         endif
643
391
644
392
         y_ij_fks_fix=y_ij_fks_fix_save
645
393
         xi_i_fks_fix=0.1d0
646
394
         ntry=1
647
395
         wgt=1d0
648
396
         do jj=1,ndim
649
397
            x(jj)=ran2()
650
398
         enddo
651
399
         new_point=.true.
652
400
         call generate_momenta(ndim,iconfig,wgt,x,p)
653
401
         do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
654
402
            wgt=1d0
655
403
            do jj=1,ndim
656
404
               x(jj)=ran2()
657
405
            enddo
658
406
            new_point=.true.
659
407
            call generate_momenta(ndim,iconfig,wgt,x,p)
660
408
            ntry=ntry+1
661
409
         enddo
662
410
         if(nsofttests.le.10)write (*,*) 'ntry',ntry
663
411
c Set xi_i_fks to zero, to correctly generate the collinear momenta for the
664
412
c configurations close to the soft-collinear limit
665
413
         xi_i_fks_fix=0.d0
666
414
         wgt=1d0
667
415
         call generate_momenta(ndim,iconfig,wgt,x,p)
668
416
         calculatedBorn=.false.
669
417
         call set_cms_stuff(0)
670
418
         calculatedBorn=.false.
671
419
672
420
c Initialise shower_S_scale to a large value, not to get spurious dead zones
673
421
         shower_S_scale=1d20
674
422
675
423
         if(ilim.eq.0)then
676
424
           call xmcsubt_wrap(p1_cnt(0,1,0),zero,y_ij_fks_ev,fxl)
677
425
         else
678
426
           call sreal(p1_cnt(0,1,0),zero,y_ij_fks_ev,fxl) 
679
427
         endif
680
428
         fxl=fxl*jac_cnt(0)
681
429
682
430
c Now generate the momenta for the original xi_i_fks=0.1, slightly shifted,
683
431
c because otherwise fresh random will be used...
684
432
         xi_i_fks_fix=0.100001d0
685
433
         wgt=1d0
686
434
         call generate_momenta(ndim,iconfig,wgt,x,p)
687
435
         calculatedBorn=.false.
688
436
689
437
         call set_cms_stuff(-100)
690
438
         call xmcsubt_wrap(p,xi_i_fks_ev,y_ij_fks_ev,fx)
691
439
         limit(1)=fx*wgt
692
440
         wlimit(1)=wgt
693
441
694
442
         do k=1,nexternal
695
443
            do l=0,3
696
444
               lxp(l,k)=p1_cnt(l,k,0)
697
445
               xp(1,l,k)=p(l,k)
698
446
            enddo
699
447
         enddo
700
448
         do l=0,3
701
449
            lxp(l,nexternal+1)=p_i_fks_cnt(l,0)
702
450
            xp(1,l,nexternal+1)=p_i_fks_ev(l)
703
451
         enddo
704
452
705
453
         do i=2,imax
706
454
            xi_i_fks_fix=xi_i_fks_fix/10d0
707
455
            wgt=1d0
708
456
            call generate_momenta(ndim,iconfig,wgt,x,p)
709
457
            calculatedBorn=.false.
710
458
            call set_cms_stuff(-100)
711
459
            call xmcsubt_wrap(p,xi_i_fks_ev,y_ij_fks_ev,fx)
712
460
            limit(i)=fx*wgt
713
461
            wlimit(i)=wgt
714
462
            do k=1,nexternal
715
463
               do l=0,3
716
464
                  xp(i,l,k)=p(l,k)
717
465
               enddo
718
466
            enddo
719
467
            do l=0,3
720
468
               xp(i,l,nexternal+1)=p_i_fks_ev(l)
721
469
            enddo
722
470
         enddo
723
471
724
472
         if(nsofttests.le.10)then
725
473
           write (*,*) 'Soft limit:'
726
474
           do i=1,imax
727
475
              call xprintout(6,limit(i),fxl)
728
476
           enddo
729
477
c
730
478
           write(80,*)'  '
731
479
           write(80,*)'****************************'
732
480
           write(80,*)'  '
733
481
           do k=1,nexternal+1
734
482
              write(80,*)''
735
483
              write(80,*)'part:',k
736
484
              do l=0,3
737
485
                 write(80,*)'comp:',l
738
486
                 do i=1,10
739
487
                    call xprintout(80,xp(i,l,k),lxp(l,k))
740
488
                 enddo
741
489
              enddo
742
490
           enddo
743
491
        else
744
492
           iflag=0
745
493
           call checkres(limit,fxl,wlimit,jac_cnt(0),xp,lxp,
746
494
     &                   iflag,imax,j,nexternal,i_fks,j_fks,iret)
747
495
           nerr=nerr+iret
748
496
        endif
749
497
750
498
      enddo
751
499
      if(nsofttests.gt.10)then
752
500
         write(*,*)'Soft tests done for (Born) config',iconfig
753
501
         write(*,*)'Failures:',nerr
754
502
         fail_frac= nerr/dble(nsofttests)
755
503
         if (fail_frac.lt.max_fail) then
756
504
             write(*,401) nFKSprocess, fail_frac
757
505
         else
758
506
             write(*,402) nFKSprocess, fail_frac
759
507
         endif
760
508
      endif
761
509
762
510
      write (*,*) ''
763
511
      write (*,*) ''
764
512
      write (*,*) ''
765
513
      
766
514
      include 'pmass.inc'
767
515
768
516
      if (pmass(j_fks).ne.0d0) then
769
517
         write (*,*) 'No collinear test for massive j_fks'
770
518
         goto 123
771
519
      endif
772
520
773
521
      softtest=.false.
774
522
      colltest=.true.
775
523
776
524
c Set rotategranny=.true. to align grandmother along the z axis, when 
777
525
c grandmother is not the c.m. system (if granny=cms, this rotation coincides
778
526
c with the identity, and the following is harmless).
779
527
c WARNING: the setting of rotategranny changes the definition of xij_aor
780
528
c in genps_fks_test.f
781
529
      rotategranny=.false.
782
530
783
531
      nerr=0
784
532
      imax=10
785
533
      do j=1,ncolltests
786
534
         call get_helicity(i_fks,j_fks)
787
535
788
536
         if(ncolltests.le.10)then
789
537
            write (*,*) ' '
790
538
            write (*,*) ' '
791
539
         endif
792
540
793
541
         y_ij_fks_fix=0.9d0
794
542
         xi_i_fks_fix=xi_i_fks_fix_save
795
543
         ntry=1
796
544
         wgt=1d0
797
545
         do jj=1,ndim
798
546
            x(jj)=ran2()
799
547
         enddo
800
548
         new_point=.true.
801
549
         call generate_momenta(ndim,iconfig,wgt,x,p)
802
550
         do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
803
551
            wgt=1d0
804
552
            do jj=1,ndim
805
553
               x(jj)=ran2()
806
554
            enddo
807
555
            new_point=.true.
808
556
            call generate_momenta(ndim,iconfig,wgt,x,p)
809
557
            ntry=ntry+1
810
558
         enddo
811
559
         calculatedBorn=.false.
812
560
         if(ncolltests.le.10)write (*,*) 'ntry',ntry
813
561
         call set_cms_stuff(1)
814
562
         if(ilim.eq.0)then
815
563
           call xmcsubt_wrap(p1_cnt(0,1,1),xi_i_fks_cnt(1),one,fxl)
816
564
         else
817
565
           call sreal(p1_cnt(0,1,1),xi_i_fks_cnt(1),one,fxl) 
818
566
         endif
819
567
         fxl=fxl*jac_cnt(1)
820
568
         call set_cms_stuff(-100)
821
569
         call xmcsubt_wrap(p,xi_i_fks_ev,y_ij_fks_ev,fx)
822
570
         limit(1)=fx*wgt
823
571
         wlimit(1)=wgt
824
572
825
573
         do k=1,nexternal
826
574
            do l=0,3
827
575
               lxp(l,k)=p1_cnt(l,k,1)
828
576
               xp(1,l,k)=p(l,k)
829
577
            enddo
830
578
         enddo
831
579
         do l=0,3
832
580
            lxp(l,nexternal+1)=p_i_fks_cnt(l,1)
833
581
            xp(1,l,nexternal+1)=p_i_fks_ev(l)
834
582
         enddo
835
583
836
584
         do i=2,imax
837
585
            y_ij_fks_fix=1-0.1d0**i
838
586
            wgt=1d0
839
587
            call generate_momenta(ndim,iconfig,wgt,x,p)
840
588
            calculatedBorn=.false.
841
589
            call set_cms_stuff(-100)
842
590
            call xmcsubt_wrap(p,xi_i_fks_ev,y_ij_fks_ev,fx)
843
591
            limit(i)=fx*wgt
844
592
            wlimit(i)=wgt
845
593
            do k=1,nexternal
846
594
               do l=0,3
847
595
                  xp(i,l,k)=p(l,k)
848
596
               enddo
849
597
            enddo
850
598
            do l=0,3
851
599
               xp(i,l,nexternal+1)=p_i_fks_ev(l)
852
600
            enddo
853
601
         enddo
854
602
         if(ncolltests.le.10)then
855
603
            write (*,*) 'Collinear limit:'
856
604
            do i=1,imax
857
605
               call xprintout(6,limit(i),fxl)
858
606
            enddo
859
607
c
860
608
            write(80,*)'  '
861
609
            write(80,*)'****************************'
862
610
            write(80,*)'  '
863
611
            do k=1,nexternal+1
864
612
               write(80,*)''
865
613
               write(80,*)'part:',k
866
614
               do l=0,3
867
615
                  write(80,*)'comp:',l
868
616
                  do i=1,10
869
617
                     call xprintout(80,xp(i,l,k),lxp(l,k))
870
618
                  enddo
871
619
               enddo
872
620
            enddo
873
621
         else
874
622
            iflag=1
875
623
            call checkres(limit,fxl,wlimit,jac_cnt(1),xp,lxp,
876
624
     &                    iflag,imax,j,nexternal,i_fks,j_fks,iret)
877
625
            nerr=nerr+iret
878
626
         endif
879
627
      enddo
880
628
      if(ncolltests.gt.10)then
881
629
         write(*,*)'Collinear tests done for (Born) config', iconfig
882
630
         write(*,*)'Failures:',nerr
883
631
         fail_frac= nerr/dble(ncolltests)
884
632
         if (fail_frac.lt.max_fail) then
885
633
             write(*,501) nfksprocess, fail_frac
886
634
         else
887
635
             write(*,502) nfksprocess, fail_frac
888
636
         endif
889
637
      endif
890
638
 123  continue
891
639
      enddo                     ! Loop over Born configurations
892
640
      enddo                     ! Loop over nFKSprocess
893
641
894
642
      return
895
643
 401  format('     Soft test ',i2,' PASSED. Fraction of failures: ',
896
644
     & f4.2) 
897
645
 402  format('     Soft test ',I2,' FAILED. Fraction of failures: ',
898
646
     & f4.2) 
899
647
 501  format('Collinear test ',i2,' PASSED. Fraction of failures: ',
900
648
     & f4.2) 
901
649
 502  format('Collinear test ',I2,' FAILED. Fraction of failures: ',
902
650
     & f4.2) 
903
651
      end
904
652
905
653
c
906
654
c
907
655
c Dummy routines
908
656
c
909
657
c
910
658
      subroutine clear_events()
911
659
      end
912
660
      subroutine initplot
913
661
      end
914
662
      subroutine store_events()
915
663
      end
916
664
      integer function n_unwgted()
917
665
      n_unwgted = 1
918
666
      end
919
667
920
668
      subroutine outfun(pp,www)
921
669
      implicit none
922
670
      include 'nexternal.inc'
923
671
      real*8 pp(0:3,nexternal),www
924
672
c
925
673
      write(*,*)'This routine should not be called here'
926
674
      stop
927
675
      end
928
0
676
929
=== added file 'Template/NLO/SubProcesses/symmetry_fks_test_ME.f.OTHER'
930
--- Template/NLO/SubProcesses/symmetry_fks_test_ME.f.OTHER	1970-01-01 00:00:00 +0000
931
+++ Template/NLO/SubProcesses/symmetry_fks_test_ME.f.OTHER	2018-09-06 19:46:41 +0000
932
@@ -0,0 +1,598 @@
933
1
      program symmetry
934
2
c*****************************************************************************
935
3
c     Given identical particles, and the configurations. This program identifies
936
4
c     identical configurations and specifies which ones can be skipped
937
5
c*****************************************************************************
938
6
      implicit none
939
7
c
940
8
c     Constants
941
9
c
942
10
      include 'genps.inc'      
943
11
      include 'nexternal.inc'
944
12
      include '../../Source/run_config.inc'
945
13
      include 'nFKSconfigs.inc'
946
14
      include 'fks_info.inc'
947
15
      include 'run.inc'
948
16
      include 'cuts.inc'
949
17
      include 'mint.inc'
950
18
      
951
19
      double precision ZERO,one
952
20
      parameter       (ZERO = 0d0)
953
21
      parameter       (one = 1d0)
954
22
      integer   maxswitch
955
23
      parameter(maxswitch=99)
956
24
c
957
25
c     Local
958
26
c
959
27
      integer iforest(2,-max_branch:-1,lmaxconfigs)
960
28
      integer mapconfig(0:lmaxconfigs)
961
29
      integer sprop(-max_branch:-1,lmaxconfigs)
962
30
      integer itree(2,-max_branch:-1)
963
31
      integer imatch
964
32
      integer i,j, k, n, nsym,l,ii,jj
965
33
      double precision diff,xi_i_fks
966
34
c$$$      double precision pmass(-max_branch:-1,lmaxconfigs)   !Propagotor mass
967
35
      double precision pmass(nexternal)
968
36
      double precision pwidth(-max_branch:-1,lmaxconfigs)  !Propagotor width
969
37
      integer pow(-max_branch:-1,lmaxconfigs)
970
38
971
39
      integer biforest(2,-max_branch:-1,lmaxconfigs)
972
40
      integer fksmother,fksgrandmother,fksaunt,compare
973
41
      integer fksconfiguration,mapbconf(0:lmaxconfigs)
974
42
      integer r2b(lmaxconfigs),b2r(lmaxconfigs)
975
43
      logical searchforgranny,is_beta_cms,is_granny_sch,topdown,non_prop
976
44
      integer nbranch,ns_channel,nt_channel
977
45
c      include "fks.inc"
978
46
      integer fks_j_from_i(nexternal,0:nexternal)
979
47
     &     ,particle_type(nexternal),pdg_type(nexternal)
980
48
      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
981
49
      double precision fxl,limit(15),wlimit(15)
982
50
      double precision lxp(0:3,nexternal+1),xp(15,0:3,nexternal+1)
983
51
      double precision fks_Sij
984
52
      double precision check,tolerance,zh,h_damp
985
53
      parameter (tolerance=1.d-4)
986
54
      integer kk,ll,bs,bs_min,bs_max,iconfig_in
987
55
988
56
      integer nsofttests,ncolltests,nerr,imax,iflag,iret
989
57
c
990
58
c     Local for generating amps
991
59
c
992
60
      double precision p(0:3,99), wgt, x(99), fx
993
61
      double complex wgt1(2)
994
62
      double precision p1(0:3,99)
995
63
      integer ninvar, ndim,  minconfig, maxconfig
996
64
      common/tosigint/ndim
997
65
      integer ncall,itmax,nconfigs,ntry
998
66
      integer ic(nexternal,maxswitch), jc(12),nswitch
999
67
      integer nmatch, ibase
1000
68
      logical mtc, even
1001
69
1002
70
      double precision totmass
1003
71
1004
72
      double precision xi_i_fks_fix_save,y_ij_fks_fix_save
1005
73
      double precision xi_i_fks_fix,y_ij_fks_fix
1006
74
      common/cxiyfix/xi_i_fks_fix,y_ij_fks_fix
1007
75
c
1008
76
c     Global
1009
77
c
1010
78
      include 'coupl.inc'
1011
79
1012
80
      logical calculatedBorn
1013
81
      common/ccalculatedBorn/calculatedBorn
1014
82
1015
83
      integer i_fks,j_fks
1016
84
      common/fks_indices/i_fks,j_fks
1017
85
1018
86
      double precision p1_cnt(0:3,nexternal,-2:2)
1019
87
      double precision wgt_cnt(-2:2)
1020
88
      double precision pswgt_cnt(-2:2)
1021
89
      double precision jac_cnt(-2:2)
1022
90
      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
1023
91
1024
92
      double precision p_born(0:3,nexternal-1)
1025
93
      common/pborn/p_born
1026
94
1027
95
      double precision xi_i_fks_ev,y_ij_fks_ev
1028
96
      double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
1029
97
      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
1030
98
1031
99
      double precision xi_i_fks_cnt(-2:2)
1032
100
      common /cxiifkscnt/xi_i_fks_cnt
1033
101
1034
102
      logical rotategranny
1035
103
      common/crotategranny/rotategranny
1036
104
1037
105
      logical softtest,colltest
1038
106
      common/sctests/softtest,colltest
1039
107
      
1040
108
      logical xexternal
1041
109
      common /toxexternal/ xexternal
1042
110
1043
111
c Particle types (=color) of i_fks, j_fks and fks_mother
1044
112
      integer i_type,j_type,m_type
1045
113
      common/cparticle_types/i_type,j_type,m_type
1046
114
1047
115
c
1048
116
c     External
1049
117
c
1050
118
      logical check_swap
1051
119
      double precision dsig,ran2
1052
120
      external dsig,ran2
1053
121
      external check_swap, fks_Sij
1054
122
1055
123
c define here the maximum fraction of failures to consider the test
1056
124
c   passed
1057
125
      double precision max_fail, fail_frac
1058
126
      parameter (max_fail=0.3d0)
1059
127
1060
128
c helicity stuff
1061
129
      integer          isum_hel
1062
130
      logical                    multi_channel
1063
131
      common/to_matrix/isum_hel, multi_channel
1064
132
1065
133
      integer fks_conf_number,fks_loop_min,fks_loop_max,fks_loop
1066
134
      INTEGER NFKSPROCESS
1067
135
      COMMON/C_NFKSPROCESS/NFKSPROCESS
1068
136
      LOGICAL  IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
1069
137
      LOGICAL  IS_A_PH(NEXTERNAL)
1070
138
      COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
1071
139
      
1072
140
      logical new_point
1073
141
      common /c_new_point/new_point
1074
142
      ! stuff for the on-shell reshuffling
1075
143
      common /to_reshuffle/x
1076
144
c      integer icomp
1077
145
c
1078
146
c     DATA
1079
147
c
1080
148
      integer tprid(-max_branch:-1,lmaxconfigs)
1081
149
      include 'born_conf.inc'
1082
150
c-----
1083
151
c  Begin Code
1084
152
c-----
1085
153
      if (fks_configs.eq.1) then
1086
154
         if (pdg_type_d(1,fks_i_d(1)).eq.-21) then
1087
155
            write (*,*) 'Process generated with [LOonly=QCD]. '/
1088
156
     $           /'No tests to do.'
1089
157
            return
1090
158
         endif
1091
159
      endif
1092
160
1093
161
      write(*,*)'Enter xi_i, y_ij to be used in coll/soft tests'
1094
162
      write(*,*)' Enter -2 to generate them randomly'
1095
163
      read(*,*)xi_i_fks_fix_save,y_ij_fks_fix_save
1096
164
1097
165
      write(*,*)'Enter number of tests for soft and collinear limits'
1098
166
      read(*,*)nsofttests,ncolltests
1099
167
1100
168
      write(*,*)'Sum over helicity (0), or random helicity (1)'
1101
169
      read(*,*) isum_hel
1102
170
1103
171
      call setrun                !Sets up run parameters
1104
172
      call setpara('param_card.dat')   !Sets up couplings and masses
1105
173
      call setcuts               !Sets up cuts 
1106
174
1107
175
c When doing hadron-hadron collision reduce the effect collision energy.
1108
176
c Note that tests are always performed at fixed energy with Bjorken x=1.
1109
177
      totmass = 0.0d0
1110
178
      include 'pmass.inc' ! make sure to set the masses after the model has been included
1111
179
      do i=nincoming+1,nexternal
1112
180
         if (is_a_j(i) .and. i.ne.nexternal) then
1113
181
            totmass = totmass + max(ptj,pmass(i))
1114
182
         elseif ((is_a_lp(i).or.is_a_lm(i)) .and. i.ne.nexternal) then
1115
183
            totmass = totmass + max(mll/2d0,mll_sf/2d0,ptl,pmass(i))
1116
184
         else
1117
185
            totmass = totmass + pmass(i)
1118
186
         endif
1119
187
      enddo
1120
188
      if (lpp(1).ne.0) ebeam(1)=max(ebeam(1)/20d0,totmass)
1121
189
      if (lpp(2).ne.0) ebeam(2)=max(ebeam(2)/20d0,totmass)
1122
190
c
1123
191
1124
192
      write (*,*) 'Give FKS configuration number ("0" loops over all)'
1125
193
      read (*,*) fks_conf_number
1126
194
1127
195
      if (fks_conf_number.eq.0) then
1128
196
         fks_loop_min=1
1129
197
         fks_loop_max=fks_configs
1130
198
      else
1131
199
         fks_loop_min=fks_conf_number
1132
200
         fks_loop_max=fks_conf_number
1133
201
      endif
1134
202
1135
203
      do fks_loop=fks_loop_min,fks_loop_max
1136
204
         nFKSprocess=fks_loop
1137
205
         write (*,*) ''
1138
206
         write (*,*) '================================================='
1139
207
         write (*,*) ''
1140
208
         write (*,*) 'NEW FKS CONFIGURATION:'
1141
209
1142
210
         call fks_inc_chooser()
1143
211
         call leshouche_inc_chooser()
1144
212
         write (*,*) 'FKS configuration number is ',nFKSprocess
1145
213
         write (*,*) 'FKS partons are: i=',i_fks,'  j=',j_fks
1146
214
         write (*,*) 'with PDGs:       i=',PDG_type(i_fks),'  j='
1147
215
     $        ,PDG_type(j_fks)
1148
216
1149
217
c
1150
218
      ndim = 55
1151
219
      ncall = 10000
1152
220
      itmax = 10
1153
221
      ninvar = 35
1154
222
      nconfigs = 1
1155
223
1156
224
c Set color types of i_fks, j_fks and fks_mother.
1157
225
      i_type=particle_type(i_fks)
1158
226
      j_type=particle_type(j_fks)
1159
227
      if (abs(i_type).eq.abs(j_type)) then
1160
228
         m_type=8
1161
229
         if ( (j_fks.le.nincoming .and.
1162
230
     &        abs(i_type).eq.3 .and. j_type.ne.i_type) .or.
1163
231
     &        (j_fks.gt.nincoming .and.
1164
232
     &        abs(i_type).eq.3 .and. j_type.ne.-i_type)) then
1165
233
            write(*,*)'Flavour mismatch #1 in setfksfactor',
1166
234
     &           i_fks,j_fks,i_type,j_type
1167
235
            stop
1168
236
         endif
1169
237
      elseif(abs(i_type).eq.3 .and. j_type.eq.8)then
1170
238
         if(j_fks.le.nincoming)then
1171
239
            m_type=-i_type
1172
240
         else
1173
241
            write (*,*) 'Error in setfksfactor: (i,j)=(q,g)'
1174
242
            stop
1175
243
         endif
1176
244
      elseif(i_type.eq.8 .and. abs(j_type).eq.3)then
1177
245
         if (j_fks.le.nincoming) then
1178
246
            m_type=j_type
1179
247
         else
1180
248
            m_type=j_type
1181
249
         endif
1182
250
      else
1183
251
         write(*,*)'Flavour mismatch #2 in setfksfactor',
1184
252
     &        i_type,j_type,m_type
1185
253
         stop
1186
254
      endif
1187
255
1188
256
1189
257
c     
1190
258
c     Get momentum configuration
1191
259
c
1192
260
1193
261
c Set xexternal to true to use the x's from external vegas in the
1194
262
c x_to_f_arg subroutine
1195
263
      xexternal=.true.
1196
264
      
1197
265
      write(*,*)'  '
1198
266
      write(*,*)'  '
1199
267
      write(*,*)'Enter graph number (iconfig), '
1200
268
     &     //"'0' loops over all graphs"
1201
269
      read(*,*)iconfig_in
1202
270
      
1203
271
      if (iconfig_in.eq.0) then
1204
272
         bs_min=1
1205
273
         bs_max=mapconfig(0)
1206
274
      elseif (iconfig_in.eq.-1) then
1207
275
         bs_min=1
1208
276
         bs_max=1
1209
277
      else
1210
278
         bs_min=iconfig_in
1211
279
         bs_max=iconfig_in
1212
280
      endif
1213
281
1214
282
      do iconfig=bs_min,bs_max  ! Born configurations
1215
283
         ichan=1
1216
284
         iconfigs(1)=iconfig
1217
285
      call setcuts
1218
286
      call setfksfactor(.false.)
1219
287
      wgt=1d0
1220
288
      ntry=1
1221
289
1222
290
      softtest=.false.
1223
291
      colltest=.false.
1224
292
1225
293
      do jj=1,ndim
1226
294
         x(jj)=ran2()
1227
295
      enddo
1228
296
      new_point=.true.
1229
297
      call generate_momenta(ndim,iconfig,wgt,x,p)
1230
298
      calculatedBorn=.false.
1231
299
      do while (( wgt.lt.0 .or. p(0,1).le.0d0 .or. p_born(0,1).le.0d0
1232
300
     &           ) .and. ntry .lt. 1000)
1233
301
         do jj=1,ndim
1234
302
            x(jj)=ran2()
1235
303
         enddo
1236
304
         new_point=.true.
1237
305
         wgt=1d0
1238
306
         call generate_momenta(ndim,iconfig,wgt,x,p)
1239
307
         calculatedBorn=.false.
1240
308
         ntry=ntry+1
1241
309
      enddo
1242
310
1243
311
      if (ntry.ge.1000) then
1244
312
         write (*,*) 'No points passed cuts...'
1245
313
         write (12,*) 'ERROR: no points passed cuts...'
1246
314
     &        //' Cannot perform ME tests properly for config',iconfig
1247
315
         exit
1248
316
      endif
1249
317
1250
318
      call sborn(p_born,wgt1)
1251
319
      
1252
320
      write (*,*) ''
1253
321
      write (*,*) ''
1254
322
      write (*,*) ''
1255
323
1256
324
      softtest=.true.
1257
325
      colltest=.false.
1258
326
      nerr=0
1259
327
      imax=10
1260
328
      do j=1,nsofttests
1261
329
      call get_helicity(i_fks,j_fks)
1262
330
1263
331
         if(nsofttests.le.10)then
1264
332
           write (*,*) ' '
1265
333
           write (*,*) ' '
1266
334
         endif
1267
335
1268
336
         y_ij_fks_fix=y_ij_fks_fix_save
1269
337
         xi_i_fks_fix=0.1d0
1270
338
         ntry=1
1271
339
         wgt=1d0
1272
340
         do jj=1,ndim
1273
341
            x(jj)=ran2()
1274
342
         enddo
1275
343
         new_point=.true.
1276
344
         call generate_momenta(ndim,iconfig,wgt,x,p)
1277
345
         do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
1278
346
            wgt=1d0
1279
347
            do jj=1,ndim
1280
348
               x(jj)=ran2()
1281
349
            enddo
1282
350
            new_point=.true.
1283
351
            call generate_momenta(ndim,iconfig,wgt,x,p)
1284
352
            ntry=ntry+1
1285
353
         enddo
1286
354
         if(nsofttests.le.10)write (*,*) 'ntry',ntry
1287
355
         calculatedBorn=.false.
1288
356
         call set_cms_stuff(0)
1289
357
         call sreal(p1_cnt(0,1,0),zero,y_ij_fks_ev,fxl) 
1290
358
         fxl=fxl*jac_cnt(0)
1291
359
1292
360
         call set_cms_stuff(-100)
1293
361
         call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx)
1294
362
         limit(1)=fx*wgt
1295
363
         wlimit(1)=wgt
1296
364
1297
365
         do k=1,nexternal
1298
366
            do l=0,3
1299
367
               lxp(l,k)=p1_cnt(l,k,0)
1300
368
               xp(1,l,k)=p(l,k)
1301
369
            enddo
1302
370
         enddo
1303
371
         do l=0,3
1304
372
            lxp(l,nexternal+1)=p_i_fks_cnt(l,0)
1305
373
            xp(1,l,nexternal+1)=p_i_fks_ev(l)
1306
374
         enddo
1307
375
1308
376
         do i=2,imax
1309
377
            xi_i_fks_fix=xi_i_fks_fix/10d0
1310
378
            wgt=1d0
1311
379
            call generate_momenta(ndim,iconfig,wgt,x,p)
1312
380
            calculatedBorn=.false.
1313
381
            call set_cms_stuff(-100)
1314
382
            call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx)
1315
383
            limit(i)=fx*wgt
1316
384
            wlimit(i)=wgt
1317
385
            do k=1,nexternal
1318
386
               do l=0,3
1319
387
                  xp(i,l,k)=p(l,k)
1320
388
               enddo
1321
389
            enddo
1322
390
            do l=0,3
1323
391
               xp(i,l,nexternal+1)=p_i_fks_ev(l)
1324
392
            enddo
1325
393
         enddo
1326
394
1327
395
         if(nsofttests.le.10)then
1328
396
           write (*,*) 'Soft limit:'
1329
397
           do i=1,imax
1330
398
              call xprintout(6,limit(i),fxl)
1331
399
           enddo
1332
400
c
1333
401
           write(80,*)'  '
1334
402
           write(80,*)'****************************'
1335
403
           write(80,*)'  '
1336
404
           do k=1,nexternal+1
1337
405
              write(80,*)''
1338
406
              write(80,*)'part:',k
1339
407
              do l=0,3
1340
408
                 write(80,*)'comp:',l
1341
409
                 do i=1,10
1342
410
                    call xprintout(80,xp(i,l,k),lxp(l,k))
1343
411
                 enddo
1344
412
              enddo
1345
413
           enddo
1346
414
        else
1347
415
           iflag=0
1348
416
           call checkres(limit,fxl,wlimit,jac_cnt(0),xp,lxp,
1349
417
     &                   iflag,imax,j,nexternal,i_fks,j_fks,iret)
1350
418
           nerr=nerr+iret
1351
419
        endif
1352
420
1353
421
      enddo
1354
422
      if(nsofttests.gt.10)then
1355
423
         write(*,*)'Soft tests done for (Born) config',iconfig
1356
424
         write(*,*)'Failures:',nerr
1357
425
         fail_frac= nerr/dble(nsofttests)
1358
426
         if (fail_frac.lt.max_fail) then
1359
427
             write(*,401) nFKSprocess, fail_frac
1360
428
         else
1361
429
             write(*,402) nFKSprocess, fail_frac
1362
430
         endif
1363
431
      endif
1364
432
1365
433
      write (*,*) ''
1366
434
      write (*,*) ''
1367
435
      write (*,*) ''
1368
436
      
1369
437
      include 'pmass.inc'
1370
438
1371
439
      if (pmass(j_fks).ne.0d0) then
1372
440
         write (*,*) 'No collinear test for massive j_fks'
1373
441
         goto 123
1374
442
      endif
1375
443
1376
444
      softtest=.false.
1377
445
      colltest=.true.
1378
446
1379
447
c Set rotategranny=.true. to align grandmother along the z axis, when 
1380
448
c grandmother is not the c.m. system (if granny=cms, this rotation coincides
1381
449
c with the identity, and the following is harmless).
1382
450
c WARNING: the setting of rotategranny changes the definition of xij_aor
1383
451
c in genps_fks_test.f
1384
452
      rotategranny=.false.
1385
453
1386
454
      nerr=0
1387
455
      imax=10
1388
456
      do j=1,ncolltests
1389
457
         call get_helicity(i_fks,j_fks)
1390
458
1391
459
         if(ncolltests.le.10)then
1392
460
            write (*,*) ' '
1393
461
            write (*,*) ' '
1394
462
         endif
1395
463
1396
464
         y_ij_fks_fix=0.9d0
1397
465
         xi_i_fks_fix=xi_i_fks_fix_save
1398
466
         ntry=1
1399
467
         wgt=1d0
1400
468
         do jj=1,ndim
1401
469
            x(jj)=ran2()
1402
470
         enddo
1403
471
         new_point=.true.
1404
472
         call generate_momenta(ndim,iconfig,wgt,x,p)
1405
473
         do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
1406
474
            wgt=1d0
1407
475
            do jj=1,ndim
1408
476
               x(jj)=ran2()
1409
477
            enddo
1410
478
            new_point=.true.
1411
479
            call generate_momenta(ndim,iconfig,wgt,x,p)
1412
480
            ntry=ntry+1
1413
481
         enddo
1414
482
         if(ncolltests.le.10)write (*,*) 'ntry',ntry
1415
483
         calculatedBorn=.false.
1416
484
         call set_cms_stuff(1)
1417
485
         call sreal(p1_cnt(0,1,1),xi_i_fks_cnt(1),one,fxl) 
1418
486
         fxl=fxl*jac_cnt(1)
1419
487
1420
488
         call set_cms_stuff(-100)
1421
489
         call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx)
1422
490
         limit(1)=fx*wgt
1423
491
         wlimit(1)=wgt
1424
492
1425
493
         do k=1,nexternal
1426
494
            do l=0,3
1427
495
               lxp(l,k)=p1_cnt(l,k,1)
1428
496
               xp(1,l,k)=p(l,k)
1429
497
            enddo
1430
498
         enddo
1431
499
         do l=0,3
1432
500
            lxp(l,nexternal+1)=p_i_fks_cnt(l,1)
1433
501
            xp(1,l,nexternal+1)=p_i_fks_ev(l)
1434
502
         enddo
1435
503
1436
504
         do i=2,imax
1437
505
            y_ij_fks_fix=1-0.1d0**i
1438
506
            wgt=1d0
1439
507
            call generate_momenta(ndim,iconfig,wgt,x,p)
1440
508
            calculatedBorn=.false.
1441
509
            call set_cms_stuff(-100)
1442
510
            call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx)
1443
511
            limit(i)=fx*wgt
1444
512
            wlimit(i)=wgt
1445
513
            do k=1,nexternal
1446
514
               do l=0,3
1447
515
                  xp(i,l,k)=p(l,k)
1448
516
               enddo
1449
517
            enddo
1450
518
            do l=0,3
1451
519
               xp(i,l,nexternal+1)=p_i_fks_ev(l)
1452
520
            enddo
1453
521
         enddo
1454
522
         if(ncolltests.le.10)then
1455
523
            write (*,*) 'Collinear limit:'
1456
524
            do i=1,imax
1457
525
               call xprintout(6,limit(i),fxl)
1458
526
            enddo
1459
527
c
1460
528
            write(80,*)'  '
1461
529
            write(80,*)'****************************'
1462
530
            write(80,*)'  '
1463
531
            do k=1,nexternal+1
1464
532
               write(80,*)''
1465
533
               write(80,*)'part:',k
1466
534
               do l=0,3
1467
535
                  write(80,*)'comp:',l
1468
536
                  do i=1,10
1469
537
                     call xprintout(80,xp(i,l,k),lxp(l,k))
1470
538
                  enddo
1471
539
               enddo
1472
540
            enddo
1473
541
         else
1474
542
            iflag=1
1475
543
            call checkres(limit,fxl,wlimit,jac_cnt(1),xp,lxp,
1476
544
     &                    iflag,imax,j,nexternal,i_fks,j_fks,iret)
1477
545
            nerr=nerr+iret
1478
546
         endif
1479
547
      enddo
1480
548
      if(ncolltests.gt.10)then
1481
549
         write(*,*)'Collinear tests done for (Born) config', iconfig
1482
550
         write(*,*)'Failures:',nerr
1483
551
         fail_frac= nerr/dble(ncolltests)
1484
552
         if (fail_frac.lt.max_fail) then
1485
553
             write(*,501) nFKSprocess, fail_frac
1486
554
         else
1487
555
             write(*,502) nFKSprocess, fail_frac
1488
556
         endif
1489
557
      endif
1490
558
1491
559
 123  continue
1492
560
1493
561
      enddo                     ! Loop over Born configurations
1494
562
      enddo                     ! Loop over nFKSprocess
1495
563
1496
564
1497
565
      return
1498
566
 401  format('     Soft test ',i2,' PASSED. Fraction of failures: ',
1499
567
     & f4.2) 
1500
568
 402  format('     Soft test ',I2,' FAILED. Fraction of failures: ',
1501
569
     & f4.2) 
1502
570
 501  format('Collinear test ',i2,' PASSED. Fraction of failures: ',
1503
571
     & f4.2) 
1504
572
 502  format('Collinear test ',I2,' FAILED. Fraction of failures: ',
1505
573
     & f4.2) 
1506
574
      end
1507
575
1508
576
c
1509
577
c
1510
578
c Dummy routines
1511
579
c
1512
580
c
1513
581
      subroutine clear_events()
1514
582
      end
1515
583
      subroutine initplot
1516
584
      end
1517
585
      subroutine store_events()
1518
586
      end
1519
587
      integer function n_unwgted()
1520
588
      n_unwgted = 1
1521
589
      end
1522
590
1523
591
      subroutine outfun(pp,www)
1524
592
      implicit none
1525
593
      include 'nexternal.inc'
1526
594
      real*8 pp(0:3,nexternal),www
1527
595
c
1528
596
      write(*,*)'This routine should not be called here'
1529
597
      stop
1530
598
      end
1531
0
599
1532
=== added file 'Template/NLO/SubProcesses/transform_os.f'
1533
--- Template/NLO/SubProcesses/transform_os.f	1970-01-01 00:00:00 +0000
1534
+++ Template/NLO/SubProcesses/transform_os.f	2018-09-06 19:46:41 +0000
1535
@@ -0,0 +1,1281 @@
1536
1
      subroutine transform_os_spect(p,q,ip,jp,kp,mass_i,mass_j,mass_k,mass_ij,stat)
1537
2
************************************************************************
1538
3
*     Authors: Dorival Goncalves & Marco Zaro                           *
1539
4
*     Given momenta p(nu,nexternal) produce q(nu,external) with the    *
1540
5
*     momentum qij on shell. The reshuffling is done as in the         *                             
1541
6
*     final-final CS dipole subtraction. kp is the spectator, ip and   *
1542
7
*     jp are the the decays products from the possible OS resonance    *
1543
8
************************************************************************
1544
9
      implicit none 
1545
10
      include 'nexternal.inc'
1546
11
      include 'coupl.inc'
1547
12
C-----Arguments      
1548
13
      double precision p(0:3,nexternal),q(0:3,nexternal)
1549
14
      double precision mass_i,mass_j,mass_k,mass_ij
1550
15
      integer ip,jp,kp
1551
16
      integer stat
1552
17
C-----Local
1553
18
      integer i,j,nu
1554
19
      double precision qtot(0:3),qsq,pij(0:3),qij(0:3),qz(0:3),dot,threedot
1555
20
      double precision lambda_tr,msq_ij,msq_k,msq_i,msq_j,a,b,c,d
1556
21
      external dot,threedot,lambda_tr
1557
22
      double precision zero, one, two
1558
23
      parameter (zero=0d0,one=1d0,two=2d0)
1559
24
 
1560
25
C---------------                                                                                   
1561
26
C     BEGIN CODE                                                                    
1562
27
C---------------    
1563
28
      stat=0 ! everything is fine for the moment
1564
29
1565
30
      msq_i=mass_i*mass_i
1566
31
      msq_k=mass_k*mass_k
1567
32
      msq_j=mass_j*mass_j
1568
33
      msq_ij=mass_ij*mass_ij
1569
34
1570
35
      do nu=0,3
1571
36
         qtot(nu)=p(nu,ip)+p(nu,jp)+p(nu,kp)
1572
37
         pij(nu)=p(nu,ip)+p(nu,jp)
1573
38
      enddo
1574
39
      qsq=dot(qtot,qtot)
1575
40
1576
41
C----- check that the reshuffling can be done; otherwise return stat=1
1577
42
      if (dsqrt(qsq).lt.mass_ij+mass_k) then
1578
43
        stat=1
1579
44
        return
1580
45
      endif
1581
46
      
1582
47
      do nu=0,3
1583
48
C-----spectator kp            
1584
49
         q(nu,kp)=dsqrt(lambda_tr(qsq,msq_ij,msq_k))/
1585
50
     &        dsqrt(lambda_tr(qsq,dot(pij,pij),msq_k))
1586
51
     &        *(p(nu,kp)-dot(qtot,p(0,kp))/qsq*qtot(nu))
1587
52
     &        +(qsq+msq_k-msq_ij)/two/qsq*qtot(nu)
1588
53
         
1589
54
C-----OS particle
1590
55
         qij(nu)=qtot(nu)-q(nu,kp)
1591
56
      enddo
1592
57
      
1593
58
C-  i'(or j')  momemntum: chose the one which has mass=0
1594
59
C 2 i' dot qij is fixed
1595
60
C assume i' direction is the same as i one
1596
61
      if (mass_i.eq.0) then
1597
62
        a = (msq_ij -msq_j) / 2d0 / 
1598
63
     $   (qij(0) * dsqrt(threedot(p(0,ip),p(0,ip))) -
1599
64
     $    threedot(p(0,ip),qij))
1600
65
        do nu = 1, 3
1601
66
          q(nu,ip) = a * p(nu,ip)
1602
67
        enddo
1603
68
        q(0,ip) = dsqrt(threedot(q(0,ip), q(0,ip)))
1604
69
        do nu = 0, 3
1605
70
          q(nu,jp) = qij(nu) - q(nu,ip)
1606
71
        enddo
1607
72
      else if (mass_j.eq.0) then
1608
73
        a = (msq_ij -msq_i) / 2d0 / 
1609
74
     $   (qij(0) * dsqrt(threedot(p(0,jp),p(0,jp))) -
1610
75
     $    threedot(p(0,jp),qij))
1611
76
        do nu = 1, 3
1612
77
          q(nu,jp) = a * p(nu,jp)
1613
78
        enddo
1614
79
        q(0,jp) = dsqrt(threedot(q(0,jp), q(0,jp)))
1615
80
        do nu = 0, 3
1616
81
          q(nu,ip) = qij(nu) - q(nu,jp)
1617
82
        enddo
1618
83
      else 
1619
84
        write(*,*) 'ERROR IN TRANSFORM_OS_SPECT: I and J both massive', mass_i, mass_j
1620
85
        stop
1621
86
      endif
1622
87
1623
88
      do nu=0,3
1624
89
C-----Keeping the rest of the momenta without any changes
1625
90
         do j=1,nexternal
1626
91
            if(j.ne.kp .and. j.ne.ip .and. j.ne.jp)then
1627
92
               q(nu,j)=p(nu,j)
1628
93
            endif
1629
94
         enddo
1630
95
      enddo
1631
96
      
1632
97
C-----Consistency checks
1633
98
C--------mass shell condition for spectator
1634
99
      if (dabs(dot(q(0,kp),q(0,kp))-dot(p(0,kp),p(0,kp)))
1635
100
     $ .gt. 1d-4 * max(dot(p(0,kp),p(0,kp)), 1d0)) then
1636
101
        write(*,*) 'ERROR IN TRANSFORM_OS_SPECT: K NOT KEPT ON SHELL'
1637
102
        write(*,*) 'MSQ before', dot(p(0,kp),p(0,kp))
1638
103
        write(*,*) 'MSQ after ', dot(q(0,kp),q(0,kp))
1639
104
        stop
1640
105
      endif
1641
106
1642
107
      ! check the momenta before returning
1643
108
      call OS_check_momenta(p, q, ip, jp)
1644
109
1645
110
      return     
1646
111
      end
1647
112
          
1648
113
1649
114
      REAL*8 function lambda_tr(x,y,z)
1650
115
C-----triangular function
1651
116
      implicit none
1652
117
      real*8 x,y,z
1653
118
      lambda_tr=x**2+y**2+z**2-2d0*x*y-2d0*x*z-2d0*y*z
1654
119
      return
1655
120
      end
1656
121
1657
122
1658
123
      double precision function lambda2(a,b,c)
1659
124
      implicit none
1660
125
      double precision a,b,c
1661
126
      if (a.le.0d0 .or. abs(b+c).gt.abs(a) .or. abs(b-c).gt.abs(a)) then
1662
127
         write (*,*) 'Error #1 in lambda2: inputs not consistent',a,b,c
1663
128
         stop 1
1664
129
      endif
1665
130
      lambda2=sqrt(1d0-(b+c)**2/a**2)*sqrt(1d0-(b-c)**2/a**2)
1666
131
      return
1667
132
      end
1668
133
1669
134
1670
135
      subroutine transform_os_genps(p,q,ip,jp,mass_i,mass_j,mass_ij,mom_perm,stat)
1671
136
************************************************************************
1672
137
*     Authors: Marco Zaro                                              *
1673
138
*     after rescaling the center-of-mass energy so that the reshuffling*
1674
139
*     can be done, it uses a reshuffling as described in sect 5.2 of   *
1675
140
*     0709.2092                                                        *
1676
141
*     Reshuffling routines have been coded by Rikkert Frederix         *
1677
142
************************************************************************
1678
143
      implicit none 
1679
144
      include 'nexternal.inc'
1680
145
C-----Arguments      
1681
146
      double precision p(0:3,nexternal), q(0:3,nexternal)
1682
147
      double precision mass_i, mass_j, mass_ij
1683
148
      integer ip, jp
1684
149
      integer mom_perm(nexternal)
1685
150
      integer stat 
1686
151
C-----Local
1687
152
      integer i, j, nu
1688
153
      double precision qtot(0:3), qsq, pij(0:3), qij(0:3), qi(0:3), qj(0:3)
1689
154
      double precision preco(0:3), qreco(0:3), ptmp(0:3)
1690
155
      double precision pcom(0:3,nexternal), qcom(0:3,nexternal), pboost(0:3)
1691
156
      double precision shat
1692
157
      double precision etot, ztot
1693
158
      double precision msq_ij, msq_k, msq_i, msq_j, msq_reco, resc
1694
159
      double precision mreco, memitt
1695
160
      double precision a, b
1696
161
      double precision dot, threedot, sumdot
1697
162
      external dot, threedot, sumdot
1698
163
      double precision rescale_init
1699
164
      double precision pmass(nexternal), sum_mass, ximass
1700
165
      !changing ximass may lead to catastrophic failures
1701
166
      parameter (ximass = 1d0) 
1702
167
      double precision p_shat(0:3, nexternal)
1703
168
1704
169
      include 'coupl.inc'
1705
170
      include 'run.inc'
1706
171
      double precision zero
1707
172
      parameter (zero=0d0)
1708
173
! variables to be passed to the subroutine that puts momenta on-shell
1709
174
      integer next
1710
175
      parameter (next=2*nexternal-3)
1711
176
      double precision p_os_in(0:4,next), p_os_out(0:4,next)
1712
177
      integer decay_tree(2,-next:-1,2)
1713
178
      integer i_os, skip
1714
179
1715
180
C--------------- 
1716
181
C     BEGIN CODE                                                                   
1717
182
C---------------    
1718
183
      include 'pmass.inc'
1719
184
1720
185
      stat = 0
1721
186
      msq_i=mass_i*mass_i
1722
187
      msq_j=mass_j*mass_j
1723
188
      msq_ij=mass_ij*mass_ij
1724
189
1725
190
c ----reconstruct the momentum pij
1726
191
      do nu=0,3
1727
192
        pij(nu)=p(nu,ip) + p(nu,jp)
1728
193
      enddo
1729
194
C the center of mass energy
1730
195
      shat = 2d0 * dot(p(0,1),p(0,2))
1731
196
C the sum of all the other masses (except for those in the reshuffling
1732
197
      sum_mass = 0d0
1733
198
      do i = nincoming+1, nexternal
1734
199
C        remember thatthe order of particles in the FS of the OS 
1735
200
C        counterterm may differ from the one in the rest of the code
1736
201
        if (mom_perm(i).ne.ip.and.mom_perm(i).ne.jp)
1737
202
     $   sum_mass = sum_mass + pmass(i)
1738
203
      enddo
1739
204
1740
205
C rescale shat (keep it below the collider energy) ; 
1741
206
      rescale_init=(sum_mass+mass_ij)/(sum_mass+min(ximass*mass_ij,dsqrt(dot(pij,pij))))
1742
207
      shat=min(shat*rescale_init**2, (ebeam(1)+ebeam(2))**2*0.99d0)
1743
208
1744
209
      if (dsqrt(shat).lt.sum_mass+mass_ij) then
1745
210
        write(*,*) 'ERROR IN TRANSFORM_OS_GENPS: INCONSITENT SHAT',
1746
211
     $   dsqrt(shat),sum_mass+mass_ij
1747
212
        stop 1
1748
213
      endif
1749
214
1750
215
      if (rescale_init.eq.1d0) then
1751
216
        ! if shat is the same, just copy the momenta
1752
217
        do i = 1, nexternal
1753
218
          do nu = 0, 3
1754
219
            p_shat(nu, i) = p(nu, i)
1755
220
          enddo
1756
221
        enddo
1757
222
      else
1758
223
        ! otherwise, call the PS generation routine to generate the momenta 
1759
224
        !  with the new shat
1760
225
        call generate_momenta_shat(shat,mom_perm,p_shat)
1761
226
      endif
1762
227
1763
228
      if (p_shat(0,1).le.0d0) then
1764
229
        write(*,*) 'RETURNING 1'
1765
230
        stat=1
1766
231
        return
1767
232
      endif
1768
233
1769
234
C boost the momenta in the partonic C.o.M. frame
1770
235
      do nu=0,3
1771
236
        pboost(nu) = p_shat(nu,1) + p_shat(nu,2)
1772
237
      enddo
1773
238
      do i=1,nexternal
1774
239
        call invboostx(p_shat(0,i), pboost, pcom(0,i))
1775
240
      enddo
1776
241
1777
242
      ! use the OS_resonance_and_decay_reshuffle function to put stuff
1778
243
      ! on-shell
1779
244
1780
245
      ! first, copy the argunment
1781
246
      ! input momentum : copy all particles but ip and jp, which will 
1782
247
      ! be put last
1783
248
      i_os = 0
1784
249
      do i = 1, nexternal
1785
250
        if (i.eq.ip.or.i.eq.jp) cycle
1786
251
        i_os = i_os + 1
1787
252
        do j = 0, 3 
1788
253
          p_os_in(j, i_os) = pcom(j, i)
1789
254
        enddo
1790
255
        ! it is perhaps a dirty solution, but should be ok
1791
256
        p_os_in(4,i_os) = dsqrt(max(dot(pcom(0,j),pcom(0,j)),0d0))
1792
257
      enddo
1793
258
1794
259
      ! safety check
1795
260
      if (i_os.ne.nexternal-2) then
1796
261
        write(*,*) 'ERROR IN TRANSFORM_OS_GENPS: INCONSITENT I_OS',
1797
262
     $     i_os, nexternal - 2
1798
263
        stop 1
1799
264
      endif
1800
265
1801
266
      ! then, add in order, the resonance (i+j), pi and pj
1802
267
      do nu = 0, 3
1803
268
        p_os_in(nu, nexternal - 1) = pcom(nu,ip) + pcom(nu,jp)
1804
269
        p_os_in(nu, nexternal) = pcom(nu,ip)
1805
270
        p_os_in(nu, nexternal + 1) = pcom(nu,jp)
1806
271
      enddo
1807
272
      p_os_in(4, nexternal - 1) = dsqrt(sumdot(pcom(0,ip), pcom(0,jp), 1d0))
1808
273
      p_os_in(4, nexternal) = mass_i
1809
274
      p_os_in(4, nexternal + 1) = mass_j
1810
275
1811
276
      ! finally, the decay tree
1812
277
      decay_tree(1, -1, 1) = nexternal
1813
278
      decay_tree(2, -1, 1) = nexternal + 1
1814
279
      decay_tree(1, -1, 2) = nexternal
1815
280
      decay_tree(2, -1, 2) = nexternal + 1
1816
281
1817
282
CC      subroutine OS_resonance_and_decay_reshuffle(npart,ibeta,decay_tree
1818
283
CC     $     ,Mbeta,p_os_in,p_os_out)
1819
284
      !  X  npart=nexternal (it is the number of particles before the
1820
285
      !           extra emission, so neternal-1 + the resonance
1821
286
      !  X  ibeta=nexternal -1 (the position of the resonance in p_os_in)
1822
287
      !  X  mbeta=mass_ij
1823
288
      call OS_resonance_and_decay_reshuffle(nexternal,nexternal-1,decay_tree
1824
289
     $     ,mass_ij,p_os_in,p_os_out)
1825
290
1826
291
      ! copy back from p_os_out to qcom
1827
292
      i_os = 0
1828
293
      do i = 1, nexternal
1829
294
        if (i.eq.ip.or.i.eq.jp) cycle
1830
295
        i_os = i_os + 1
1831
296
        do j = 0, 3 
1832
297
          qcom(j, i) = p_os_out(j, i_os)
1833
298
        enddo
1834
299
      enddo
1835
300
      ! as before, i is in p_os_out at position nexternal, j at nexternal+1 
1836
301
      do j = 0, 3 
1837
302
        qcom(j, ip) = p_os_out(j, nexternal)
1838
303
        qcom(j, jp) = p_os_out(j, nexternal + 1)
1839
304
      enddo
1840
305
1841
306
      ! finally boost back to the lab-frame
1842
307
      do i=1, nexternal
1843
308
        call boostx(qcom(0,i),pboost,q(0,i))
1844
309
      enddo
1845
310
1846
311
      ! check the momenta before returning
1847
312
      call OS_check_momenta(p, q, ip, jp)
1848
313
1849
314
      return     
1850
315
      end
1851
316
1852
317
1853
318
      subroutine transform_os_final(p,q,ip,jp,mass_i,mass_j,mass_ij,stat)
1854
319
************************************************************************
1855
320
*     Authors: Marco Zaro                                              *
1856
321
*     Given momenta p(nu,nexternal) produce q(nu,external) with the    *
1857
322
*     momentum qij on shell. The reshuffling is done by compensating   *
1858
323
*     on all the other final-state particles.                          *
1859
324
*     As in the case of transform_os_init, the angles in the           *
1860
325
*     pi+pj rest frame are not changed                                 *
1861
326
************************************************************************
1862
327
      implicit none 
1863
328
      include 'nexternal.inc'
1864
329
C-----Arguments      
1865
330
      double precision p(0:3,nexternal), q(0:3,nexternal)
1866
331
      double precision mass_i, mass_j, mass_ij
1867
332
      integer ip, jp
1868
333
      integer stat 
1869
334
C-----Local
1870
335
      integer i, j, nu
1871
336
      double precision qtot(0:3), qsq, pij(0:3), qij(0:3), qi(0:3), qj(0:3)
1872
337
      double precision preco(0:3), qreco(0:3), ptmp(0:3)
1873
338
      double precision pcom(0:3,nexternal), qcom(0:3,nexternal), pboost(0:3)
1874
339
      double precision shat
1875
340
      double precision etot, ztot
1876
341
      double precision msq_ij, msq_k, msq_i, msq_j, msq_reco, resc
1877
342
      double precision a, b
1878
343
      double precision dot, threedot
1879
344
      external dot, threedot
1880
345
      double precision rescale_init
1881
346
C--------------- 
1882
347
C     BEGIN CODE                                                                   
1883
348
C---------------    
1884
349
1885
350
      stat = 0
1886
351
      msq_i=mass_i*mass_i
1887
352
      msq_j=mass_j*mass_j
1888
353
      msq_ij=mass_ij*mass_ij
1889
354
1890
355
C the center of mass energy
1891
356
      shat = 2d0 * dot(p(0,1),p(0,2))
1892
357
1893
358
C boost the momenta in the partonic C.o.M. frame
1894
359
      do nu=0,3
1895
360
        pboost(nu) = p(nu,1) + p(nu,2)
1896
361
      enddo
1897
362
      do i=1,nexternal
1898
363
        call invboostx(p(0,i), pboost, pcom(0,i))
1899
364
      enddo
1900
365
1901
366
c reconstruct the recoil system (all FS particles which are not i and j)
1902
367
      do nu=0,3
1903
368
        preco(nu)=0d0
1904
369
      enddo
1905
370
      do i=nincoming+1,nexternal
1906
371
        if (i.eq.jp.or.i.eq.ip) cycle
1907
372
        do nu=0,3
1908
373
          preco(nu)=preco(nu)+pcom(nu,i)
1909
374
        enddo
1910
375
      enddo
1911
376
1912
377
      msq_reco=dot(preco,preco)
1913
378
1914
379
c ----reconstruct the momentum pij
1915
380
      do nu=0,3
1916
381
        pij(nu)=pcom(nu,ip) + pcom(nu,jp)
1917
382
      enddo
1918
383
1919
384
C check energy conditions
1920
385
      if (abs(preco(0)-sqrt(shat)/2d0*(1d0-(dot(pij,pij)-msq_reco)/shat))/shat.gt.1d-4) then
1921
386
          write(*,*) 'ERROR IN TRANSFORM_OS_FINAL: ENERGY #1'
1922
387
          write(*,*) 'shat', shat, sqrt(shat)
1923
388
          write(*,*) preco(0),sqrt(shat)/2d0*(1d0-(dot(pij,pij)-msq_reco)/shat)
1924
389
          stop
1925
390
      endif
1926
391
      if (abs(pij(0)-sqrt(shat)/2d0*(1d0+(dot(pij,pij)-msq_reco)/shat))/shat.gt.1d-4) then
1927
392
          write(*,*) 'ERROR IN TRANSFORM_OS_FINAL: ENERGY #2'
1928
393
          write(*,*) pij(0),sqrt(shat)/2d0*(1d0-(dot(pij,pij)-msq_reco)/shat)
1929
394
          stop
1930
395
      endif
1931
396
1932
397
C check if the reshuffling can be done (ie if there is enough energy).
1933
398
C otherwise, rescale shat and the initial momenta; 
1934
399
      if (dsqrt(shat).lt.mass_ij+dsqrt(msq_reco)) then
1935
400
        rescale_init=(dsqrt(msq_reco)+mass_ij)/(dsqrt(msq_reco)+dsqrt(dot(pij,pij)))
1936
401
        do nu = 0, 3
1937
402
          pcom(nu,1) = pcom(nu,1)*rescale_init
1938
403
          pcom(nu,2) = pcom(nu,2)*rescale_init
1939
404
        enddo
1940
405
        shat=shat*rescale_init**2
1941
406
      endif
1942
407
C alternatively one can put the momenta on the largest invariant mass that can be
1943
408
C generated (this is discarded for now)
1944
409
CC      if (dsqrt(shat).lt.mass_ij+dsqrt(msq_reco)) then
1945
410
CC          stat=1
1946
411
CC          msq_ij = max(0.99d0 * (dsqrt(shat)-dsqrt(msq_reco))**2, (mass_i+mass_j)**2)
1947
412
CC      endif
1948
413
1949
414
C the reshuffled momenta, qij and qreco, will have energy components
1950
415
C wchic correspond to qij having the reshuffled mass and qreco keeping
1951
416
C its invariant mass
1952
417
      qij(0) = dsqrt(shat)/2d0*(1d0+(msq_ij-msq_reco)/shat)
1953
418
      qreco(0) = dsqrt(shat)/2d0*(1d0-(msq_ij-msq_reco)/shat)
1954
419
C the other components have the same direction as pij, preco, and must
1955
420
C satisfy the mass-shell conditions qij^2 = m_ij^2, qreco^2=m_reco^2
1956
421
      do nu=1,3
1957
422
        qij(nu) = pij(nu)/dsqrt(threedot(pij,pij))*dsqrt(qij(0)*qij(0)-msq_ij)
1958
423
        qreco(nu) = preco(nu)/dsqrt(threedot(preco,preco))*dsqrt(qreco(0)*qreco(0)-msq_reco)
1959
424
      enddo
1960
425
1961
426
C now let us turn to the i and j particles and to the recoiling ones
1962
427
C *** IJ
1963
428
1964
429
c now go to the pij rest frame
1965
430
      call invboostx(pcom(0,ip), pij, qi)
1966
431
      call invboostx(pcom(0,jp), pij, qj)
1967
432
1968
433
C check that momenta are back to back
1969
434
      do nu=1,3
1970
435
        if (dabs(qi(nu)+qj(nu))/(qi(0)+qj(0)).gt.1d-4) then
1971
436
          write(*,*) 'ERROR IN TRANSFORM_OS_FINAL: QI QJ not B2B'
1972
437
          write(*,*) 'QI', qi, dsqrt(dot(qi,qi)), mass_i
1973
438
          write(*,*) 'QJ', qj, dsqrt(dot(qj,qj)), mass_j
1974
439
          stop
1975
440
        endif
1976
441
      enddo
1977
442
1978
443
C putting qi, qj on shell amounts in rescaling their spatial 
1979
444
c  component by a factor 
1980
445
      resc = dsqrt(
1981
446
     $ (((msq_ij - msq_i - msq_j)/2d0)**2 - msq_i * msq_j) / 
1982
447
     $ (threedot(qi,qi) * msq_ij))
1983
448
      do nu=1,3
1984
449
        qi(nu) = qi(nu) * resc
1985
450
        qj(nu) = qj(nu) * resc
1986
451
      enddo
1987
452
      qi(0) = dsqrt(msq_i + threedot(qi,qi))
1988
453
      qj(0) = dsqrt(msq_j + threedot(qj,qj))
1989
454
1990
455
c check that the sum of the energies is mij (use a low tolerance here)
1991
456
      if (dabs((dsqrt(msq_ij) - qi(0) - qj(0))/dsqrt(msq_ij)).gt.1d-2) then
1992
457
        write(*,*) 'ERROR IN TRANSFORM_OS_FINAL: QI QJ not on shell'
1993
458
        write(*,*) 'DIFF', dabs(dsqrt(msq_ij) - qi(0) - qj(0)), dabs(dsqrt(msq_ij) - qi(0) - qj(0))/dsqrt(msq_ij)
1994
459
        write(*,*) 'MIJ', dsqrt(msq_ij), stat
1995
460
        write(*,*) 'QI', qi, dsqrt(dot(qi,qi)), mass_i
1996
461
        write(*,*) 'QJ', qj, dsqrt(dot(qj,qj)), mass_j
1997
462
        stop
1998
463
      endif
1999
464
2000
465
C     boost qi,qj in the lab frame with qij
2001
466
      call boostx(qi, qij, qcom(0,ip))
2002
467
      call boostx(qj, qij, qcom(0,jp))
2003
468
2004
469
C check that momenta sum to qij
2005
470
      do nu=0,3
2006
471
        if (dabs(qcom(nu,ip)+qcom(nu,jp)-qij(nu))/(qij(0)).gt.1d-3) then
2007
472
          write(*,*) 'ERROR IN TRANSFORM_OS_FINAL: QIJ not conserved',nu
2008
473
          write(*,*) 'QI', qcom(0,ip), qcom(1,ip), qcom(2,ip), qcom(3,ip), dsqrt(dot(qcom(0,ip),qcom(0,ip)))
2009
474
          write(*,*) 'QJ', qcom(0,jp), qcom(1,jp), qcom(2,jp), qcom(3,jp), dsqrt(dot(qcom(0,jp),qcom(0,jp)))
2010
475
          write(*,*) 'QIJ', qij, dsqrt(dot(qij,qij)), dsqrt(msq_ij), stat
2011
476
          stop
2012
477
        endif
2013
478
      enddo
2014
479
2015
480
C *** other recoiling particles
2016
481
C boost them to the preco rest frame and then back to
2017
482
C the lab frame using qreco
2018
483
      do i=1, nexternal
2019
484
        if (i.eq.ip.or.i.eq.jp) cycle
2020
485
        call invboostx(pcom(0,i), preco, ptmp)
2021
486
        call boostx(ptmp, qreco, qcom(0,i))
2022
487
      enddo
2023
488
2024
489
C *** finally the initial state particles, do nothing
2025
490
      do i=1, nincoming
2026
491
        do nu=0,3
2027
492
          qcom(nu,i)=pcom(nu,i)
2028
493
        enddo
2029
494
      enddo
2030
495
2031
496
C *** boost the momenta back to the original frame
2032
497
      do i=1, nexternal
2033
498
        call boostx(qcom(0,i),pboost,q(0,i))
2034
499
      enddo
2035
500
2036
501
      ! check the momenta before returning
2037
502
      call OS_check_momenta(p, q, ip, jp)
2038
503
2039
504
      return     
2040
505
      end
2041
506
2042
507
2043
508
      subroutine transform_os_init(p,q,ip,jp,mass_i,mass_j,mass_ij)
2044
509
************************************************************************
2045
510
*     Authors: Marco Zaro                                              *
2046
511
*     Given momenta p(nu,nexternal) produce q(nu,external) with the    *
2047
512
*     momentum qij on shell. The reshuffling is done keeping the       * 
2048
513
*     3 momentum of pi+pj fixed, and changing the energy (the initial  *
2049
514
*     state momenta are also changed. Furthermore, the angles in the   *
2050
515
*     pi+pj rest frame are not changed                                 *
2051
516
************************************************************************
2052
517
      implicit none 
2053
518
      include 'nexternal.inc'
2054
519
C-----Arguments      
2055
520
      double precision p(0:3,nexternal), q(0:3,nexternal)
2056
521
      double precision mass_i, mass_j, mass_ij
2057
522
      integer ip, jp
2058
523
C-----Local
2059
524
      integer i, j, nu
2060
525
      double precision qtot(0:3), qsq, pij(0:3), qij(0:3), qi(0:3), qj(0:3)
2061
526
      double precision etot, ztot
2062
527
      double precision msq_ij, msq_k, msq_i, msq_j, resc
2063
528
      double precision a, b
2064
529
      double precision dot, threedot
2065
530
      external dot, threedot
2066
531
 
2067
532
C--------------- 
2068
533
C     BEGIN CODE                                                                    
2069
534
C---------------    
2070
535
      msq_i=mass_i*mass_i
2071
536
      msq_j=mass_j*mass_j
2072
537
      msq_ij=mass_ij*mass_ij
2073
538
2074
539
c ----reconstruct the momentum pij
2075
540
      do nu=0,3
2076
541
        pij(nu)=p(nu,ip) + p(nu,jp)
2077
542
      enddo
2078
543
2079
544
c and go to the pij rest frame
2080
545
      call invboostx(p(0,ip), pij, qi)
2081
546
      call invboostx(p(0,jp), pij, qj)
2082
547
2083
548
C check that momenta are back to back
2084
549
      do nu=1,3
2085
550
        if (dabs(qi(nu)+qj(nu))/(qi(0)+qj(0)).gt.1d-4) then
2086
551
          write(*,*) 'ERROR IN TRANSFORM_OS_INIT: QI QJ not B2B'
2087
552
          write(*,*) 'QI', qi, dsqrt(dot(qi,qi)), mass_i
2088
553
          write(*,*) 'QJ', qj, dsqrt(dot(qj,qj)), mass_j
2089
554
          stop
2090
555
        endif
2091
556
      enddo
2092
557
2093
558
C putting qi, qj on shell amounts in rescaling their spatial 
2094
559
c  component by a factor 
2095
560
      resc = dsqrt(
2096
561
     $ (((msq_ij - msq_i - msq_j)/2d0)**2 - msq_i * msq_j) / 
2097
562
     $ (threedot(qi,qi) * msq_ij))
2098
563
      do nu=1,3
2099
564
        qi(nu) = qi(nu) * resc
2100
565
        qj(nu) = qj(nu) * resc
2101
566
      enddo
2102
567
      qi(0) = dsqrt(msq_i + threedot(qi,qi))
2103
568
      qj(0) = dsqrt(msq_j + threedot(qj,qj))
2104
569
2105
570
c check that the sum of the energies is mij
2106
571
      if (dabs((mass_ij - qi(0) - qj(0))/mass_ij).gt.1d-4) then
2107
572
        write(*,*) 'ERROR IN TRANSFORM_OS_INIT: QI QJ not on shell'
2108
573
        write(*,*) 'MIJ', mass_ij
2109
574
        write(*,*) 'QI', qi, dsqrt(dot(qi,qi)), mass_i
2110
575
        write(*,*) 'QJ', qj, dsqrt(dot(qj,qj)), mass_j
2111
576
        stop
2112
577
      endif
2113
578
      
2114
579
c ----qij has the same spatial components of pij
2115
580
      do nu=1,3
2116
581
        qij(nu) = pij(nu)
2117
582
      enddo
2118
583
      qij(0) = dsqrt(msq_ij + threedot(qij,qij))
2119
584
2120
585
C     boost qi,qj in the lab frame
2121
586
      call boostx(qi, qij, q(0,ip))
2122
587
      call boostx(qj, qij, q(0,jp))
2123
588
2124
589
C check that momenta sum to qij
2125
590
      do nu=0,3
2126
591
        if (dabs(q(nu,ip)+q(nu,jp)-qij(nu))/(qij(0)).gt.1d-4) then
2127
592
          write(*,*) 'ERROR IN TRANSFORM_OS_INIT: QIJ not conserved'
2128
593
          write(*,*) 'QI', q(0,ip), q(1,ip), q(2,ip), q(3,ip)
2129
594
          write(*,*) 'QJ', q(0,jp), q(1,jp), q(2,jp), q(3,jp)
2130
595
          write(*,*) 'QIJ', qij
2131
596
          stop
2132
597
        endif
2133
598
      enddo
2134
599
2135
600
C-----Keep the rest of the FS momenta without any changes
2136
601
      etot = 0d0
2137
602
      ztot = 0d0
2138
603
      do j=nincoming+1,nexternal
2139
604
        do nu=0,3
2140
605
          if(j.ne.ip .and. j.ne.jp)then
2141
606
            q(nu,j) = p(nu,j)
2142
607
          endif
2143
608
        enddo
2144
609
        etot = etot + q(0,j)
2145
610
        ztot = ztot + q(3,j)
2146
611
      enddo
2147
612
2148
613
C initial state momenta: one knows the sum and the difference of
2149
614
C energies (sum of z components of FS momenta)
2150
615
      q(0,1) = (etot + ztot)/2d0
2151
616
      q(1,1) = 0d0
2152
617
      q(2,1) = 0d0
2153
618
      q(3,1) = dsign(q(0,1), p(3,1)) 
2154
619
2155
620
      q(0,2) = (etot - ztot)/2d0
2156
621
      q(1,2) = 0d0
2157
622
      q(2,2) = 0d0
2158
623
      q(3,2) = dsign(q(0,2), p(3,2)) 
2159
624
2160
625
      ! check the momenta before returning
2161
626
      call OS_check_momenta(p, q, ip, jp)
2162
627
2163
628
      return     
2164
629
      end
2165
630
2166
631
2167
632
      subroutine OS_check_momenta(p, q, ip, jp)
2168
633
      ! performs some consistency checks on the momenta
2169
634
      implicit none
2170
635
      include 'nexternal.inc'
2171
636
      double precision p(0:3,nexternal), q(0:3,nexternal)
2172
637
      integer ip, jp
2173
638
      double precision a, b
2174
639
      integer i, j
2175
640
      double precision dot
2176
641
2177
642
      if (nincoming.ne.2) then
2178
643
        write(*,*) 'ERROR IN OS_CHECK_MOMENTA:, nincoming != 2 not'//
2179
644
     $   ' implemented', nincoming
2180
645
        stop
2181
646
      endif
2182
647
2183
648
C--------mass shell conditions
2184
649
      if (dabs(dot(q(0,ip),q(0,ip))-dot(p(0,ip),p(0,ip)))
2185
650
     $ .gt. 1d-3 * max(dot(p(0,ip),p(0,ip)), 1d0)) then
2186
651
        write(*,*) 'ERROR IN OS_CHECK_MOMENTA: I NOT KEPT ON SHELL'
2187
652
        write(*,*) 'MSQ before', dot(p(0,ip),p(0,ip))
2188
653
        write(*,*) 'MSQ after ', dot(q(0,ip),q(0,ip))
2189
654
        stop
2190
655
      endif
2191
656
      if (dabs(dot(q(0,jp),q(0,jp))-dot(p(0,jp),p(0,jp)))
2192
657
     $ .gt. 1d-3 * max(dot(p(0,jp),p(0,jp)), 1d0)) then
2193
658
        write(*,*) 'ERROR IN OS_CHECK_MOMENTA: J NOT KEPT ON SHELL'
2194
659
        write(*,*) 'MSQ before', dot(p(0,jp),p(0,jp))
2195
660
        write(*,*) 'MSQ after ', dot(q(0,jp),q(0,jp))
2196
661
        stop
2197
662
      endif
2198
663
2199
664
C--------momentum conservation
2200
665
      do i = 0,3
2201
666
        a = 0d0
2202
667
        b = 0d0
2203
668
        do j = 1, nexternal
2204
669
          b = max(b, dabs(q(i,j)))
2205
670
          if (j.le.nincoming) then
2206
671
            a = a - q(i,j)
2207
672
          else
2208
673
            a = a + q(i,j)
2209
674
          endif
2210
675
        enddo
2211
676
        if (dabs(a)/b.gt.1d-6) then
2212
677
          write(*,*) 'ERROR IN OS_CHECK_MOMENTA: MOMENTUM CONSERVATION',
2213
678
     $      i, dabs(a), b
2214
679
          do j = 1, nexternal
2215
680
            write(*,*) q(0,j), q(1,j), q(2,j), q(3,j), dsqrt(dot(q(0,j), q(0,j)))
2216
681
          enddo
2217
682
          stop
2218
683
        endif
2219
684
      enddo
2220
685
2221
686
      return
2222
687
      end
2223
688
2224
689
2225
690
      subroutine transform_os_ident(p,q)
2226
691
************************************************************************
2227
692
*     Authors: Marco Zaro                                              *
2228
693
*     No reshuffling is performed in this case                         *
2229
694
************************************************************************
2230
695
      implicit none 
2231
696
      include 'nexternal.inc'
2232
697
C-----Arguments      
2233
698
      double precision p(0:3,nexternal), q(0:3,nexternal)
2234
699
2235
700
      integer i,j
2236
701
2237
702
      do j = 1, nexternal
2238
703
        do i = 0, 3
2239
704
          q(i,j) = p(i,j)
2240
705
        enddo
2241
706
      enddo 
2242
707
2243
708
      return
2244
709
      end
2245
710
2246
711
2247
712
      subroutine generate_momenta_shat(shat, mom_perm, p_reord)
2248
713
C a wrapper to the standard generate_momneta subroutine
2249
714
      implicit none
2250
715
      include 'nexternal.inc'
2251
716
      double precision shat
2252
717
      integer mom_perm(nexternal)
2253
718
      double precision p_reord(0:3,nexternal)
2254
719
      double precision p(0:3,nexternal), wgt
2255
720
      double precision p_dum(0:3,nexternal)
2256
721
      include 'mint.inc'
2257
722
      ! stuff for the on-shell reshuffling
2258
723
      double precision x(99)
2259
724
      common /to_reshuffle/x
2260
725
      integer         ndim
2261
726
      common/tosigint/ndim
2262
727
      integer i, j
2263
728
      logical force_shat
2264
729
      double precision shat_genps
2265
730
      common/to_shat_genps/force_shat, shat_genps
2266
731
2267
732
      integer i_fks,j_fks
2268
733
      common/fks_indices/i_fks,j_fks
2269
734
      double precision xi_i_fks_ev,y_ij_fks_ev
2270
735
      double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
2271
736
      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
2272
737
      double precision sumdot
2273
738
! switch flags for soft/coll tests off in any case
2274
739
      logical softtest,colltest,softtest_save,colltest_save
2275
740
      common/sctests/softtest,colltest
2276
741
2277
742
      double precision shat_target, shat_in_up, shat_in_down
2278
743
      double precision shat_out_up, shat_out_down, shat_out_half
2279
744
      integer icount
2280
745
      double precision shat_tolerance
2281
746
      parameter(shat_tolerance=1d-8)
2282
747
2283
748
      shat_target=shat
2284
749
      ! reduce shat when j is initial state
2285
750
      if (j_fks.lt.nincoming) shat_target=shat*(1d0-xi_i_fks_ev)
2286
751
      shat_genps=shat_target
2287
752
2288
753
      force_shat=.true.
2289
754
      colltest_save=colltest
2290
755
      softtest_save=softtest
2291
756
      colltest=.false.
2292
757
      softtest=.false.
2293
758
2294
759
      call generate_momenta(ndim,iconfig,wgt,x,p)
2295
760
      ! this loop is necessary because, the phase space always generate
2296
761
      ! the born shat, and in the case of ISR (j<=2) the real-emission
2297
762
      ! shat also depends on xi_fks, which also changes (changing shat_born
2298
763
      ! changes the boundaries of xi)
2299
764
      ! normally, just a couple (2) iterations are enough to get shat
2300
765
      ! accurage at less than 1e-8 relative accuracy 
2301
766
      do while (abs(shat-sumdot(p(0,1),p(0,2),1d0))/shat.gt.
2302
767
     $            shat_tolerance)
2303
768
        shat_genps = shat_genps*shat/sumdot(p(0,1),p(0,2),1d0)
2304
769
        call generate_momenta(ndim,iconfig,wgt,x,p)
2305
770
      enddo
2306
771
2307
772
C remember to reorder the momenta, as the order of particles in the FS
2308
773
C of the OS counterterm may differ from the one in the rest of the code
2309
774
      do j = 1, nexternal
2310
775
        do i = 0, 3
2311
776
          p_reord(i, j) = p(i, mom_perm(j))
2312
777
        enddo
2313
778
      enddo
2314
779
2315
780
      ! the second call is to avoid border effects
2316
781
      force_shat=.false.
2317
782
      colltest=colltest_save
2318
783
      softtest=softtest_save
2319
784
      call generate_momenta(ndim,iconfig,wgt,x,p_dum)
2320
785
      ! also call set_cms_stuff
2321
786
      call set_cms_stuff(-100)
2322
787
      return
2323
788
      end 
2324
789
2325
790
          
2326
791
2327
792
      subroutine invboostx(p,q , pboost)
2328
793
c
2329
794
c This subroutine performs the Lorentz boost of a four-momentum.  The
2330
795
c momenta p and q are assumed to be given in the same frame.pboost is
2331
796
c the momentum p boosted to the q rest frame.  q must be a
2332
797
c timelike momentum.
2333
798
c it is the inverse of boostx
2334
799
c
2335
800
c input:
2336
801
c       real    p(0:3)         : four-momentum p in the same frame as q
2337
802
c       real    q(0:3)         : four-momentum q 
2338
803
c
2339
804
c output:
2340
805
c       real    pboost(0:3)    : four-momentum p in the boosted frame
2341
806
c
2342
807
      implicit none
2343
808
      double precision p(0:3),q(0:3),pboost(0:3),pq,qq,m,lf
2344
809
2345
810
      double precision rZero
2346
811
      parameter( rZero = 0.0d0 )
2347
812
2348
813
c#ifdef HELAS_CHECK
2349
814
c      integer stdo
2350
815
c      parameter( stdo = 6 )
2351
816
c      double precision pp
2352
817
c#endif
2353
818
c
2354
819
      qq = q(1)**2+q(2)**2+q(3)**2
2355
820
2356
821
c#ifdef HELAS_CHECK
2357
822
c      if (abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero) then
2358
823
c         write(stdo,*)
2359
824
c     &        ' helas-error : p(0:3) in boostx is zero momentum'
2360
825
c      endif
2361
826
c      if (abs(q(0))+qq.eq.rZero) then
2362
827
c         write(stdo,*)
2363
828
c     &        ' helas-error : q(0:3) in boostx is zero momentum'
2364
829
c      endif
2365
830
c      if (p(0).le.rZero) then
2366
831
c         write(stdo,*)
2367
832
c     &        ' helas-warn  : p(0:3) in boostx has not positive energy'
2368
833
c         write(stdo,*)
2369
834
c     &        '             : p(0) = ',p(0)
2370
835
c      endif
2371
836
c      if (q(0).le.rZero) then
2372
837
c         write(stdo,*)
2373
838
c     &        ' helas-error : q(0:3) in boostx has not positive energy'
2374
839
c         write(stdo,*)
2375
840
c     &        '             : q(0) = ',q(0)
2376
841
c      endif
2377
842
c      pp=p(0)**2-p(1)**2-p(2)**2-p(3)**2
2378
843
c      if (pp.lt.rZero) then
2379
844
c         write(stdo,*)
2380
845
c     &        ' helas-warn  : p(0:3) in boostx is spacelike'
2381
846
c         write(stdo,*)
2382
847
c     &        '             : p**2 = ',pp
2383
848
c      endif
2384
849
c      if (q(0)**2-qq.le.rZero) then
2385
850
c         write(stdo,*)
2386
851
c     &        ' helas-error : q(0:3) in boostx is not timelike'
2387
852
c         write(stdo,*)
2388
853
c     &        '             : q**2 = ',q(0)**2-qq
2389
854
c      endif
2390
855
c      if (qq.eq.rZero) then
2391
856
c         write(stdo,*)
2392
857
c     &   ' helas-warn  : q(0:3) in boostx has zero spacial components'
2393
858
c      endif
2394
859
c#endif
2395
860
2396
861
      if ( qq.ne.rZero ) then
2397
862
         pq = p(1)*q(1)+p(2)*q(2)+p(3)*q(3)
2398
863
         m = dsqrt(max(q(0)**2-qq,1d-99))
2399
864
         lf = (-(q(0)-m)*pq/qq+p(0))/m
2400
865
         pboost(0) = (p(0)*q(0)-pq)/m
2401
866
         pboost(1) =  p(1)-q(1)*lf
2402
867
         pboost(2) =  p(2)-q(2)*lf
2403
868
         pboost(3) =  p(3)-q(3)*lf
2404
869
      else
2405
870
         pboost(0) = p(0)
2406
871
         pboost(1) = p(1)
2407
872
         pboost(2) = p(2)
2408
873
         pboost(3) = p(3)
2409
874
      endif
2410
875
c
2411
876
      return
2412
877
      end
2413
878
2414
879
2415
880
      subroutine get_pdf_flux_ratio(p, q, pdf_ratio, flux_ratio)
2416
881
      ! p is off-shell, q is the reshuffled (on-shell) momentum
2417
882
      implicit none
2418
883
      include 'nexternal.inc'
2419
884
      double precision p(0:3,nexternal), q(0:3,nexternal)
2420
885
      double precision pdf_ratio, flux_ratio
2421
886
      include 'run.inc' ! to acces the Bjorken x's
2422
887
      double precision xbk_save(2)
2423
888
      double precision tiny
2424
889
      parameter (tiny=1e-6)
2425
890
      double precision dot, dlum
2426
891
2427
892
      pdf_ratio = 0d0
2428
893
      flux_ratio = 0d0
2429
894
2430
895
C if initial momenta are identical, just return 1. for both
2431
896
      if (abs(p(0,1)-q(0,1)/p(0,1)+q(0,1)).lt.tiny.and.
2432
897
     $    abs(p(0,2)-q(0,2)/p(0,2)+q(0,2)).lt.tiny) then
2433
898
        pdf_ratio = 1d0
2434
899
        flux_ratio = 1d0
2435
900
        return
2436
901
      endif
2437
902
2438
903
c the xbk corresponding to the os kinematics are choosen in order
2439
904
c to have the same ratio as the original ones
2440
905
      xbk_save(1) = xbk(1)
2441
906
      xbk_save(2) = xbk(2)
2442
907
2443
908
C check for the energy in the P_OS array
2444
909
      if (q(0,1).gt.ebeam(1).or.q(0,2).gt.ebeam(2)) goto 999
2445
910
2446
911
c Reweight parton luminosities
2447
912
      pdf_ratio = 1d0 / dlum()
2448
913
2449
914
      xbk(1) = dsqrt(dot(q(0,1), q(0,2)) / dot(p(0,1), p(0,2))) * xbk(1)
2450
915
      xbk(2) = dsqrt(dot(q(0,1), q(0,2)) / dot(p(0,1), p(0,2))) * xbk(2)
2451
916
2452
917
C restore the bjorken x; in any case flux_ratio is zero, therefore 
2453
918
C the whole on-shell subtraction term will be set to zero as well
2454
919
      if (xbk(1).gt.1d0.or.xbk(2).gt.1d0) goto 999
2455
920
2456
921
      flux_ratio = (xbk_save(1)*xbk_save(2)) / (xbk(1)*xbk(2))
2457
922
      pdf_ratio = pdf_ratio * dlum()
2458
923
2459
924
  999 continue
2460
925
C finally, restore the bjorken X's to the original values
2461
926
      xbk(1) = xbk_save(1)
2462
927
      xbk(2) = xbk_save(2)
2463
928
      return
2464
929
      end
2465
930
2466
931
     
2467
932
      subroutine get_bw_ratio(p, mom_mass, mom_wdth, idau1, idau2, ibw, bw_ratio) 
2468
933
      ! compute the ratio of BW functions
2469
934
      ! ibw==0, return 1.
2470
935
      ! ibw==1, standard BW
2471
936
      ! ibw==2, running BW
2472
937
      implicit none
2473
938
      include 'nexternal.inc'
2474
939
      double precision p(0:3, nexternal)
2475
940
      double precision mom_mass, mom_wdth
2476
941
      integer idau1, idau2, ibw
2477
942
      double precision bw_ratio
2478
943
2479
944
      double precision m2_reco
2480
945
      double precision sumdot
2481
946
2482
947
      m2_reco = sumdot(p(0,idau1),p(0,idau2),1d0)
2483
948
2484
949
      if (ibw.eq.0) then
2485
950
        ! nothing
2486
951
        bw_ratio = 1d0
2487
952
      elseif (ibw.eq.1) then
2488
953
        ! standard BW
2489
954
        bw_ratio = (mom_mass * mom_wdth)**2 / ((m2_reco - mom_mass**2)**2 + (mom_mass * mom_wdth)**2)
2490
955
      elseif (ibw.eq.2) then
2491
956
        ! running BW
2492
957
        bw_ratio = m2_reco * mom_wdth**2 / ((m2_reco - mom_mass**2)**2 + m2_reco * mom_wdth**2)
2493
958
      else
2494
959
        write(*,*) 'ERROR in get_bw_ratio: not implemented', ibw
2495
960
        stop
2496
961
      endif
2497
962
2498
963
      return
2499
964
      end
2500
965
2501
966
2502
967
2503
968
      subroutine OS_resonance_and_decay_reshuffle(npart,ibeta,decay_tree
2504
969
     $     ,Mbeta,p_after,p_out)
2505
970
c
2506
971
c     Authors: Rikkert Frederix
2507
972
c
2508
973
cc Given the inputs
2509
974
c   npart : number of particles before shower (this includes
2510
975
c     resonances). The number of particles after shower should
2511
976
c     be one greater than this.
2512
977
c   ibeta : the label for the resonance which invariant mass needs to
2513
978
c     be conserved. (Label corresponds to particle in p_after).
2514
979
c   decay_tree : structure of 1->2 splittings for decays of the
2515
980
c     resonance 'ibeta'. Arguments are
2516
981
c       1st argument, 1:2, labels repsective daughter
2517
982
c       2nd argument, -next:-1, mother particle.
2518
983
c       3rd argument, 1:2. If 1=2, daughter is final state OR daughter is
2519
984
c         also s-channel particle but not written in event file. If 1!=2, 1
2520
985
c         gives internal labeling (which is negative), 2 gives label of the
2521
986
c         actual resonance in momenta list (which is positive).
2522
987
c   Mbeta : Inv. mass of beta before shower
2523
988
c   p_after : momenta after one emission (i.e. after shower)
2524
989
c this subroutine returns the p_out momenta, determined by reshuffling
2525
990
c p_after such that the invariant mass of beta becomes the same as Mbeta.
2526
991
      implicit none
2527
992
      include 'nexternal.inc'
2528
993
      integer next
2529
994
      parameter (next=2*nexternal-3)
2530
995
      integer npart,decay_tree(2,-next:-1,2),ibeta,i,j,d
2531
996
      double precision Mbeta,p_after(0:4,next),p_out(0:4
2532
997
     $     ,next),p_after_b(0:4,next),pbetab(0:3)
2533
998
      double precision dot
2534
999
c determine reshuffled momenta of resonance and all particles not part
2535
1000
c of the decay of the resonance (EQ.6).
2536
1001
      call OS_reshuffle_resonance(npart,ibeta,Mbeta,p_after
2537
1002
     $     ,p_out)
2538
1003
c From the resonance momenta and the momenta after shower reshuffle all
2539
1004
c the decay products of the resonace.
2540
1005
c First, boost all decays to restframe of beta (EQ.17-24 and 31).
2541
1006
      do j=1,3
2542
1007
         pbetab(j)=-p_after(j,ibeta)
2543
1008
      enddo
2544
1009
      pbetab(0)=p_after(0,ibeta)
2545
1010
      do i=-1,-next,-1
2546
1011
         do j=1,2
2547
1012
            d=decay_tree(j,i,2)
2548
1013
            if (d.lt.0) cycle ! skip s-channels not in event file
2549
1014
            if (d.eq.0) exit  ! done
2550
1015
            call boostm(p_after(0,d),pbetab,p_after(4,ibeta),p_after_b(0
2551
1016
     $           ,d))
2552
1017
            p_after_b(4,d)=p_after(4,d)
2553
1018
         enddo
2554
1019
      enddo
2555
1020
c loop over the decay_tree to reshuffle all the momenta of the decay
2556
1021
c     products to be consistent with the reshuffled resonance momentum
2557
1022
      call OS_reshuffle_decays(ibeta,decay_tree,p_after_b(0,1),p_out(0,1))
2558
1023
c finally, boost all decays back from beta restframe to c.m. frame
2559
1024
c (EQ.26-27 and 40)
2560
1025
      do i=-1,-next,-1
2561
1026
         do j=1,2
2562
1027
            d=decay_tree(j,i,2)
2563
1028
            if (d.lt.0) cycle ! skip s-channels not in event file
2564
1029
            if (d.eq.0) exit  ! done
2565
1030
            call boostm(p_out(0,d),p_out(0,ibeta),p_out(4,ibeta)
2566
1031
     $           ,p_out(0,d))
2567
1032
         enddo
2568
1033
      enddo
2569
1034
      return
2570
1035
      end
2571
1036
2572
1037
2573
1038
2574
1039
      subroutine OS_reshuffle_decays(ibeta,decay_tree,ppa,ppout)
2575
1040
c Reshuffle the momenta of all the particles in the decay-tree of beta.
2576
1041
      implicit none
2577
1042
      include 'nexternal.inc'
2578
1043
      integer next
2579
1044
      parameter (next=2*nexternal-3)
2580
1045
      integer decay_tree(2,-next:-1,2),ibeta,idec,d1,d2,i,j
2581
1046
      double precision ppa(0:4,next), ppout(0:4,next)
2582
1047
      double precision pa(0:4,-next:next),pout(0:4,-next:next),pa_d(0:4
2583
1048
     $     ,2),pout_d(0:4,2)
2584
1049
c First, fill the all the 'decay' momenta in the 1->2 decays (apart from
2585
1050
c the momentum of beta (hence, stop the loop at -2)). Like in genps_fks,
2586
1051
c these have negative labels in 'pa' and 'pout'. If some of the decays
2587
1052
c of beta are themselves resonances, some of these momenta should be
2588
1053
c identical to momenta already in 'pa'. Simply ignore this here, but
2589
1054
c update the ones in pout below.
2590
1055
      do i = 1, next
2591
1056
        do j = 0, 4
2592
1057
          pa(j, i) = ppa(j, i)
2593
1058
          pout(j, i) = ppout(j, i)
2594
1059
        enddo
2595
1060
      enddo
2596
1061
2597
1062
      do idec=-next,-2
2598
1063
         d1=decay_tree(1,idec,1)
2599
1064
         d2=decay_tree(2,idec,1)
2600
1065
         if (d1.eq.0 .and. d2.eq.0) then
2601
1066
            cycle
2602
1067
         elseif (d1.eq.0 .or. d2.eq.0) then
2603
1068
            write (*,*)'Error #1 in OS_reshuffle_decays: both daughters'/
2604
1069
     $           /' should exist or be absent',idec,d1,d2
2605
1070
            stop 1
2606
1071
         endif
2607
1072
         do j=0,3
2608
1073
            pa(j,idec)=pa(j,d1)+pa(j,d2)
2609
1074
         enddo
2610
1075
         pa(4,idec)=sqrt(pa(0,idec)**2-pa(1,idec)**2-pa(2,idec)**2-pa(3
2611
1076
     $        ,idec)**2)
2612
1077
      enddo
2613
1078
c Now loop over all the 1->2 decays, and do their reshuffling.      
2614
1079
      do idec=-1,-next,-1
2615
1080
         d1=decay_tree(1,idec,1)
2616
1081
         d2=decay_tree(2,idec,1)
2617
1082
         if (d1.eq.0 .or. d2.eq.0) exit ! done when no more decays
2618
1083
         do j=0,4
2619
1084
            pa_d(j,1)=pa(j,d1)
2620
1085
            pa_d(j,2)=pa(j,d2)
2621
1086
         enddo
2622
1087
         if (idec.eq.-1) then   ! decay of beta
2623
1088
            call OS_reshuffle_one_decay(pout(4,ibeta),pa_d,pout_d)
2624
1089
         else
2625
1090
           write (*,*)'Error #2 in OS_reshuffle_decays: not implemented',idec
2626
1091
           stop 1
2627
1092
CC            call OS_reshuffle_decay_of_decay(pout(0,decay_tree(1,idec
2628
1093
CC     $           ,2)),pa_d,pout_d)
2629
1094
         endif
2630
1095
         do j=0,4
2631
1096
            pout(j,d1)=pout_d(j,1)
2632
1097
            pout(j,d2)=pout_d(j,2)
2633
1098
         enddo
2634
1099
c     if some of the decays of beta are themselves resonances to be
2635
1100
c     written in the event file, update those momenta as well.
2636
1101
         do i=1,2
2637
1102
            if (decay_tree(i,idec,1).ne.decay_tree(i,idec,2)) then
2638
1103
               do j=0,4
2639
1104
                  pout(j,decay_tree(i,idec,2))=pout(j,decay_tree(i,idec
2640
1105
     $                 ,1))
2641
1106
               enddo
2642
1107
            endif
2643
1108
         enddo
2644
1109
      enddo
2645
1110
      do i = 1, next
2646
1111
        do j = 0, 4
2647
1112
          ppout(j, i) = pout(j, i)
2648
1113
        enddo
2649
1114
      enddo
2650
1115
2651
1116
      return
2652
1117
      end
2653
1118
2654
1119
2655
1120
      subroutine OS_reshuffle_resonance(npart,ibeta,M,p_a,p_out)
2656
1121
c Given the momenta of beta after shower and the recoil momenta,
2657
1122
c reshuffle the former such that its invariant mass coincides with the
2658
1123
c one from before showering
2659
1124
      implicit none
2660
1125
      include 'nexternal.inc'
2661
1126
      integer npart, ibeta
2662
1127
      integer next
2663
1128
      parameter (next=2*nexternal-3)
2664
1129
      double precision q(0:4),p_rec(0:4),M,p_a(0:4,next),p_out(0:4,next)
2665
1130
      double precision expybst,shybst,chybst,chybstmo,xdir(3)
2666
1131
      integer i,j
2667
1132
      double precision lambda2,rho
2668
1133
      external lambda2,rho
2669
1134
      double precision tiny,vtiny,vvtiny
2670
1135
      parameter (tiny=1d-5,vtiny=1d-8,vvtiny=1d-14)
2671
1136
c determine recoil momentum and c.m. momentum (EQ.2-3)
2672
1137
      if (nincoming.eq.2) then
2673
1138
         do j=0,3
2674
1139
            q(j)=p_a(j,1)+p_a(j,2)
2675
1140
         enddo
2676
1141
         q(4)=q(0)
2677
1142
      else
2678
1143
         do j=0,3
2679
1144
            q(j)=p_a(j,1)
2680
1145
         enddo
2681
1146
         q(4)=q(0)
2682
1147
      endif
2683
1148
c     check that we are in the c.m. frame
2684
1149
      do j=1,3
2685
1150
         if (abs(q(j))/q(0).gt.vvtiny) then
2686
1151
            write (*,*) 'Error #2 in OS_reshuffle_resonance: '/
2687
1152
     $           /'Not in C.M. frame',q
2688
1153
            stop 1
2689
1154
         endif
2690
1155
      enddo
2691
1156
      do j=0,3
2692
1157
         p_rec(j)=q(j)-p_a(j,ibeta)
2693
1158
      enddo
2694
1159
      p_rec(4)=sqrt(max(p_rec(0)**2-p_rec(1)**2-p_rec(2)**2-p_rec(3)**2
2695
1160
     $     ,0d0))
2696
1161
c determine the boost -- EQ.(8):
2697
1162
      if (p_rec(4)/q(0).gt.vtiny) then
2698
1163
         expybst=(p_rec(0)+rho(p_rec))/(2*q(0)*p_rec(4)**2) * (q(4)**2
2699
1164
     $        +p_rec(4)**2-M**2-q(4)**2*lambda2(q(0),p_rec(4),M))
2700
1165
         expybst=1/expybst
2701
1166
      else
2702
1167
         expybst=q(0)*(p_rec(0)+rho(p_rec))/(q(4)**2-M**2)*
2703
1168
     $            (1d0+M**2*p_rec(4)**2/(q(4)**2-M**2)**2)
2704
1169
      endif
2705
1170
      if (expybst.lt.0d0) then
2706
1171
         write (*,*) 'Error #3 in OS_reshuffle_resonance: '/
2707
1172
     $        /'Lorentz boost factor is negative',expybst
2708
1173
         stop 1
2709
1174
      elseif (expybst.lt.vvtiny) then
2710
1175
         write (*,*) 'Error #3 in OS_reshuffle_resonance: '/
2711
1176
     $        /'Lorentz boost too extreme',expybst
2712
1177
         stop 1
2713
1178
      endif
2714
1179
      shybst=(expybst-1/expybst)/2.d0
2715
1180
      chybst=(expybst+1/expybst)/2.d0
2716
1181
      chybstmo=chybst-1.d0
2717
1182
c and the boost direction:
2718
1183
      do j=1,3
2719
1184
         xdir(j)=p_a(j,ibeta)/rho(p_a(0,ibeta))
2720
1185
      enddo
2721
1186
c copy the initial-state momenta as they are
2722
1187
      do i=1,nincoming
2723
1188
        do j=0,3
2724
1189
          p_out(j,i)=p_a(j,i)
2725
1190
        enddo
2726
1191
      enddo
2727
1192
c compute all the new momenta EQ.(6):
2728
1193
      do i=nincoming+1,npart+1
2729
1194
         if (i.eq.ibeta) then
2730
1195
            call boostwdir2(chybst,shybst,chybstmo,xdir,p_rec,
2731
1196
     &           p_out(0,i))
2732
1197
            do j=0,3
2733
1198
               p_out(j,i)=q(j)-p_out(j,i)
2734
1199
            enddo
2735
1200
            p_out(4,i)=sqrt(max(p_out(0,i)**2-p_out(1,i)**2
2736
1201
     $           -p_out(2,i)**2-p_out(3,i)**2,0d0))
2737
1202
         else
2738
1203
c     here we fill all p_out, including the decay products. They will be
2739
1204
c     overwritten later.
2740
1205
            call boostwdir2(chybst,shybst,chybstmo,xdir,p_a(0,i),
2741
1206
     &           p_out(0,i))
2742
1207
            p_out(4,i)=p_a(4,i)
2743
1208
         endif
2744
1209
      enddo
2745
1210
c Check that mass is okay:
2746
1211
      if (abs(p_out(4,ibeta)-M)/(p_out(4,ibeta)+M).gt.tiny) then
2747
1212
         write (*,*) 'Error #1 in OS_reshuffle_resonance: '/
2748
1213
     $        /'masses not identical after reshuffle',p_out(4,ibeta),M
2749
1214
         stop 1
2750
1215
      endif
2751
1216
      return
2752
1217
      end
2753
1218
2754
1219
2755
1220
      subroutine OS_reshuffle_one_decay(M,pa,pout)
2756
1221
c Given the two decay products (after shower) in 'pa', reshuffles them
2757
1222
c to make them consistent with decaying particle with mass 'M' and
2758
1223
c returns them in 'pout'. The two daughter momenta are assumed to be
2759
1224
c back-to-back, i.e., we are in the restframe of their mother.
2760
1225
      implicit none
2761
1226
      double precision M,pa(0:4,2),pout(0:4,2),fac1,fac2
2762
1227
      integer i,j
2763
1228
      double precision lambda2,rho
2764
1229
      external lambda2,rho
2765
1230
      double precision vtiny
2766
1231
      parameter (vtiny=1d-12)
2767
1232
2768
1233
      double precision sumdot
2769
1234
      do j=1,3
2770
1235
         if (abs(pa(j,1)+pa(j,2))/(pa(0,1)+pa(0,2)).gt.vtiny) then
2771
1236
            write (*,*) 'Error #1 in OS_reshuffle_one_decay: '/
2772
1237
     $           /'Not in C.M. frame of mother'
2773
1238
            write (*,*) 'daughter 1:',(pa(i,1),i=0,4)
2774
1239
            write (*,*) 'daughter 1:',(pa(i,2),i=0,4)
2775
1240
            write (*,*) 'sum       :',(pa(i,1)+pa(i,2),i=0,4)
2776
1241
            stop 1
2777
1242
         endif
2778
1243
      enddo
2779
1244
      fac1=(pa(4,1)**2-pa(4,2)**2)/M**2
2780
1245
      fac2=M/2d0*lambda2(M,pa(4,1),pa(4,2))/rho(pa(0,1))
2781
1246
      do j=0,3
2782
1247
         if (j.eq.0) then
2783
1248
            pout(j,1)=M/2d0 *(1d0+fac1)
2784
1249
            pout(j,2)=M/2d0 *(1d0-fac1)
2785
1250
         else
2786
1251
            pout(j,1)=fac2*pa(j,1)
2787
1252
            pout(j,2)=fac2*pa(j,2)
2788
1253
         endif
2789
1254
      enddo
2790
1255
      pout(4,1)=pa(4,1)
2791
1256
      pout(4,2)=pa(4,2)
2792
1257
      return
2793
1258
      end
2794
1259
2795
1260
2796
1261
      subroutine write_momenta(p)
2797
1262
      implicit none
2798
1263
      include 'nexternal.inc'
2799
1264
      double precision p(0:3,nexternal)
2800
1265
      integer i
2801
1266
      do i = 1, nexternal
2802
1267
        write(*,*) i, p(0,i), p(1,i), p(2,i), p(3,i)
2803
1268
      enddo
2804
1269
      return
2805
1270
      end
2806
1271
2807
1272
      subroutine write_momenta4(p)
2808
1273
      implicit none
2809
1274
      include 'nexternal.inc'
2810
1275
      double precision p(0:4,nexternal)
2811
1276
      integer i
2812
1277
      do i = 1, nexternal
2813
1278
        write(*,*) i, p(0,i), p(1,i), p(2,i), p(3,i), p(4,i)
2814
1279
      enddo
2815
1280
      return
2816
1281
      end
2817
0
1282
2818
=== modified file 'madgraph/fks/fks_base.py'
2819
--- madgraph/fks/fks_base.py	2016-02-24 13:44:29 +0000
2820
+++ madgraph/fks/fks_base.py	2018-09-06 19:46:41 +0000
2821
@@ -50,14 +50,16 @@
2822
50
        if not 'OLP' in self.keys():
50
        if not 'OLP' in self.keys():
2823
51
            self['OLP'] = 'MadLoop'
51
            self['OLP'] = 'MadLoop'
2824
52
            self['ncores_for_proc_gen'] = 0
52
            self['ncores_for_proc_gen'] = 0
2825
53
2826
53
    
54
    
2827
54
    def get_sorted_keys(self):
55
    def get_sorted_keys(self):
2828
55
        """Return particle property names as a nicely sorted list."""
56
        """Return particle property names as a nicely sorted list."""
2829
56
        keys = super(FKSMultiProcess, self).get_sorted_keys()
57
        keys = super(FKSMultiProcess, self).get_sorted_keys()
2830
57
        keys += ['born_processes', 'real_amplitudes', 'real_pdgs', 'has_isr', 
58
        keys += ['born_processes', 'real_amplitudes', 'real_pdgs', 'has_isr', 
2832
58
                 'has_fsr', 'OLP', 'ncores_for_proc_gen']
59
                 'has_fsr', 'OLP', 'has_os_div', 'ncores_for_proc_gen']
2833
59
        return keys
60
        return keys
2834
60
61
2835
62
2836
61
    def filter(self, name, value):
63
    def filter(self, name, value):
2837
62
        """Filter for valid leg property values."""
64
        """Filter for valid leg property values."""
2838
63
65
2839
@@ -87,6 +89,7 @@
2840
87
                    "%s is not a valid value for ncores_for_proc_gen " % str(value)
89
                    "%s is not a valid value for ncores_for_proc_gen " % str(value)
2841
88
                                                     
90
                                                     
2842
89
        return super(FKSMultiProcess,self).filter(name, value)
91
        return super(FKSMultiProcess,self).filter(name, value)
2843
92
2844
90
    
93
    
2845
91
    def __init__(self, procdef=None, options={}):
94
    def __init__(self, procdef=None, options={}):
2846
92
        """Initializes the original multiprocess, then generates the amps for the 
95
        """Initializes the original multiprocess, then generates the amps for the 
2847
@@ -94,7 +97,6 @@
2848
94
        Real amplitudes are stored in real_amplitudes according on the pdgs of their
97
        Real amplitudes are stored in real_amplitudes according on the pdgs of their
2849
95
        legs (stored in pdgs, so that they need to be generated only once and then reicycled
98
        legs (stored in pdgs, so that they need to be generated only once and then reicycled
2850
96
        """
99
        """
2851
97
2852
98
        #swhich the other loggers off
100
        #swhich the other loggers off
2853
99
        loggers_off = [logging.getLogger('madgraph.diagram_generation'), 
101
        loggers_off = [logging.getLogger('madgraph.diagram_generation'), 
2854
100
                       logging.getLogger('madgraph.loop_diagram_generation')]
102
                       logging.getLogger('madgraph.loop_diagram_generation')]
2855
@@ -108,6 +110,13 @@
2856
108
            olp = options['OLP']
110
            olp = options['OLP']
2857
109
            del options['OLP']
111
            del options['OLP']
2858
110
112
2859
113
        # on-shell removal
2860
114
        if 'remove_os' in options.keys():
2861
115
            self['remove_os'] = options['remove_os']
2862
116
            del options['remove_os']
2863
117
        else:
2864
118
            self['remove_os'] = False
2865
119
2866
111
        ncores_for_proc_gen = 0
120
        ncores_for_proc_gen = 0
2867
112
        # ncores_for_proc_gen has the following meaning
121
        # ncores_for_proc_gen has the following meaning
2868
113
        #   0 : do things the old way
122
        #   0 : do things the old way
2869
@@ -169,8 +178,7 @@
2870
169
        for i, amp in enumerate(amps):
178
        for i, amp in enumerate(amps):
2871
170
            logger.info("Generating FKS-subtracted matrix elements for born process%s (%d / %d)" \
179
            logger.info("Generating FKS-subtracted matrix elements for born process%s (%d / %d)" \
2872
171
                % (amp['process'].nice_string(print_weighted=False).replace(\
180
                % (amp['process'].nice_string(print_weighted=False).replace(\
2875
172
                                                                 'Process', ''),
181
                    'Process', ''), i + 1, len(amps)))
2874
173
                 i + 1, len(amps)))
2876
174
182
2877
175
            born = FKSProcess(amp, ncores_for_proc_gen = self['ncores_for_proc_gen'])
183
            born = FKSProcess(amp, ncores_for_proc_gen = self['ncores_for_proc_gen'])
2878
176
            self['born_processes'].append(born)
184
            self['born_processes'].append(born)
2879
@@ -185,6 +193,9 @@
2880
185
            for born in self['born_processes']:
193
            for born in self['born_processes']:
2881
186
                for real in born.real_amps:
194
                for real in born.real_amps:
2882
187
                    real.find_fks_j_from_i(born_pdg_list)
195
                    real.find_fks_j_from_i(born_pdg_list)
2883
196
                    if self['remove_os']:
2884
197
                        real.find_os_divergences()
2885
198
2886
188
            if amps:
199
            if amps:
2887
189
                if self['process_definitions'][0].get('NLO_mode') == 'all':
200
                if self['process_definitions'][0].get('NLO_mode') == 'all':
2888
190
                    self.generate_virtuals()
201
                    self.generate_virtuals()
2889
@@ -217,6 +228,8 @@
2890
217
228
2891
218
        self['has_isr'] = any([proc.isr for proc in self['born_processes']])
229
        self['has_isr'] = any([proc.isr for proc in self['born_processes']])
2892
219
        self['has_fsr'] = any([proc.fsr for proc in self['born_processes']])
230
        self['has_fsr'] = any([proc.fsr for proc in self['born_processes']])
2893
231
        self['has_os_div'] = any([real.os_amplitudes for proc in self['born_processes'] for real in proc.real_amps])
2894
232
2895
220
233
2896
221
    def add(self, other):
234
    def add(self, other):
2897
222
        """combines self and other, extending the lists of born/real amplitudes"""
235
        """combines self and other, extending the lists of born/real amplitudes"""
2898
@@ -227,19 +240,23 @@
2899
227
        self['pdgs'].extend(other['pdgs'])
240
        self['pdgs'].extend(other['pdgs'])
2900
228
        self['has_isr'] = self['has_isr'] or other['has_isr']
241
        self['has_isr'] = self['has_isr'] or other['has_isr']
2901
229
        self['has_fsr'] = self['has_fsr'] or other['has_fsr']
242
        self['has_fsr'] = self['has_fsr'] or other['has_fsr']
2902
243
        self['has_os_div'] = self['has_os_div'] or other['has_os_div']
2903
230
        self['OLP'] = other['OLP']
244
        self['OLP'] = other['OLP']
2904
231
        self['ncores_for_proc_gen'] = other['ncores_for_proc_gen']
245
        self['ncores_for_proc_gen'] = other['ncores_for_proc_gen']
2905
232
246
2906
247
2907
233
    def get_born_amplitudes(self):
248
    def get_born_amplitudes(self):
2908
234
        """return an amplitudelist with the born amplitudes"""
249
        """return an amplitudelist with the born amplitudes"""
2909
235
        return diagram_generation.AmplitudeList([born.born_amp \
250
        return diagram_generation.AmplitudeList([born.born_amp \
2910
236
                for born in self['born_processes']])
251
                for born in self['born_processes']])
2911
237
252
2912
253
2913
238
    def get_virt_amplitudes(self):
254
    def get_virt_amplitudes(self):
2914
239
        """return an amplitudelist with the virt amplitudes"""
255
        """return an amplitudelist with the virt amplitudes"""
2915
240
        return diagram_generation.AmplitudeList([born.virt_amp \
256
        return diagram_generation.AmplitudeList([born.virt_amp \
2916
241
                for born in self['born_processes'] if born.virt_amp])
257
                for born in self['born_processes'] if born.virt_amp])
2917
242
258
2918
259
2919
243
    def get_real_amplitudes(self):
260
    def get_real_amplitudes(self):
2920
244
        """return an amplitudelist with the real amplitudes"""
261
        """return an amplitudelist with the real amplitudes"""
2921
245
        return self.get('real_amplitudes')
262
        return self.get('real_amplitudes')
2922
@@ -286,7 +303,7 @@
2923
286
    -- fks_infos (list containing the possible fks configs for a given process
303
    -- fks_infos (list containing the possible fks configs for a given process
2924
287
    -- amplitude 
304
    -- amplitude 
2925
288
    -- is_to_integrate
305
    -- is_to_integrate
2927
289
    -- leg permutation<<REMOVED!.
306
    -- os_amplitudes
2928
290
    """
307
    """
2929
291
    
308
    
2930
292
    def __init__(self, born_proc, leglist, ij, ijglu,
309
    def __init__(self, born_proc, leglist, ij, ijglu,
2931
@@ -295,6 +312,15 @@
2932
295
        Stores the fks informations into the list of dictionaries fks_infos
312
        Stores the fks informations into the list of dictionaries fks_infos
2933
296
        """      
313
        """      
2934
297
        self.fks_infos = []
314
        self.fks_infos = []
2935
315
        # amplitudes needed for the subtraction of on-shell divergences
2936
316
        self.os_amplitudes = []
2937
317
        # id's of the intermediate resonances and of decay products 
2938
318
        # each element is a list [mother, daughter1, daughter2]
2939
319
        self.os_ids = []
2940
320
        # list of diagram numbers for the diagrams which are also
2941
321
        # in the OS counterterms
2942
322
        self.os_diagrams = []
2943
323
2944
298
        for leg in leglist:
324
        for leg in leglist:
2945
299
            if leg.get('fks') == 'i':
325
            if leg.get('fks') == 'i':
2946
300
                i_fks = leg.get('number')
326
                i_fks = leg.get('number')
2947
@@ -350,6 +376,131 @@
2948
350
        return self.amplitude
376
        return self.amplitude
2949
351
377
2950
352
378
2951
379
    def find_os_divergences(self):
2952
380
        """this function looks for possible on shell contributions 
2953
381
        to be removed.
2954
382
        It is called by the FKSMultiProcess, when the remove_os
2955
383
        option is set to True.
2956
384
        In order to be agnostic on mass hierarchies all splittings
2957
385
        1->2 3 are investigated, with m1 != 0, m2, m3 != m1
2958
386
        since these resonances are introduced at the real-emission
2959
387
        level, one must have m2=0 or m3=0
2960
388
        """
2961
389
        model = self.process['model']
2962
390
        # take account of the orders for the on shell processes
2963
391
        weighted_order = self.process['orders']['WEIGHTED']
2964
392
2965
393
        # focus only on final state legs
2966
394
        final_legs = [copy.copy(l) for l in self.process['legs'] if l['state']]
2967
395
        for leg_2 in final_legs:
2968
396
            for leg_3 in [l for l in final_legs if l['number'] > leg_2['number']]:
2969
397
                # one of the two legs must be massless
2970
398
                if not leg_2['massless'] and not leg_3['massless']:
2971
399
                    continue
2972
400
                # prepare the leglist for the 'on shell' process, which should
2973
401
                # not contain leg_2 and leg_3, but should contain their mother particle
2974
402
                # if it exists
2975
403
                other_legs = [copy.copy(l) for l in self.process['legs'] if \
2976
404
                        l != leg_2 and l != leg_3]
2977
405
                assert(len(other_legs) == (len(self.process['legs']) - 2))
2978
406
                leg_2_part = model.get('particle_dict')[leg_2['id']]
2979
407
                leg_3_part = model.get('particle_dict')[leg_3['id']]
2980
408
                interactions = [inte for inte in model.get('interaction_dict').values() \
2981
409
                                if len(inte['particles']) == 3 and \
2982
410
                                    leg_2_part in inte['particles'] and \
2983
411
                                    leg_3_part in inte['particles']]
2984
412
2985
413
                for inte in interactions:
2986
414
                    particles = [copy.copy(p) for p in inte['particles']]
2987
415
                    try:
2988
416
                        particles.remove(leg_2_part)
2989
417
                        particles.remove(leg_3_part)
2990
418
                    except ValueError:
2991
419
                        # this is when leg_2 and leg_3 are the same particle
2992
420
                        # and it appears only once in the interacion, so
2993
421
                        # the interaction has to be skipped
2994
422
                        continue
2995
423
                    leg_1_part = particles[0]
2996
424
                    # check that it is massive and its mass it is different from
2997
425
                    # leg_2 and leg_3
2998
426
                    if leg_1_part['mass'].lower() == 'zero' or \
2999
427
                       leg_1_part['mass'] == leg_2_part['mass'] or \
3000
428
                       leg_1_part['mass'] == leg_3_part['mass']:
3001
429
                        continue
3002
430
                    # this should be the final particle (take the antiparticle as
3003
431
                    # it has to go "into" the interaction)
3004
432
3005
433
                    leg_1 = MG.Leg({'state' : True,
3006
434
                                    'id' : leg_1_part.get_anti_pdg_code(),
3007
435
                                    'number': leg_2['number']})
3008
436
3009
437
                    os_legs = [copy.copy(l) for l in other_legs]
3010
438
                    os_legs.insert(leg_2['number'] - 1, leg_1)
3011
439
                    assert(len(os_legs) == (len(self.process['legs']) - 1))
3012
440
                    # count the occurences of leg 1 in the final state legs
3013
441
                    # only one of them has to be decayed
3014
442
                    nleg_1 = [l['id'] for l in os_legs].count(leg_1['id'])
3015
443
                    # construct the decay chain and the process
3016
444
                    # definition
3017
445
                    leg_1_decay = MG.Leg({'id': leg_1['id'], 'state': False})
3018
446
                    leg_2_decay = MG.Leg({'id': leg_2['id'], 'state': True})
3019
447
                    leg_3_decay = MG.Leg({'id': leg_3['id'], 'state': True})
3020
448
                    decay_chain_legs = MG.LegList(\
3021
449
                                       [leg_1_decay, leg_2_decay, leg_3_decay])
3022
450
                    decay_chain = MG.Process(\
3023
451
                                  {'model': model,
3024
452
                                   'legs': MG.LegList(decay_chain_legs),
3025
453
                                   'is_decay_chain': True})
3026
454
3027
455
                    # construct the 'trivial' decay chain to be used when leg_1
3028
456
                    # occurs more than once in the final state legs
3029
457
                    leg_1_decayed = MG.Leg({'id': leg_1['id'], 'state': True})
3030
458
                    trivial_decay_chain_legs = MG.LegList(\
3031
459
                                       [leg_1_decay, leg_1_decayed])
3032
460
                    trivial_decay_chain = MG.Process(\
3033
461
                                  {'model': model,
3034
462
                                   'legs': MG.LegList(trivial_decay_chain_legs),
3035
463
                                   'is_decay_chain': True})
3036
464
                    
3037
465
                    decay_chains = MG.ProcessList([decay_chain] + \
3038
466
                                        [trivial_decay_chain] * (nleg_1 - 1))
3039
467
3040
468
                    for leg in os_legs:
3041
469
                        leg['number'] = os_legs.index(leg) + 1
3042
470
                    # the orders in os_procdef refer only to the production process
3043
471
                    # so the orders of the splitting have to be subtracted
3044
472
                    prod_weighted_order = weighted_order - \
3045
473
                            sum([v * model.get('order_hierarchy')[o] \
3046
474
                                 for o, v in inte['orders'].items()])
3047
475
                    os_procdef =  MG.Process(\
3048
476
                                 {'model': model,
3049
477
                                  'legs': MG.LegList(os_legs),
3050
478
                                  'decay_chains': decay_chains,
3051
479
                                  'orders': {'WEIGHTED': prod_weighted_order}})
3052
480
                    # now generate the amplitude. 
3053
481
                    # Do nothing if any InvalidCmd is raised (e.g. charge not conserved)
3054
482
                    # or if no diagrams are there
3055
483
                    # set the logger to CRITICAL in order not to warn about 1 -> 1
3056
484
                    # (trivial) decay chains
3057
485
                    
3058
486
                    loglevel = logging.getLogger('madgraph.diagram_generation').level
3059
487
                    logging.getLogger('madgraph.diagram_generation').setLevel(logging.CRITICAL)
3060
488
                    try:
3061
489
                        os_amp = diagram_generation.DecayChainAmplitude(os_procdef)
3062
490
                    except InvalidCmd:
3063
491
                        continue
3064
492
                    logging.getLogger('madgraph.diagram_generation').setLevel(loglevel)
3065
493
                    
3066
494
                    if not all([amp['diagrams'] for amp in os_amp['amplitudes']]):
3067
495
                        continue
3068
496
                    logger.info('Process %s has been generated for on-shell subtraction'
3069
497
                            % os_procdef.input_string())
3070
498
                    self.os_amplitudes.append(os_amp)
3071
499
                    self.os_ids.append([leg_1['id'], leg_2['id'], leg_3['id']])
3072
500
                    self.os_diagrams.append(fks_common.find_os_diagrams(\
3073
501
                            self.amplitude, [leg_1, leg_2, leg_3]))
3074
502
3075
503
3076
353
    def find_fks_j_from_i(self, born_pdg_list): #test written
504
    def find_fks_j_from_i(self, born_pdg_list): #test written
3077
354
        """Returns a dictionary with the entries i : [j_from_i], if the born pdgs are in 
505
        """Returns a dictionary with the entries i : [j_from_i], if the born pdgs are in 
3078
355
        born_pdg_list"""
506
        born_pdg_list"""
3079
@@ -386,6 +537,7 @@
3080
386
                    'get_leg_i should only be called before combining processes')
537
                    'get_leg_i should only be called before combining processes')
3081
387
        return self.process.get('legs')[self.fks_infos[0]['i'] - 1]
538
        return self.process.get('legs')[self.fks_infos[0]['i'] - 1]
3082
388
539
3083
540
3084
389
    def get_leg_j(self): #test written
541
    def get_leg_j(self): #test written
3085
390
        """Returns leg corresponding to j fks.
542
        """Returns leg corresponding to j fks.
3086
391
        An error is raised if the fks_infos list has more than one entry"""
543
        An error is raised if the fks_infos list has more than one entry"""
3087
@@ -501,8 +653,9 @@
3088
501
653
3089
502
654
3090
503
    def generate_real_amplitudes(self, pdg_list, real_amp_list):
655
    def generate_real_amplitudes(self, pdg_list, real_amp_list):
3093
504
        """generates the real amplitudes for all the real emission processes, using pdgs and real_amps
656
        """generates the real amplitudes for all the real emission processes, 
3094
505
        to avoid multiple generation of the same amplitude"""
657
        using pdgs and real_amps to avoid multiple generation of the same amplitude
3095
658
        """
3096
506
659
3097
507
        for amp in self.real_amps:
660
        for amp in self.real_amps:
3098
508
            try:
661
            try:
3099
@@ -514,7 +667,8 @@
3100
514
667
3101
515
    def combine_real_amplitudes(self):
668
    def combine_real_amplitudes(self):
3102
516
        """combines real emission processes if the pdgs are the same, combining the lists 
669
        """combines real emission processes if the pdgs are the same, combining the lists 
3104
517
        of fks_infos"""
670
        of fks_infos
3105
671
        """
3106
518
        pdgs = []
672
        pdgs = []
3107
519
        real_amps = []
673
        real_amps = []
3108
520
        old_real_amps = copy.copy(self.real_amps)
674
        old_real_amps = copy.copy(self.real_amps)
3109
521
675
3110
=== modified file 'madgraph/fks/fks_common.py'
3111
--- madgraph/fks/fks_common.py	2016-08-16 12:38:17 +0000
3112
+++ madgraph/fks/fks_common.py	2018-09-06 19:46:41 +0000
3113
@@ -47,6 +47,22 @@
3114
47
        return [((leg.get('id'), leg.get('number')), leg.get('number'))]
47
        return [((leg.get('id'), leg.get('number')), leg.get('number'))]
3115
48
48
3116
49
49
3117
50
def find_os_diagrams(amp, legs):
3118
51
    """ return the diagram number of the diagrams which correspond to the production
3119
52
    x decay (legs[0] -> legs[1] -> legs[2]
3120
53
    """
3121
54
    mom, dau1, dau2 = legs
3122
55
    os_diagrams = []
3123
56
3124
57
    for i, diag in enumerate(amp['diagrams']):
3125
58
        for vert in diag['vertices']:
3126
59
            ids = [abs(l['id']) for l in vert['legs']]
3127
60
            if dau1 in vert['legs'] and dau2 in vert['legs'] and abs(mom['id']) in ids:
3128
61
                os_diagrams.append(i)
3129
62
3130
63
    return os_diagrams
3131
64
3132
65
3133
50
def link_rb_configs(born_amp, real_amp, i, j, ij):
66
def link_rb_configs(born_amp, real_amp, i, j, ij):
3134
51
    """finds the real configurations that match the born ones, i.e.  for
67
    """finds the real configurations that match the born ones, i.e.  for
3135
52
    each born configuration, the real configuration that has the ij ->
68
    each born configuration, the real configuration that has the ij ->
3136
53
69
3137
=== modified file 'madgraph/fks/fks_helas_objects.py'
3138
--- madgraph/fks/fks_helas_objects.py	2017-06-16 12:58:46 +0000
3139
+++ madgraph/fks/fks_helas_objects.py	2018-09-06 19:46:41 +0000
3140
@@ -202,10 +202,11 @@
3141
202
    def get_sorted_keys(self):
202
    def get_sorted_keys(self):
3142
203
        """Return particle property names as a nicely sorted list."""
203
        """Return particle property names as a nicely sorted list."""
3143
204
        keys = super(FKSHelasMultiProcess, self).get_sorted_keys()
204
        keys = super(FKSHelasMultiProcess, self).get_sorted_keys()
3145
205
        keys += ['real_matrix_elements', ['has_isr'], ['has_fsr'], 
205
        keys += ['real_matrix_elements', ['has_isr'], ['has_fsr'], ['has_os_div'],
3146
206
                 'used_lorentz', 'used_couplings', 'max_configs', 'max_particles', 'processes']
206
                 'used_lorentz', 'used_couplings', 'max_configs', 'max_particles', 'processes']
3147
207
        return keys
207
        return keys
3148
208
208
3149
209
3150
209
    def filter(self, name, value):
210
    def filter(self, name, value):
3151
210
        """Filter for valid leg property values."""
211
        """Filter for valid leg property values."""
3152
211
212
3153
@@ -213,6 +214,7 @@
3154
213
            if not isinstance(value, helas_objects.HelasMultiProcess):
214
            if not isinstance(value, helas_objects.HelasMultiProcess):
3155
214
                raise self.PhysicsObjectError, \
215
                raise self.PhysicsObjectError, \
3156
215
                        "%s is not a valid list for real_matrix_element " % str(value)                             
216
                        "%s is not a valid list for real_matrix_element " % str(value)                             
3157
217
3158
216
    
218
    
3159
217
    def __init__(self, fksmulti, loop_optimized = False, gen_color =True, decay_ids =[]):
219
    def __init__(self, fksmulti, loop_optimized = False, gen_color =True, decay_ids =[]):
3160
218
        """Initialization from a FKSMultiProcess"""
220
        """Initialization from a FKSMultiProcess"""
3161
@@ -410,12 +412,14 @@
3162
410
412
3163
411
        self['has_isr'] = fksmulti['has_isr']
413
        self['has_isr'] = fksmulti['has_isr']
3164
412
        self['has_fsr'] = fksmulti['has_fsr']
414
        self['has_fsr'] = fksmulti['has_fsr']
3165
415
        self['has_os_div'] = fksmulti['has_os_div']
3166
413
416
3167
414
        logger.info('... Done')
417
        logger.info('... Done')
3168
415
418
3169
416
        for i, logg in enumerate(loggers_off):
419
        for i, logg in enumerate(loggers_off):
3170
417
            logg.setLevel(old_levels[i])
420
            logg.setLevel(old_levels[i])
3171
418
        
421
        
3172
422
3173
419
    def get_used_lorentz(self):
423
    def get_used_lorentz(self):
3174
420
        """Return a list of (lorentz_name, conjugate, outgoing) with
424
        """Return a list of (lorentz_name, conjugate, outgoing) with
3175
421
        all lorentz structures used by this HelasMultiProcess."""
425
        all lorentz structures used by this HelasMultiProcess."""
3176
@@ -429,6 +433,7 @@
3177
429
        return self['used_lorentz']
433
        return self['used_lorentz']
3178
430
434
3179
431
435
3180
436
3181
432
    def get_used_couplings(self):
437
    def get_used_couplings(self):
3182
433
        """Return a list with all couplings used by this
438
        """Return a list with all couplings used by this
3183
434
        HelasMatrixElement."""
439
        HelasMatrixElement."""
3184
@@ -656,9 +661,9 @@
3185
656
                          optimized_output = loop_optimized)
661
                          optimized_output = loop_optimized)
3186
657
            else: 
662
            else: 
3187
658
                self.virt_matrix_element = None
663
                self.virt_matrix_element = None
3188
659
#            self.color_links_info = fksproc.find_color_links()
3189
660
            self.color_links = []
664
            self.color_links = []
3190
661
665
3191
666
3192
662
    def set_color_links(self):
667
    def set_color_links(self):
3193
663
        """this function computes and returns the color links, it should be called
668
        """this function computes and returns the color links, it should be called
3194
664
        after the initialization and the setting of the color basis"""
669
        after the initialization and the setting of the color basis"""
3195
@@ -673,6 +678,18 @@
3196
673
                                    self.born_matrix_element.get('base_amplitude')),
678
                                    self.born_matrix_element.get('base_amplitude')),
3197
674
                                color_links_info)    
679
                                color_links_info)    
3198
675
680
3199
681
3200
682
    def get_os_ids(self):
3201
683
        """Returns the list of the fks infos for all processes in the format
3202
684
        {n_me, pdgs, fks_info}, where n_me is the number of real_matrix_element the configuration
3203
685
        belongs to"""
3204
686
        os_ids = []
3205
687
        for real in self.real_processes:
3206
688
            # append only the mother particle, i.e. the 1st particle in each list of ids
3207
689
            os_ids += [ids[0] for ids in real.os_ids]
3208
690
        return set(os_ids)
3209
691
3210
692
3211
676
    def get_fks_info_list(self):
693
    def get_fks_info_list(self):
3212
677
        """Returns the list of the fks infos for all processes in the format
694
        """Returns the list of the fks infos for all processes in the format
3213
678
        {n_me, pdgs, fks_info}, where n_me is the number of real_matrix_element the configuration
695
        {n_me, pdgs, fks_info}, where n_me is the number of real_matrix_element the configuration
3214
@@ -703,17 +720,21 @@
3215
703
        """the get function references to born_matrix_element"""
720
        """the get function references to born_matrix_element"""
3216
704
        return self.born_matrix_element.get(key)
721
        return self.born_matrix_element.get(key)
3217
705
    
722
    
3218
723
3219
706
    def get_used_lorentz(self):
724
    def get_used_lorentz(self):
3220
707
        """the get_used_lorentz function references to born, reals
725
        """the get_used_lorentz function references to born, reals
3221
708
        and virtual matrix elements"""
726
        and virtual matrix elements"""
3222
709
        lorentz_list = self.born_matrix_element.get_used_lorentz()
727
        lorentz_list = self.born_matrix_element.get_used_lorentz()
3223
710
        for real in self.real_processes:
728
        for real in self.real_processes:
3224
711
            lorentz_list.extend(real.matrix_element.get_used_lorentz())
729
            lorentz_list.extend(real.matrix_element.get_used_lorentz())
3225
730
            for os_real in real.os_matrix_elements:
3226
731
                lorentz_list.extend(os_real.get_used_lorentz())
3227
712
        if self.virt_matrix_element:
732
        if self.virt_matrix_element:
3228
713
            lorentz_list.extend(self.virt_matrix_element.get_used_lorentz())
733
            lorentz_list.extend(self.virt_matrix_element.get_used_lorentz())
3229
714
734
3230
715
        return list(set(lorentz_list))
735
        return list(set(lorentz_list))
3231
716
    
736
    
3232
737
3233
717
    def get_used_couplings(self):
738
    def get_used_couplings(self):
3234
718
        """the get_used_couplings function references to born, reals
739
        """the get_used_couplings function references to born, reals
3235
719
        and virtual matrix elements"""
740
        and virtual matrix elements"""
3236
@@ -721,6 +742,9 @@
3237
721
        for real in self.real_processes:
742
        for real in self.real_processes:
3238
722
            coupl_list.extend([c for c in\
743
            coupl_list.extend([c for c in\
3239
723
                        real.matrix_element.get_used_couplings()])
744
                        real.matrix_element.get_used_couplings()])
3240
745
            for os_real in real.os_matrix_elements:
3241
746
                coupl_list.extend([c for c in\
3242
747
                        os_real.get_used_couplings()])
3243
724
        if self.virt_matrix_element:
748
        if self.virt_matrix_element:
3244
725
            coupl_list.extend(self.virt_matrix_element.get_used_couplings())
749
            coupl_list.extend(self.virt_matrix_element.get_used_couplings())
3245
726
        return coupl_list    
750
        return coupl_list    
3246
@@ -735,6 +759,7 @@
3247
735
            nexternal += 1
759
            nexternal += 1
3248
736
        return (nexternal, ninitial)
760
        return (nexternal, ninitial)
3249
737
    
761
    
3250
762
3251
738
    def __eq__(self, other):
763
    def __eq__(self, other):
3252
739
        """the equality between two FKSHelasProcesses is defined up to the 
764
        """the equality between two FKSHelasProcesses is defined up to the 
3253
740
        color links"""
765
        color links"""
3254
@@ -756,16 +781,20 @@
3255
756
        else: 
781
        else: 
3256
757
            return False
782
            return False
3257
758
    
783
    
3259
759
    def add_process(self, other): #test written, ppwj
784
3260
785
    def add_process(self, other): #test written
3261
760
        """adds processes from born and reals of other to itself. Note that 
786
        """adds processes from born and reals of other to itself. Note that 
3262
761
        corresponding real processes may not be in the same order. This is 
787
        corresponding real processes may not be in the same order. This is 
3263
762
        taken care of by constructing the list of self_reals."""
788
        taken care of by constructing the list of self_reals."""
3264
789
        # add the born and the virtuals
3265
763
        self.born_matrix_element.get('processes').extend(
790
        self.born_matrix_element.get('processes').extend(
3266
764
                other.born_matrix_element.get('processes'))
791
                other.born_matrix_element.get('processes'))
3267
765
        if self.virt_matrix_element and other.virt_matrix_element:
792
        if self.virt_matrix_element and other.virt_matrix_element:
3268
766
            self.virt_matrix_element.get('processes').extend(
793
            self.virt_matrix_element.get('processes').extend(
3269
767
                    other.virt_matrix_element.get('processes'))
794
                    other.virt_matrix_element.get('processes'))
3270
768
        self_reals = [real.matrix_element for real in self.real_processes]
795
        self_reals = [real.matrix_element for real in self.real_processes]
3271
796
        #
3272
797
        # now the real-emission processes
3273
769
        for oth_real in other.real_processes:
798
        for oth_real in other.real_processes:
3274
770
            this_real = self.real_processes[self_reals.index(oth_real.matrix_element)]
799
            this_real = self.real_processes[self_reals.index(oth_real.matrix_element)]
3275
771
            #need to store pdg lists rather than processes in order to keep mirror processes different
800
            #need to store pdg lists rather than processes in order to keep mirror processes different
3276
@@ -776,10 +805,19 @@
3277
776
                if oth_pdgs not in this_pdgs:
805
                if oth_pdgs not in this_pdgs:
3278
777
                    this_real.matrix_element['processes'].append(oth_proc)
806
                    this_real.matrix_element['processes'].append(oth_proc)
3279
778
                    this_pdgs.append(oth_pdgs)
807
                    this_pdgs.append(oth_pdgs)
3280
808
            #
3281
809
            # finally on-shell subtractions
3282
810
            for oth_on_shell in oth_real.os_matrix_elements:
3283
811
                this_on_shell = this_real.os_matrix_elements[
3284
812
                        this_real.os_matrix_elements.index(oth_on_shell)]
3285
813
                this_pdgs = [[leg['id'] for leg in proc['legs']] \
3286
814
                    for proc in this_on_shell['processes']]
3287
815
                for oth_proc in oth_on_shell['processes']:
3288
816
                    oth_pdgs = [leg['id'] for leg in oth_proc['legs']]
3289
817
                    if oth_pdgs not in this_pdgs:
3290
818
                        this_on_shell['processes'].append(oth_proc)
3291
819
                        this_pdgs.append(oth_pdgs)
3292
779
820
3293
780
 #                       if p not in self.real_processes[\
3294
781
 #                       self_reals.index(oth_real.matrix_element)].matrix_element['processes']])
3295
782
            
3296
783
    
821
    
3297
784
class FKSHelasRealProcess(object): #test written
822
class FKSHelasRealProcess(object): #test written
3298
785
    """class to generate the Helas calls for a FKSRealProcess
823
    """class to generate the Helas calls for a FKSRealProcess
3299
@@ -792,7 +830,8 @@
3300
792
    -- fks_j_from_i
830
    -- fks_j_from_i
3301
793
    -- matrix element
831
    -- matrix element
3302
794
    -- is_to_integrate
832
    -- is_to_integrate
3304
795
    -- leg permutation<<REMOVED"""
833
    -- os_matrix_elements
3305
834
    """
3306
796
    
835
    
3307
797
    def __init__(self, fksrealproc=None, real_me_list = [], real_amp_list =[], **opts):
836
    def __init__(self, fksrealproc=None, real_me_list = [], real_amp_list =[], **opts):
3308
798
        """constructor, starts from a fksrealproc and then calls the
837
        """constructor, starts from a fksrealproc and then calls the
3309
@@ -800,13 +839,16 @@
3310
800
        Sets i/j fks and the permutation.
839
        Sets i/j fks and the permutation.
3311
801
        real_me_list and real_amp_list are the lists of pre-generated matrix elements in 1-1 
840
        real_me_list and real_amp_list are the lists of pre-generated matrix elements in 1-1 
3312
802
        correspondance with the amplitudes"""
841
        correspondance with the amplitudes"""
3314
803
        
842
3315
843
        self.os_matrix_elements = []
3316
804
        if fksrealproc != None:
844
        if fksrealproc != None:
3317
805
            self.isfinite = False
845
            self.isfinite = False
3318
806
            self.colors = fksrealproc.colors
846
            self.colors = fksrealproc.colors
3319
807
            self.charges = fksrealproc.charges
847
            self.charges = fksrealproc.charges
3320
808
            self.fks_infos = fksrealproc.fks_infos
848
            self.fks_infos = fksrealproc.fks_infos
3321
809
            self.is_to_integrate = fksrealproc.is_to_integrate
849
            self.is_to_integrate = fksrealproc.is_to_integrate
3322
850
            self.os_ids = fksrealproc.os_ids
3323
851
            self.os_diagrams = fksrealproc.os_diagrams
3324
810
852
3325
811
            # real_me_list is a list in the old NLO generation mode;
853
            # real_me_list is a list in the old NLO generation mode;
3326
812
            # in the new one it is a matrix element
854
            # in the new one it is a matrix element
3327
@@ -823,7 +865,6 @@
3328
823
                self.matrix_element = real_me_list
865
                self.matrix_element = real_me_list
3329
824
866
3330
825
            else:
867
            else:
3331
826
3332
827
                if real_me_list and real_amp_list:
868
                if real_me_list and real_amp_list:
3333
828
                    self.matrix_element = copy.deepcopy(real_me_list[real_amp_list.index(fksrealproc.amplitude)])
869
                    self.matrix_element = copy.deepcopy(real_me_list[real_amp_list.index(fksrealproc.amplitude)])
3334
829
                    self.matrix_element['processes'] = copy.deepcopy(self.matrix_element['processes'])
870
                    self.matrix_element['processes'] = copy.deepcopy(self.matrix_element['processes'])
3335
@@ -837,9 +878,14 @@
3336
837
                    self.matrix_element.set('color_matrix',
878
                    self.matrix_element.set('color_matrix',
3337
838
                                     color_amp.ColorMatrix(
879
                                     color_amp.ColorMatrix(
3338
839
                                        self.matrix_element.get('color_basis')))
880
                                        self.matrix_element.get('color_basis')))
3339
881
                for os_amp in fksrealproc.os_amplitudes:
3340
882
                    me = helas_objects.HelasDecayChainProcess(os_amp)
3341
883
                    self.os_matrix_elements.append(me.combine_decay_chain_processes()[0])
3342
840
            #self.fks_j_from_i = fksrealproc.find_fks_j_from_i()
884
            #self.fks_j_from_i = fksrealproc.find_fks_j_from_i()
3343
885
           
3344
841
            self.fks_j_from_i = fksrealproc.fks_j_from_i
886
            self.fks_j_from_i = fksrealproc.fks_j_from_i
3345
842
887
3346
888
3347
843
    def get_nexternal_ninitial(self):
889
    def get_nexternal_ninitial(self):
3348
844
        """Refers to the matrix_element function"""
890
        """Refers to the matrix_element function"""
3349
845
        return self.matrix_element.get_nexternal_ninitial()
891
        return self.matrix_element.get_nexternal_ninitial()
3350
846
892
3351
=== modified file 'madgraph/interface/amcatnlo_interface.py'
3352
--- madgraph/interface/amcatnlo_interface.py	2018-08-28 15:23:19 +0000
3353
+++ madgraph/interface/amcatnlo_interface.py	2018-09-06 19:46:41 +0000
3354
@@ -455,6 +455,7 @@
3355
455
        #self.options['group_subprocesses'] = 'False'
455
        #self.options['group_subprocesses'] = 'False'
3356
456
        collect_mirror_procs = False
456
        collect_mirror_procs = False
3357
457
        ignore_six_quark_processes = self.options['ignore_six_quark_processes']
457
        ignore_six_quark_processes = self.options['ignore_six_quark_processes']
3358
458
        remove_os = self.options['remove_os']
3359
458
        if ',' in line:
459
        if ',' in line:
3360
459
            myprocdef, line = mg_interface.MadGraphCmd.extract_decay_chain_process(self,line)
460
            myprocdef, line = mg_interface.MadGraphCmd.extract_decay_chain_process(self,line)
3361
460
            if myprocdef.are_decays_perturbed():
461
            if myprocdef.are_decays_perturbed():
3362
@@ -494,6 +495,7 @@
3363
494
        # this is the options dictionary to pass to the FKSMultiProcess
495
        # this is the options dictionary to pass to the FKSMultiProcess
3364
495
        fks_options = {'OLP': self.options['OLP'],
496
        fks_options = {'OLP': self.options['OLP'],
3365
496
                       'ignore_six_quark_processes': self.options['ignore_six_quark_processes'],
497
                       'ignore_six_quark_processes': self.options['ignore_six_quark_processes'],
3366
498
                       'remove_os': self.options['remove_os'],
3367
497
                       'ncores_for_proc_gen': self.ncores_for_proc_gen}
499
                       'ncores_for_proc_gen': self.ncores_for_proc_gen}
3368
498
        try:
500
        try:
3369
499
            self._fks_multi_proc.add(fks_base.FKSMultiProcess(myprocdef,fks_options))
501
            self._fks_multi_proc.add(fks_base.FKSMultiProcess(myprocdef,fks_options))
3370
@@ -656,7 +658,7 @@
3371
656
            #_curr_matrix_element is a FKSHelasMultiProcess Object 
658
            #_curr_matrix_element is a FKSHelasMultiProcess Object 
3372
657
            self._fks_directories = []
659
            self._fks_directories = []
3373
658
            proc_charac = self._curr_exporter.proc_characteristic 
660
            proc_charac = self._curr_exporter.proc_characteristic 
3375
659
            for charac in ['has_isr', 'has_fsr', 'has_loops']:
661
            for charac in ['has_isr', 'has_fsr', 'has_loops', 'has_os_div']:
3376
660
                proc_charac[charac] = self._curr_matrix_elements[charac]
662
                proc_charac[charac] = self._curr_matrix_elements[charac]
3377
661
663
3378
662
            # prepare for the generation
664
            # prepare for the generation
3379
663
665
3380
=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
3381
--- madgraph/interface/amcatnlo_run_interface.py	2018-08-27 14:51:23 +0000
3382
+++ madgraph/interface/amcatnlo_run_interface.py	2018-09-06 19:46:41 +0000
3383
@@ -1623,7 +1623,8 @@
3384
1623
            options = options.__dict__
1623
            options = options.__dict__
3385
1624
            self.check_launch(argss, options)
1624
            self.check_launch(argss, options)
3386
1625
1625
3388
1626
        
1626
        self.get_characteristics(pjoin(self.me_dir, 'SubProcesses', 'proc_characteristics'))
3389
1627
3390
1627
        if 'run_name' in options.keys() and options['run_name']:
1628
        if 'run_name' in options.keys() and options['run_name']:
3391
1628
            self.run_name = options['run_name']
1629
            self.run_name = options['run_name']
3392
1629
            # if a dir with the given run_name already exists
1630
            # if a dir with the given run_name already exists
3393
1630
1631
3394
=== modified file 'madgraph/interface/common_run_interface.py'
3395
--- madgraph/interface/common_run_interface.py	2018-07-22 11:40:04 +0000
3396
+++ madgraph/interface/common_run_interface.py	2018-09-06 19:46:41 +0000
3397
@@ -31,7 +31,12 @@
3398
31
import sys
31
import sys
3399
32
import time
32
import time
3400
33
import traceback
33
import traceback
3401
34
<<<<<<< TREE
3402
34
import urllib
35
import urllib
3403
36
=======
3404
37
import copy
3405
38
import sets
3406
39
>>>>>>> MERGE-SOURCE
3407
35
import glob
40
import glob
3408
36
import StringIO
41
import StringIO
3409
37
42
3410
@@ -913,8 +918,14 @@
3411
913
918
3412
914
919
3413
915
            if amcatnlo and not keepwidth:
920
            if amcatnlo and not keepwidth:
3414
921
                # check if one needs to keep track of width (if removing os divergences)
3415
922
3416
916
                # force particle in final states to have zero width
923
                # force particle in final states to have zero width
3417
917
                pids = self.get_pid_final_initial_states()
924
                pids = self.get_pid_final_initial_states()
3418
925
                # keep also track of the particles that go on-shell, for which
3419
926
                # the width has to be treated in a special manner when
3420
927
                # doing os subtraction
3421
928
                os_pids = self.get_pid_onshell()
3422
918
                # check those which are charged under qcd
929
                # check those which are charged under qcd
3423
919
                if not MADEVENT and pjoin(self.me_dir,'bin','internal') not in sys.path:
930
                if not MADEVENT and pjoin(self.me_dir,'bin','internal') not in sys.path:
3424
920
                        sys.path.insert(0,pjoin(self.me_dir,'bin','internal'))
931
                        sys.path.insert(0,pjoin(self.me_dir,'bin','internal'))
3425
@@ -929,22 +940,29 @@
3426
929
940
3427
930
                import ufomodel as ufomodel
941
                import ufomodel as ufomodel
3428
931
                zero = ufomodel.parameters.ZERO
942
                zero = ufomodel.parameters.ZERO
3429
943
                pids.update(os_pids)
3430
932
                no_width = [p for p in ufomodel.all_particles
944
                no_width = [p for p in ufomodel.all_particles
3431
933
                        if (str(p.pdg_code) in pids or str(-p.pdg_code) in pids)
945
                        if (str(p.pdg_code) in pids or str(-p.pdg_code) in pids)
3432
934
                           and p.color != 1 and p.width != zero]
946
                           and p.color != 1 and p.width != zero]
3433
947
3434
935
                done = []
948
                done = []
3435
949
                # for on_shell removal, one needs to store the width of particles
3436
950
                # in some ad-hoc variables
3437
951
                decay_to_keep = []
3438
936
                for part in no_width:
952
                for part in no_width:
3439
937
                    if abs(part.pdg_code) in done:
953
                    if abs(part.pdg_code) in done:
3440
938
                        continue
954
                        continue
3441
939
                    done.append(abs(part.pdg_code))
955
                    done.append(abs(part.pdg_code))
3442
940
                    param = param_card['decay'].get((part.pdg_code,))
956
                    param = param_card['decay'].get((part.pdg_code,))
3443
957
                    if self.proc_characteristics['has_os_div']:
3444
958
                        decay_to_keep.append((part.get('width'), copy.copy(param)))
3445
941
959
3446
942
                    if  param.value != 0:
960
                    if  param.value != 0:
3447
943
                        logger.info('''For gauge cancellation, the width of \'%s\' has been set to zero.'''\
961
                        logger.info('''For gauge cancellation, the width of \'%s\' has been set to zero.'''\
3448
944
                                    % part.name,'$MG:BOLD')
962
                                    % part.name,'$MG:BOLD')
3449
945
                        param.value = 0
963
                        param.value = 0
3450
946
964
3452
947
            param_card.write_inc_file(outfile, ident_card, default)
965
            param_card.write_inc_file(outfile, ident_card, default, decay_to_keep)
3453
948
966
3454
949
    def get_model(self):
967
    def get_model(self):
3455
950
        """return the model related to this process"""
968
        """return the model related to this process"""
3456
@@ -3024,6 +3042,19 @@
3457
3024
3042
3458
3025
        return pids
3043
        return pids
3459
3026
3044
3460
3045
3461
3046
    ############################################################################
3462
3047
    def get_pid_onshell(self):
3463
3048
        """Find the pid of all particles in the intermediate on-sheel partices"""
3464
3049
        pids = set()
3465
3050
        os_ids_lines = [l.strip() for l \
3466
3051
                in open(pjoin(self.me_dir,'SubProcesses', 'os_ids.mg')) if l]
3467
3052
        for l in os_ids_lines:
3468
3053
            ids = l.split()[1:]
3469
3054
            pids.update(set(ids))
3470
3055
3471
3056
        return pids
3472
3057
3473
3027
    ############################################################################
3058
    ############################################################################
3474
3028
    def get_pdf_input_filename(self):
3059
    def get_pdf_input_filename(self):
3475
3029
        """return the name of the file which is used by the pdfset"""
3060
        """return the name of the file which is used by the pdfset"""
3476
3030
3061
3477
=== modified file 'madgraph/interface/madgraph_interface.py'
3478
--- madgraph/interface/madgraph_interface.py	2018-08-10 13:04:34 +0000
3479
+++ madgraph/interface/madgraph_interface.py	2018-09-06 19:46:41 +0000
3480
@@ -1398,7 +1398,8 @@
3481
1398
        if len(args) == 1 and args[0] in ['complex_mass_scheme',\
1398
        if len(args) == 1 and args[0] in ['complex_mass_scheme',\
3482
1399
                                          'loop_optimized_output',\
1399
                                          'loop_optimized_output',\
3483
1400
                                          'loop_color_flows',\
1400
                                          'loop_color_flows',\
3485
1401
                                          'low_mem_multicore_nlo_generation']:
1401
                                          'low_mem_multicore_nlo_generation',\
3486
1402
                                          'remove_os']:
3487
1402
            args.append('True')
1403
            args.append('True')
3488
1403
1404
3489
1404
        if len(args) > 2 and '=' == args[1]:
1405
        if len(args) > 2 and '=' == args[1]:
3490
@@ -1447,7 +1448,7 @@
3491
1447
            if not args[1].isdigit():
1448
            if not args[1].isdigit():
3492
1448
                raise self.InvalidCmd('%s values should be a integer' % args[0])
1449
                raise self.InvalidCmd('%s values should be a integer' % args[0])
3493
1449
1450
3495
1450
        if args[0] in ['loop_optimized_output', 'loop_color_flows', 'low_mem_multicore_nlo_generation']:
1451
        if args[0] in ['loop_optimized_output', 'loop_color_flows', 'low_mem_multicore_nlo_generation', 'remove_os']:
3496
1451
            try:
1452
            try:
3497
1452
                args[1] = banner_module.ConfigFile.format_variable(args[1], bool, args[0])
1453
                args[1] = banner_module.ConfigFile.format_variable(args[1], bool, args[0])
3498
1453
            except Exception:
1454
            except Exception:
3499
@@ -2470,7 +2471,7 @@
3500
2470
2471
3501
2471
        if len(args) == 2:
2472
        if len(args) == 2:
3502
2472
            if args[1] in ['group_subprocesses', 'complex_mass_scheme',\
2473
            if args[1] in ['group_subprocesses', 'complex_mass_scheme',\
3504
2473
                           'loop_optimized_output', 'loop_color_flows',\
2474
                           'loop_optimized_output', 'remove_os', 'loop_color_flows',\
3505
2474
                           'low_mem_multicore_nlo_generation']:
2475
                           'low_mem_multicore_nlo_generation']:
3506
2475
                return self.list_completion(text, ['False', 'True', 'default'])
2476
                return self.list_completion(text, ['False', 'True', 'default'])
3507
2476
            elif args[1] in ['ignore_six_quark_processes']:
2477
            elif args[1] in ['ignore_six_quark_processes']:
3508
@@ -2863,6 +2864,7 @@
3509
2863
                          'gauge':'unitary',
2864
                          'gauge':'unitary',
3510
2864
                          'stdout_level':None,
2865
                          'stdout_level':None,
3511
2865
                          'loop_optimized_output':True,
2866
                          'loop_optimized_output':True,
3512
2867
                          'remove_os':False,
3513
2866
                          'loop_color_flows':False,
2868
                          'loop_color_flows':False,
3514
2867
                          'max_npoint_for_channel': 0, # 0 means automaticly adapted
2869
                          'max_npoint_for_channel': 0, # 0 means automaticly adapted
3515
2868
                          'default_unset_couplings': 99 # 99 means infinity
2870
                          'default_unset_couplings': 99 # 99 means infinity
3516
2869
2871
3517
=== modified file 'madgraph/iolibs/export_fks.py'
3518
--- madgraph/iolibs/export_fks.py	2017-12-06 15:01:35 +0000
3519
+++ madgraph/iolibs/export_fks.py	2018-09-06 19:46:41 +0000
3520
@@ -611,6 +611,7 @@
3521
611
                     'fks_Sij.f',
611
                     'fks_Sij.f',
3522
612
                     'fks_powers.inc',
612
                     'fks_powers.inc',
3523
613
                     'fks_singular.f',
613
                     'fks_singular.f',
3524
614
                     'transform_os.f',
3525
614
                     'veto_xsec.f',
615
                     'veto_xsec.f',
3526
615
                     'veto_xsec.inc',
616
                     'veto_xsec.inc',
3527
616
                     'weight_lines.f',
617
                     'weight_lines.f',
3528
@@ -679,6 +680,17 @@
3529
679
        files.append_to_file(filename,
680
        files.append_to_file(filename,
3530
680
                             self.write_subproc,
681
                             self.write_subproc,
3531
681
                             borndir)
682
                             borndir)
3532
683
<<<<<<< TREE
3533
684
=======
3534
685
3535
686
        # Add the os_ids to os_ids.mg
3536
687
        os_ids = matrix_element.get_os_ids()
3537
688
        filename = 'os_ids.mg'
3538
689
        files.append_to_file(filename,
3539
690
                             self.write_os_ids,
3540
691
                             borndir,
3541
692
                             os_ids)
3542
693
>>>>>>> MERGE-SOURCE
3543
682
            
694
            
3544
683
        os.chdir(cwd)
695
        os.chdir(cwd)
3545
684
        # Generate info page
696
        # Generate info page
3546
@@ -927,6 +939,7 @@
3547
927
                self.opt['madanalysis5_path'], pjoin(self.dir_path,'Cards'),
939
                self.opt['madanalysis5_path'], pjoin(self.dir_path,'Cards'),
3548
928
                levels =['hadron'])
940
                levels =['hadron'])
3549
929
941
3550
942
3551
930
    def write_real_from_born_configs(self, writer, matrix_element, fortran_model):
943
    def write_real_from_born_configs(self, writer, matrix_element, fortran_model):
3552
931
        """Writes the real_from_born_configs.inc file that contains
944
        """Writes the real_from_born_configs.inc file that contains
3553
932
        the mapping to go for a given born configuration (that is used
945
        the mapping to go for a given born configuration (that is used
3554
@@ -1403,15 +1416,22 @@
3555
1403
implicit none
1416
implicit none
3556
1404
include 'nexternal.inc'
1417
include 'nexternal.inc'
3557
1405
double precision p(0:3, nexternal)
1418
double precision p(0:3, nexternal)
3559
1406
double precision wgt
1419
double precision wgt, wgt_os
3560
1407
integer nfksprocess
1420
integer nfksprocess
3561
1408
common/c_nfksprocess/nfksprocess
1421
common/c_nfksprocess/nfksprocess
3562
1409
"""
1422
"""
3563
1423
        # subtract here the on shell matrix-elements if any is there
3564
1410
        for n, info in enumerate(matrix_element.get_fks_info_list()):
1424
        for n, info in enumerate(matrix_element.get_fks_info_list()):
3565
1425
            os_lines = ''
3566
1426
            for i, os_me in \
3567
1427
              enumerate(matrix_element.real_processes[info['n_me'] - 1].os_matrix_elements):
3568
1428
                os_lines += '\n call smatrix_%d_os_%d_wrapper(p, wgt_os)\n wgt = wgt - wgt_os' \
3569
1429
                        % (info['n_me'] , i + 1)
3570
1430
3571
1411
            file += \
1431
            file += \
3572
1412
"""if (nfksprocess.eq.%(n)d) then
1432
"""if (nfksprocess.eq.%(n)d) then
3575
1413
call smatrix_%(n_me)d(p, wgt)
1433
call smatrix_%(n_me)d(p, wgt) %(os_lines)s
3576
1414
else""" % {'n': n + 1, 'n_me' : info['n_me']}
1434
else""" % {'n': n + 1, 'n_me' : info['n_me'], 'os_lines': os_lines}
3577
1415
1435
3578
1416
        if matrix_element.real_processes:
1436
        if matrix_element.real_processes:
3579
1417
            file += \
1437
            file += \
3580
@@ -1439,11 +1459,11 @@
3581
1439
        as well as for all the real emission processes"""
1459
        as well as for all the real emission processes"""
3582
1440
1460
3583
1441
        filename = 'born.ps'
1461
        filename = 'born.ps'
3584
1462
        model = matrix_element.born_matrix_element.get('processes')[0].get('model')
3585
1442
        plot = draw.MultiEpsDiagramDrawer(matrix_element.born_matrix_element.\
1463
        plot = draw.MultiEpsDiagramDrawer(matrix_element.born_matrix_element.\
3586
1443
                                    get('base_amplitude').get('diagrams'),
1464
                                    get('base_amplitude').get('diagrams'),
3587
1444
                                    filename,
1465
                                    filename,
3590
1445
                                    model=matrix_element.born_matrix_element.\
1466
                                    model=model,
3589
1446
                                    get('processes')[0].get('model'),
3591
1447
                                    amplitude=True, diagram_type='born')
1467
                                    amplitude=True, diagram_type='born')
3592
1448
        plot.draw()
1468
        plot.draw()
3593
1449
1469
3594
@@ -1452,22 +1472,125 @@
3595
1452
            plot = draw.MultiEpsDiagramDrawer(fksreal.matrix_element.\
1472
            plot = draw.MultiEpsDiagramDrawer(fksreal.matrix_element.\
3596
1453
                                        get('base_amplitude').get('diagrams'),
1473
                                        get('base_amplitude').get('diagrams'),
3597
1454
                                        filename,
1474
                                        filename,
3600
1455
                                        model=fksreal.matrix_element.\
1475
                                        model=model,
3599
1456
                                        get('processes')[0].get('model'),
3601
1457
                                        amplitude=True, diagram_type='real')
1476
                                        amplitude=True, diagram_type='real')
3602
1458
            plot.draw()
1477
            plot.draw()
3603
1478
            for nos, os_me in enumerate(fksreal.os_matrix_elements):
3604
1479
                suffix = '%d_os_%d' % (n + 1, nos + 1)
3605
1480
                filename = 'matrix_%s.ps' % suffix
3606
1481
                plot = draw.MultiEpsDiagramDrawer(os_me.\
3607
1482
                                        get('base_amplitude').get('diagrams'),
3608
1483
                                        filename,
3609
1484
                                        model=model,
3610
1485
                                        amplitude=True, diagram_type='real')
3611
1486
                plot.draw()
3612
1459
1487
3613
1460
1488
3614
1461
    def write_real_matrix_elements(self, matrix_element, fortran_model):
1489
    def write_real_matrix_elements(self, matrix_element, fortran_model):
3616
1462
        """writes the matrix_i.f files which contain the real matrix elements""" 
1490
        """writes the matrix_i.f files which contain the real matrix elements
3617
1491
        and the matrix_i_os_j.f which contain eventual on shell subtraction
3618
1492
        terms""" 
3619
1463
1493
3620
1464
1494
3621
1465
1495
3622
1466
        for n, fksreal in enumerate(matrix_element.real_processes):
1496
        for n, fksreal in enumerate(matrix_element.real_processes):
3623
1467
            filename = 'matrix_%d.f' % (n + 1)
1497
            filename = 'matrix_%d.f' % (n + 1)
3624
1468
            self.write_matrix_element_fks(writers.FortranWriter(filename),
1498
            self.write_matrix_element_fks(writers.FortranWriter(filename),
3627
1469
                                            fksreal.matrix_element, n + 1, 
1499
                                          fksreal.matrix_element, n + 1, 
3628
1470
                                            fortran_model)
1500
                                          fortran_model, 
3629
1501
                                          os_info = {'diags': fksreal.os_diagrams, 'ids': fksreal.os_ids})
3630
1502
3631
1503
            for nos, os_me in enumerate(fksreal.os_matrix_elements):
3632
1504
                suffix = '%d_os_%d' % (n + 1, nos + 1)
3633
1505
                filename = 'matrix_%s.f' % suffix
3634
1506
                self.write_matrix_element_fks(writers.FortranWriter(filename),
3635
1507
                                            os_me, suffix, fortran_model, 
3636
1508
                                            os_info = {'diags': [], 'ids': fksreal.os_ids})
3637
1509
                filename = 'wrapper_matrix_%s.f' % suffix
3638
1510
                self.write_os_wrapper(writers.FortranWriter(filename),
3639
1511
                        fksreal.matrix_element, os_me, suffix, fortran_model)
3640
1512
3641
1513
3642
1514
    def write_os_wrapper(self, writer, real_me, os_me, suffix, fortran_model):
3643
1515
        """write the wrapper for the on shell subtraction matrix-elements
3644
1516
        which takes care of reordering the momenta and of knowing which is the 
3645
1517
        mother particle"""
3646
1518
        replace_dict = {}
3647
1519
        replace_dict['suffix'] = suffix
3648
1520
3649
1521
        # find the permutation of the final state legs to map real_me onto os_me. 
3650
1522
        # look only at final state legs (initial state legs are not touched)
3651
1523
        real_ids = [l['id'] for l in real_me.get_base_amplitude()['process']['legs'] if l['state']]
3652
1524
        os_ids = os_me.get_base_amplitude()['process'].get_final_ids_after_decay()
3653
1525
        nexternal,ninitial = real_me.get_nexternal_ninitial()
3654
1526
        permutation = []
3655
1527
        #initial state legs are trivial
3656
1528
        for i in range(ninitial):
3657
1529
            permutation.append(i)
3658
1530
        for os_id in os_ids:
3659
1531
            permutation.append(ninitial + real_ids.index(os_id))
3660
1532
            # don't remove from the list, otherwise the position is
3661
1533
            # not correcly returned, just replace it by an 'x' 
3662
1534
            real_ids[real_ids.index(os_id)] = 'x'
3663
1535
        replace_dict['mom_perm'] = ', '.join([str(pp + 1) for pp in permutation])
3664
1536
        # find decay mother and daughter id's
3665
1537
        mother = [l['id'] \
3666
1538
                for l in os_me.get_base_amplitude()['process']['decay_chains'][0]['legs'] \
3667
1539
                if not l['state']]
3668
1540
        daughters = [l['id'] \
3669
1541
                for l in os_me.get_base_amplitude()['process']['decay_chains'][0]['legs'] \
3670
1542
                if l['state']]
3671
1543
        if not (len(mother) == 1 and len(daughters) == 2):
3672
1544
            raise fks_common.FKSProcessError(
3673
1545
                    'Incorrect number of mother(s) and daughters: %d, %d' % \
3674
1546
                            (len(mother), len(daughters)))
3675
1547
        
3676
1548
        model = os_me.get_base_amplitude()['process']['model']
3677
1549
3678
1550
        replace_dict['mom_external'] = {True: '.true.', False: '.false.'}[\
3679
1551
                    mother[0] in os_ids or \
3680
1552
                    model.get_particle(mother[0]).get_anti_pdg_code() in os_ids]
3681
1553
3682
1554
        # mother and daughter masses and widths
3683
1555
        replace_dict['mom_mass'] = model.get_particle(mother[0])['mass']
3684
1556
        replace_dict['mom_wdth'] = model.get_particle(mother[0])['width']
3685
1557
        replace_dict['dau1_mass'] = model.get_particle(daughters[0])['mass']
3686
1558
        replace_dict['dau2_mass'] = model.get_particle(daughters[1])['mass']
3687
1559
        # position of daughter in the array of momenta (the one of the decayed process)
3688
1560
        # count the ovvurrence of the daughters into the final state:
3689
1561
        counts = [0, 0]
3690
1562
        for idau, dau in enumerate(daughters):
3691
1563
            for idd in os_ids:
3692
1564
                if idd == dau:
3693
1565
                    counts[idau] +=1
3694
1566
3695
1567
        if counts == [1,1]:
3696
1568
            # if daughters are unique, find them in the os_ids list
3697
1569
            replace_dict['idau1'] = os_ids.index(daughters[0]) + ninitial + 1
3698
1570
            replace_dict['idau2'] = os_ids.index(daughters[1]) + ninitial + 1
3699
1571
        else:
3700
1572
            # otherwise, assign the position of the mother and the next one
3701
1573
            real_ids = [l['id'] for l in os_me.get_base_amplitude()['process']['legs'] if l['state']]
3702
1574
            replace_dict['idau1'] = real_ids.index(mother[0]) + ninitial + 1
3703
1575
            replace_dict['idau2'] = real_ids.index(mother[0]) + ninitial + 2
3704
1576
3705
1577
        # find the spectator (needed by the function which put momenta on-shell)
3706
1578
        # by default choose the first final state particle which is not a daughter
3707
1579
        for i in range(ninitial,nexternal):
3708
1580
            if i + 1 not in [replace_dict['idau1'], replace_dict['idau2']]:
3709
1581
                spectator = i + 1
3710
1582
                break
3711
1583
        replace_dict['ispect'] = spectator 
3712
1584
        replace_dict['spect_mass'] = model.get_particle(spectator)['mass']
3713
1585
3714
1586
        # finally write out the file
3715
1587
        file = open(os.path.join(_file_path, \
3716
1588
                          'iolibs/template_files/os_wrapper_fks.inc')).read()
3717
1589
        file = file % replace_dict
3718
1590
        
3719
1591
        # Write the file
3720
1592
        writer.writelines(file)
3721
1593
3722
1471
1594
3723
1472
    def write_pdf_calls(self, matrix_element, fortran_model):
1595
    def write_pdf_calls(self, matrix_element, fortran_model):
3724
1473
        """writes the parton_lum_i.f files which contain the real matrix elements.
1596
        """writes the parton_lum_i.f files which contain the real matrix elements.
3725
@@ -2425,7 +2548,7 @@
3726
2425
    # write_matrix_element_fks
2548
    # write_matrix_element_fks
3727
2426
    #===============================================================================
2549
    #===============================================================================
3728
2427
    #test written
2550
    #test written
3730
2428
    def write_matrix_element_fks(self, writer, matrix_element, n, fortran_model):
2551
    def write_matrix_element_fks(self, writer, matrix_element, n, fortran_model, os_info={}):
3731
2429
        """Export a matrix element to a matrix.f file in MG4 madevent format"""
2552
        """Export a matrix element to a matrix.f file in MG4 madevent format"""
3732
2430
    
2553
    
3733
2431
        if not matrix_element.get('processes') or \
2554
        if not matrix_element.get('processes') or \
3734
@@ -2439,7 +2562,7 @@
3735
2439
        writers.FortranWriter.downcase = False
2562
        writers.FortranWriter.downcase = False
3736
2440
    
2563
    
3737
2441
        replace_dict = {}
2564
        replace_dict = {}
3739
2442
        replace_dict['N_me'] = n
2565
        replace_dict['N_me'] = str(n)
3740
2443
    
2566
    
3741
2444
        # Extract version number and date from VERSION file
2567
        # Extract version number and date from VERSION file
3742
2445
        info_lines = self.get_mg5_info_lines()
2568
        info_lines = self.get_mg5_info_lines()
3743
@@ -2471,6 +2594,8 @@
3744
2471
        replace_dict['ngraphs'] = ngraphs
2594
        replace_dict['ngraphs'] = ngraphs
3745
2472
    
2595
    
3746
2473
        # Extract ncolor
2596
        # Extract ncolor
3747
2597
        if not matrix_element.get('color_basis'):
3748
2598
            matrix_element.process_color()
3749
2474
        ncolor = max(1, len(matrix_element.get('color_basis')))
2599
        ncolor = max(1, len(matrix_element.get('color_basis')))
3750
2475
        replace_dict['ncolor'] = ncolor
2600
        replace_dict['ncolor'] = ncolor
3751
2476
    
2601
    
3752
@@ -2482,6 +2607,24 @@
3753
2482
        helas_calls = fortran_model.get_matrix_element_calls(\
2607
        helas_calls = fortran_model.get_matrix_element_calls(\
3754
2483
                    matrix_element)
2608
                    matrix_element)
3755
2484
        replace_dict['helas_calls'] = "\n".join(helas_calls)
2609
        replace_dict['helas_calls'] = "\n".join(helas_calls)
3756
2610
3757
2611
        # if there are os_diagrams, these should be set to zero in the ME if
3758
2612
        # diagram removal without interference is done
3759
2613
        if os_info and os_info['diags']:
3760
2614
            # this is for the resonant diagrams in the full real emission ME
3761
2615
            os_diagrams = os_info['diags']
3762
2616
            os_ids = os_info['ids']
3763
2617
            replace_dict['helas_calls'] = \
3764
2618
                self.change_width_in_os_diagrams(replace_dict['helas_calls'], os_diagrams, os_ids)
3765
2619
            replace_dict['helas_calls'] += '\n' + \
3766
2620
                    self.get_os_diagrams_lines(matrix_element, os_diagrams, os_ids)
3767
2621
3768
2622
        elif os_info and os_info['ids'] and not os_info['diags']:
3769
2623
            os_ids = os_info['ids']
3770
2624
            # this is for the OS subtraction counterterms, 
3771
2625
            # in this case replace all occurrences of the particle width
3772
2626
            replace_dict['helas_calls'] = \
3773
2627
                self.change_width_in_os_diagrams(replace_dict['helas_calls'], [], os_ids)
3774
2485
    
2628
    
3775
2486
        # Extract nwavefuncs (important to place after get_matrix_element_calls
2629
        # Extract nwavefuncs (important to place after get_matrix_element_calls
3776
2487
        # so that 'me_id' is set)
2630
        # so that 'me_id' is set)
3777
@@ -2514,6 +2657,49 @@
3778
2514
        return len(filter(lambda call: call.find('#') != 0, helas_calls)), ncolor
2657
        return len(filter(lambda call: call.find('#') != 0, helas_calls)), ncolor
3779
2515
2658
3780
2516
2659
3781
2660
    def get_os_diagrams_lines(self, matrix_element, os_diagrams, os_ids):
3782
2661
        """ add the lines which set to zero the diagrams used with
3783
2662
        diagram-removal techniques
3784
2663
        """
3785
2664
        particle_dict = self.model.get('particle_dict') 
3786
2665
3787
2666
        text = 'if (iossubtr.eq.1) then\n'
3788
2667
        for diags, ids in zip(os_diagrams, os_ids): 
3789
2668
            text += 'if (%s.gt.(%s+%s)) then\n' % \
3790
2669
                    tuple([particle_dict[idd].get('mass') for idd in ids])
3791
2670
            for diag in diags:
3792
2671
                for amp in matrix_element['diagrams'][diag]['amplitudes']:
3793
2672
                    text+= 'amp(%d) = dcmplx(0d0,0d0)\n' % amp['number']
3794
2673
            text += 'endif\n'
3795
2674
        text += 'endif\n'
3796
2675
3797
2676
        return text
3798
2677
3799
2678
3800
2679
    def change_width_in_os_diagrams(self, helas_calls, os_diagrams, os_ids):
3801
2680
        """change the name of the width used in diagrams with internal resonances, so
3802
2681
        that the width in those diagrams is not set to zero
3803
2682
        """
3804
2683
        diagrams_text = helas_calls.split('# Amplitude')
3805
2684
        if os_diagrams:
3806
2685
            for diags, ids in zip(os_diagrams, os_ids):
3807
2686
                part_width = self.model.get('particle_dict')[ids[0]].get('width')
3808
2687
3809
2688
                for diag in diags:
3810
2689
                    if part_width + '_keep' not in diagrams_text:
3811
2690
                        diagrams_text[diag] = diagrams_text[diag].replace(part_width, part_width + '_keep')
3812
2691
            return '# Amplitude'.join(diagrams_text)
3813
2692
        else:
3814
2693
            new_helas_calls = copy.copy(helas_calls)
3815
2694
            for ids in os_ids:
3816
2695
                part_width = self.model.get('particle_dict')[ids[0]].get('width')
3817
2696
                if part_width + '_keep' not in new_helas_calls:
3818
2697
                    new_helas_calls = new_helas_calls.replace(part_width, part_width + '_keep')
3819
2698
            return new_helas_calls
3820
2699
3821
2700
3822
2701
3823
2702
3824
2517
    #===============================================================================
2703
    #===============================================================================
3825
2518
    # write_pdf_file
2704
    # write_pdf_file
3826
2519
    #===============================================================================
2705
    #===============================================================================
3827
@@ -3353,6 +3539,21 @@
3828
3353
        return True
3539
        return True
3829
3354
3540
3830
3355
3541
3831
3542
    #===========================================================================
3832
3543
    # write_os_ids
3833
3544
    #===========================================================================
3834
3545
    def write_os_ids(self, writer, folder, os_ids):
3835
3546
        """Append the os_ids to the os_ids file"""
3836
3547
3837
3548
        # Write line to file
3838
3549
        content = ''
3839
3550
        if os_ids:
3840
3551
            content+= '%s: %s\n' % (folder, ' '.join(['%d' % v for v in os_ids])) 
3841
3552
        writer.write(content)
3842
3553
3843
3554
        return True
3844
3555
3845
3556
3846
3356
3557
3847
3357
3558
3848
3358
3559
3849
3359
3560
3850
=== modified file 'madgraph/iolibs/export_v4.py'
3851
--- madgraph/iolibs/export_v4.py	2018-06-07 13:31:42 +0000
3852
+++ madgraph/iolibs/export_v4.py	2018-09-06 19:46:41 +0000
3853
@@ -5733,6 +5733,14 @@
3854
5733
        if widths:
5733
        if widths:
3855
5734
            fsock.writelines('double precision '+','.join(widths)+'\n')
5734
            fsock.writelines('double precision '+','.join(widths)+'\n')
3856
5735
            fsock.writelines('common/widths/ '+','.join(widths)+'\n\n')
5735
            fsock.writelines('common/widths/ '+','.join(widths)+'\n\n')
3857
5736
            if 'remove_os' in self.opt.keys() and self.opt['remove_os']:
3858
5737
                # add the block of physical widths, which should contain the 
3859
5738
                # values of the width before they are forced to be zero
3860
5739
                # in amcatnlo
3861
5740
                phys_widths = [w + '_keep' for w in widths]
3862
5741
                fsock.writelines('double precision '+','.join(phys_widths)+'\n')
3863
5742
                fsock.writelines('common/keep_widths/ '+','.join(phys_widths)+'\n\n')
3864
5743
3865
5736
            if self.opt['mp']:
5744
            if self.opt['mp']:
3866
5737
                mp_fsock_same_name.writelines(self.mp_real_format+' '+\
5745
                mp_fsock_same_name.writelines(self.mp_real_format+' '+\
3867
5738
                                                          ','.join(widths)+'\n')
5746
                                                          ','.join(widths)+'\n')
3868
5739
5747
3869
=== added file 'madgraph/iolibs/template_files/os_wrapper_fks.inc'
3870
--- madgraph/iolibs/template_files/os_wrapper_fks.inc	1970-01-01 00:00:00 +0000
3871
+++ madgraph/iolibs/template_files/os_wrapper_fks.inc	2018-09-06 19:46:41 +0000
3872
@@ -0,0 +1,137 @@
3873
1
subroutine smatrix_%(suffix)s_wrapper(p, wgt_os)
3874
2
C this function computes the "on-shell subtraction term", which is subtracted 
3875
3
C by the full real-emission matrix-element when some on-shell aprticles appear.
3876
4
implicit none
3877
5
include 'nexternal.inc'
3878
6
C Arguments
3879
7
double precision p(0:3, nexternal)
3880
8
double precision wgt_os
3881
9
C internal variables
3882
10
double precision p_os(0:3, nexternal)
3883
11
double precision p_reord(0:3, nexternal)
3884
12
double precision wgt
3885
13
double precision mom_mass, mom_wdth, dau1_mass, dau2_mass, spect_mass
3886
14
C for the OS subtraction
3887
15
integer iossubtr
3888
16
common /to_os_reshuf/iossubtr
3889
17
integer idau1, idau2, ispect
3890
18
parameter (idau1 = %(idau1)d)
3891
19
parameter (idau2 = %(idau2)d)
3892
20
parameter (ispect = %(ispect)d)
3893
21
integer mom_perm(nexternal), mom_perm_pass(nexternal), ii, i, j
3894
22
data mom_perm / %(mom_perm)s /
3895
23
include 'coupl.inc'
3896
24
double precision ZERO
3897
25
parameter (ZERO = 0d0)
3898
26
logical firsttime
3899
27
data firsttime /.true./
3900
28
double precision pdfratio, bwratio, fluxratio
3901
29
integer ibw
3902
30
integer stat
3903
31
integer os_countall, os_countbad
3904
32
common /to_os_count/ os_countall, os_countbad
3905
33
C Functions 
3906
34
double precision dot
3907
35
3908
36
os_countall = os_countall + 1
3909
37
3910
38
C do nothing for diagram removal without interference (iossubtr==1) or when no subtraction is performed (iossubtr==0), amplitudes are set to zero directly in the matrix_X.f
3911
39
wgt_os = 0d0
3912
40
if (iossubtr.eq.1.or.iossubtr.eq.0) return
3913
41
3914
42
mom_wdth = %(mom_wdth)s_keep
3915
43
3916
44
mom_mass = %(mom_mass)s
3917
45
dau1_mass = %(dau1_mass)s
3918
46
dau2_mass = %(dau2_mass)s
3919
47
spect_mass = %(spect_mass)s
3920
48
3921
49
C if daughters are heavier than the mother, nothing has to be done
3922
50
if ((dau1_mass+dau2_mass).gt.(mom_mass)) return
3923
51
3924
52
3925
53
C first reorder the momenta
3926
54
do j = 1, nexternal
3927
55
do i = 0, 3
3928
56
p_reord(i, j) = p(i, mom_perm(j))
3929
57
enddo
3930
58
enddo
3931
59
3932
60
C consistency check to assure the reshuffling was consistent
3933
61
if (firsttime) then
3934
62
if (dabs(dsqrt(dot(p_reord(0, idau1), p_reord(0, idau1))) - dau1_mass) / mom_mass .gt. 1d-3) then
3935
63
write(*,*) 'DAUGHTER 1 NOT ON SHELL', dsqrt(dot(p_reord(0, idau1), p_reord(0, idau1))), dau1_mass
3936
64
stop
3937
65
endif
3938
66
3939
67
if (dabs(dsqrt(dot(p_reord(0, idau2), p_reord(0, idau2))) - dau2_mass) / mom_mass .gt. 1d-3) then
3940
68
write(*,*) 'DAUGHTER 2 NOT ON SHELL', dsqrt(dot(p_reord(0, idau2), p_reord(0, idau2))), dau1_mass
3941
69
stop
3942
70
endif
3943
71
3944
72
firsttime = .false.
3945
73
endif
3946
74
3947
75
stat=0 !(not all the reshuffling techniques provide stat, so set it to 0 from scratch)
3948
76
3949
77
C now do the reshuffling of the momenta according to the different values of iossubtr
3950
78
C 0 and 1 have already been treated above
3951
79
if (iossubtr.eq.2) then
3952
80
C  iossubtr==2 -> DR with interferece
3953
81
  call transform_os_ident(p_reord, p_os)
3954
82
elseif (iossubtr.eq.3.or.iossubtr.eq.4) then
3955
83
C  iossubtr = 3 -> DS with reshuffling on initial state, standard BW
3956
84
C  iossubtr = 4 -> DS with reshuffling on initial state, running BW
3957
85
  call transform_os_init(p_reord, p_os, idau1, idau2, dau1_mass, dau2_mass, mom_mass)
3958
86
else if (iossubtr.eq.5.or.iossubtr.eq.6) then
3959
87
C  iossubtr = 5 -> DS with reshuffling on all FS particles, standard BW
3960
88
C  iossubtr = 6 -> DS with reshuffling on all FS particles, running BW
3961
89
call transform_os_final(p_reord, p_os, idau1, idau2, dau1_mass, dau2_mass, mom_mass, stat)
3962
90
else if (iossubtr.eq.7.or.iossubtr.eq.8) then
3963
91
C  iossubtr = 7 -> DS with option C, standard BW
3964
92
C  iossubtr = 8 -> DS with option C, running BW
3965
93
C copy the permutation and pass it to the reshuffling subroutine
3966
94
do i = 1, nexternal
3967
95
mom_perm_pass(i) = mom_perm(i)
3968
96
enddo
3969
97
call transform_os_genps(p_reord, p_os, idau1, idau2, dau1_mass, dau2_mass, mom_mass, mom_perm_pass, stat)
3970
98
else if (iossubtr.eq.9.or.iossubtr.eq.10) then
3971
99
C  iossubtr = 9 -> DS with reshuffling on spectator, standard BW
3972
100
C  iossubtr = 10-> DS with reshuffling on spectator, running BW
3973
101
  call transform_os_spect(p_reord, p_os, idau1, idau2, ispect, dau1_mass, dau2_mass, spect_mass, mom_mass, stat)
3974
102
else
3975
103
  write(*,*) 'ERROR, iossubtr not implemented', iossubtr
3976
104
  stop 1
3977
105
endif
3978
106
3979
107
C if stat != 0, the reshuffling was not possible. just exit
3980
108
if (stat.ne.0) then
3981
109
  os_countbad = os_countbad + 1
3982
110
  return
3983
111
endif
3984
112
3985
113
C now set iBW for the BW ratio
3986
114
C ibw==0, no BW
3987
115
C ibw==1, standard BW
3988
116
C ibw==2, running BW
3989
117
3990
118
if (iossubtr.eq.2) then
3991
119
  ibw=0
3992
120
else if (iossubtr.eq.3.or.iossubtr.eq.5.or.iossubtr.eq.7.or.iossubtr.eq.9) then
3993
121
  ibw=1
3994
122
else if (iossubtr.eq.4.or.iossubtr.eq.6.or.iossubtr.eq.8.or.iossubtr.eq.10) then
3995
123
  ibw=2
3996
124
endif
3997
125
3998
126
C compute BW and flux/PDF ratio
3999
127
call get_bw_ratio(p_reord, mom_mass, mom_wdth, idau1, idau2, ibw, bwratio) 
4000
128
call get_pdf_flux_ratio(p_reord, p_os, pdfratio, fluxratio)
4001
129
4002
130
C finally call the resonant matrix element with the reshuffled momenta
4003
131
call smatrix_%(suffix)s(p_os, wgt)
4004
132
C and add the weigth, compensating for the reshuffling
4005
133
wgt_os = wgt * pdfratio * bwratio * fluxratio
4006
134
4007
135
999 continue
4008
136
return
4009
137
end
4010
0
138
4011
=== modified file 'madgraph/iolibs/template_files/realmatrix_fks.inc'
4012
--- madgraph/iolibs/template_files/realmatrix_fks.inc	2016-09-08 23:15:34 +0000
4013
+++ madgraph/iolibs/template_files/realmatrix_fks.inc	2018-09-06 19:46:41 +0000
4014
@@ -1,4 +1,4 @@
4016
1
      SUBROUTINE SMATRIX_%(N_me)d(P,ANS)
1
      SUBROUTINE SMATRIX_%(N_me)s(P,ANS)
4017
2
C  
2
C  
4018
3
%(info_lines)s
3
%(info_lines)s
4019
4
C
4
C
4020
@@ -25,7 +25,7 @@
4021
25
C LOCAL VARIABLES 
25
C LOCAL VARIABLES 
4022
26
C  
26
C  
4023
27
      INTEGER IHEL,IDEN,I,T_IDENT(NCOMB)
27
      INTEGER IHEL,IDEN,I,T_IDENT(NCOMB)
4025
28
      REAL*8 MATRIX_%(N_me)d
28
      REAL*8 MATRIX_%(N_me)s
4026
29
      REAL*8 T,T_SAVE(NCOMB)
29
      REAL*8 T,T_SAVE(NCOMB)
4027
30
      SAVE T_SAVE,T_IDENT
30
      SAVE T_SAVE,T_IDENT
4028
31
      INTEGER NHEL(NEXTERNAL,NCOMB)
31
      INTEGER NHEL(NEXTERNAL,NCOMB)
4029
@@ -45,7 +45,7 @@
4030
45
          IF (NTRY.LT.2) THEN
45
          IF (NTRY.LT.2) THEN
4031
46
C for the first ps-point, check for helicities that give
46
C for the first ps-point, check for helicities that give
4032
47
C identical matrix elements
47
C identical matrix elements
4034
48
            T=MATRIX_%(N_me)d(P ,NHEL(1,IHEL))            
48
            T=MATRIX_%(N_me)s(P ,NHEL(1,IHEL))            
4035
49
            T_SAVE(IHEL)=T
49
            T_SAVE(IHEL)=T
4036
50
            T_IDENT(IHEL)=-1
50
            T_IDENT(IHEL)=-1
4037
51
            DO I=1,IHEL-1
51
            DO I=1,IHEL-1
4038
@@ -62,7 +62,7 @@
4039
62
              T=T_SAVE(T_IDENT(IHEL))
62
              T=T_SAVE(T_IDENT(IHEL))
4040
63
              T_SAVE(IHEL)=T
63
              T_SAVE(IHEL)=T
4041
64
            ELSE
64
            ELSE
4043
65
              T=MATRIX_%(N_me)d(P ,NHEL(1,IHEL))            
65
              T=MATRIX_%(N_me)s(P ,NHEL(1,IHEL))            
4044
66
              T_SAVE(IHEL)=T
66
              T_SAVE(IHEL)=T
4045
67
            ENDIF
67
            ENDIF
4046
68
          ENDIF
68
          ENDIF
4047
@@ -78,7 +78,7 @@
4048
78
      END
78
      END
4049
79
       
79
       
4050
80
       
80
       
4052
81
      REAL*8 FUNCTION MATRIX_%(N_me)d(P,NHEL)
81
      REAL*8 FUNCTION MATRIX_%(N_me)s(P,NHEL)
4053
82
C  
82
C  
4054
83
%(info_lines)s
83
%(info_lines)s
4055
84
C
84
C
4056
@@ -114,6 +114,10 @@
4057
114
      data ic /nexternal*1/
114
      data ic /nexternal*1/
4058
115
      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
115
      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
4059
116
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP(NCOLOR), W(%(wavefunctionsize)d,NWAVEFUNCS)
116
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP(NCOLOR), W(%(wavefunctionsize)d,NWAVEFUNCS)
4060
117
4061
118
C for the OS subtraction
4062
119
      integer iossubtr
4063
120
      common /to_os_reshuf/iossubtr
4064
117
C  
121
C  
4065
118
C COLOR DATA
122
C COLOR DATA
4066
119
C  
123
C  
4067
@@ -123,12 +127,12 @@
4068
123
C ----------
127
C ----------
4069
124
%(helas_calls)s
128
%(helas_calls)s
4070
125
%(jamp_lines)s
129
%(jamp_lines)s
4072
126
      MATRIX_%(N_me)d = 0.D0 
130
      MATRIX_%(N_me)s = 0.D0 
4073
127
      DO I = 1, NCOLOR
131
      DO I = 1, NCOLOR
4074
128
          ZTEMP = (0.D0,0.D0)
132
          ZTEMP = (0.D0,0.D0)
4075
129
          DO J = 1, NCOLOR
133
          DO J = 1, NCOLOR
4076
130
              ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
134
              ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
4077
131
          ENDDO
135
          ENDDO
4079
132
          MATRIX_%(N_me)d = MATRIX_%(N_me)d+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
136
          MATRIX_%(N_me)s = MATRIX_%(N_me)s+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
4080
133
      ENDDO
137
      ENDDO
4081
134
      END
138
      END
4082
135
139
4083
=== modified file 'madgraph/various/banner.py'
4084
--- madgraph/various/banner.py	2018-07-22 11:40:04 +0000
4085
+++ madgraph/various/banner.py	2018-09-06 19:46:41 +0000
4086
@@ -1426,6 +1426,7 @@
4087
1426
        self.add_param('ninitial', 0)
1426
        self.add_param('ninitial', 0)
4088
1427
        self.add_param('grouped_matrix', True)
1427
        self.add_param('grouped_matrix', True)
4089
1428
        self.add_param('has_loops', False)
1428
        self.add_param('has_loops', False)
4090
1429
        self.add_param('has_os_div', False)
4091
1429
        self.add_param('bias_module','None')
1430
        self.add_param('bias_module','None')
4092
1430
        self.add_param('max_n_matched_jets', 0)
1431
        self.add_param('max_n_matched_jets', 0)
4093
1431
        self.add_param('colored_pdgs', [1,2,3,4,5])
1432
        self.add_param('colored_pdgs', [1,2,3,4,5])
4094
@@ -3714,6 +3715,8 @@
4095
3714
        self.add_param('pdf_set_min', 244601, hidden=True)
3715
        self.add_param('pdf_set_min', 244601, hidden=True)
4096
3715
        self.add_param('pdf_set_max', 244700, hidden=True)
3716
        self.add_param('pdf_set_max', 244700, hidden=True)
4097
3716
        self.add_param('store_rwgt_info', False)
3717
        self.add_param('store_rwgt_info', False)
4098
3718
        #removal of resonances from real emissions
4099
3719
        self.add_param('iossubtr', 1)
4100
3717
        self.add_param('systematics_program', 'none', include=False, hidden=True, comment='Choose which program to use for systematics computation: none, systematics')
3720
        self.add_param('systematics_program', 'none', include=False, hidden=True, comment='Choose which program to use for systematics computation: none, systematics')
4101
3718
        self.add_param('systematics_arguments', [''], include=False, hidden=True, comment='Choose the argment to pass to the systematics command. like --mur=0.25,1,4. Look at the help of the systematics function for more details.')
3721
        self.add_param('systematics_arguments', [''], include=False, hidden=True, comment='Choose the argment to pass to the systematics command. like --mur=0.25,1,4. Look at the help of the systematics function for more details.')
4102
3719
             
3722
             
4103
3720
3723
4104
=== modified file 'models/check_param_card.py'
4105
--- models/check_param_card.py	2018-04-26 14:24:04 +0000
4106
+++ models/check_param_card.py	2018-09-06 19:46:41 +0000
4107
@@ -566,6 +566,7 @@
4108
566
                return default
566
                return default
4109
567
            raise
567
            raise
4110
568
568
4111
569
<<<<<<< TREE
4112
569
    def get_missing_block(self, identpath):
570
    def get_missing_block(self, identpath):
4113
570
        """ """
571
        """ """
4114
571
        missing = set()
572
        missing = set()
4115
@@ -612,6 +613,9 @@
4116
612
            self.write(param_card.input_path)
613
            self.write(param_card.input_path)
4117
613
        
614
        
4118
614
    def write_inc_file(self, outpath, identpath, default, need_mp=False):
615
    def write_inc_file(self, outpath, identpath, default, need_mp=False):
4119
616
=======
4120
617
    def write_inc_file(self, outpath, identpath, default, decay_to_keep = [], need_mp=False):
4121
618
>>>>>>> MERGE-SOURCE
4122
615
        """ write a fortran file which hardcode the param value"""
619
        """ write a fortran file which hardcode the param value"""
4123
616
        
620
        
4124
617
        self.secure_slha2(identpath)
621
        self.secure_slha2(identpath)
4125
@@ -643,6 +647,10 @@
4126
643
            fout.writelines(' %s = %s' % (variable, ('%e'%float(value)).replace('e','d')))
647
            fout.writelines(' %s = %s' % (variable, ('%e'%float(value)).replace('e','d')))
4127
644
            if need_mp:
648
            if need_mp:
4128
645
                fout.writelines(' mp__%s = %s_16' % (variable, value))
649
                fout.writelines(' mp__%s = %s_16' % (variable, value))
4129
650
        for width, param in decay_to_keep:
4130
651
            fout.writelines(' mdl_%s_keep = %s' % (width, ('%e'%float(param.value)).replace('e','d')))
4131
652
        
4132
653
                
4133
646
      
654
      
4134
647
    def convert_to_complex_mass_scheme(self):
655
    def convert_to_complex_mass_scheme(self):
4135
648
        """ Convert this param_card to the convention used for the complex mass scheme:
656
        """ Convert this param_card to the convention used for the complex mass scheme:
4136
@@ -855,7 +863,7 @@
4137
855
class ParamCardMP(ParamCard):
863
class ParamCardMP(ParamCard):
4138
856
    """ a param Card: list of Block with also MP definition of variables"""
864
    """ a param Card: list of Block with also MP definition of variables"""
4139
857
            
865
            
4141
858
    def write_inc_file(self, outpath, identpath, default):
866
    def write_inc_file(self, outpath, identpath, default, decay_to_keep=[]):
4142
859
        """ write a fortran file which hardcode the param value"""
867
        """ write a fortran file which hardcode the param value"""
4143
860
        
868
        
4144
861
        fout = file_writers.FortranWriter(outpath)
869
        fout = file_writers.FortranWriter(outpath)
4145
@@ -880,6 +888,8 @@
4146
880
            fout.writelines(' %s = %s' % (variable, ('%e' % value).replace('e','d')))
888
            fout.writelines(' %s = %s' % (variable, ('%e' % value).replace('e','d')))
4147
881
            fout.writelines(' %s%s = %s_16' % (self.mp_prefix, 
889
            fout.writelines(' %s%s = %s_16' % (self.mp_prefix, 
4148
882
                variable, ('%e' % value)))
890
                variable, ('%e' % value)))
4149
891
        for width, param in decay_to_keep:
4150
892
            fout.writelines(' mdl_%s_keep = %s' % (width, ('%e'%float(param.value)).replace('e','d')))
4151
883
893
4152
884
894
4153
885
  
895
  
4154
886
896
4155
=== added directory 'tests/input_files/IOTestsComparison.moved'
4156
=== added directory 'tests/input_files/IOTestsComparison.moved/IOExportFKSTest'
4157
=== added directory 'tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal'
4158
=== added file 'tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_001.f'
4159
--- tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_001.f	1970-01-01 00:00:00 +0000
4160
+++ tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_001.f	2018-09-06 19:46:41 +0000
4161
@@ -0,0 +1,156 @@
4162
1
      SUBROUTINE SB_SF_001(P1,ANS)
4163
2
C     
4164
3
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4165
4
C     By the MadGraph5_aMC@NLO Development Team
4166
5
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4167
6
C     
4168
7
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4169
8
C     AND HELICITIES
4170
9
C     FOR THE POINT IN PHASE SPACE P(0:3,NEXTERNAL-1)
4171
10
C     
4172
11
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4173
12
C     spectators: 1 2 
4174
13
4175
14
C     
4176
15
      IMPLICIT NONE
4177
16
C     
4178
17
C     CONSTANTS
4179
18
C     
4180
19
      INCLUDE 'nexternal.inc'
4181
20
      INTEGER     NCOMB
4182
21
      PARAMETER ( NCOMB=  16 )
4183
22
      INTEGER    THEL
4184
23
      PARAMETER (THEL=NCOMB*8)
4185
24
      INTEGER NGRAPHS
4186
25
      PARAMETER (NGRAPHS=   3)
4187
26
C     
4188
27
C     ARGUMENTS 
4189
28
C     
4190
29
      REAL*8 P1(0:3,NEXTERNAL-1),ANS
4191
30
C     
4192
31
C     LOCAL VARIABLES 
4193
32
C     
4194
33
      INTEGER IHEL,IDEN,J
4195
34
      REAL*8 B_SF_001
4196
35
      INTEGER IDEN_VALUES(8)
4197
36
      DATA IDEN_VALUES /256, 256, 256, 256, 256, 256, 256, 256/
4198
37
C     
4199
38
C     GLOBAL VARIABLES
4200
39
C     
4201
40
      LOGICAL GOODHEL(NCOMB,8)
4202
41
      COMMON /C_GOODHEL/ GOODHEL
4203
42
      DOUBLE PRECISION SAVEMOM(NEXTERNAL-1,2)
4204
43
      COMMON/TO_SAVEMOM/SAVEMOM
4205
44
      LOGICAL CALCULATEDBORN
4206
45
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4207
46
      INTEGER NFKSPROCESS
4208
47
      COMMON/C_NFKSPROCESS/NFKSPROCESS
4209
48
C     ----------
4210
49
C     BEGIN CODE
4211
50
C     ----------
4212
51
      IDEN=IDEN_VALUES(NFKSPROCESS)
4213
52
      IF (CALCULATEDBORN) THEN
4214
53
        DO J=1,NEXTERNAL-1
4215
54
          IF (SAVEMOM(J,1).NE.P1(0,J) .OR. SAVEMOM(J,2).NE.P1(3
4216
55
     $     ,J)) THEN
4217
56
            CALCULATEDBORN=.FALSE.
4218
57
            WRITE(*,*) 'Error in sb_sf: momenta not the same in the
4219
58
     $        born'
4220
59
            STOP
4221
60
          ENDIF
4222
61
        ENDDO
4223
62
      ELSE
4224
63
        WRITE(*,*) 'Error in sb_sf: color_linked borns should be
4225
64
     $    called only with calculatedborn = true'
4226
65
        STOP
4227
66
      ENDIF
4228
67
      ANS = 0D0
4229
68
      DO IHEL=1,NCOMB
4230
69
        IF (GOODHEL(IHEL,NFKSPROCESS)) THEN
4231
70
          ANS=ANS+B_SF_001(P1,IHEL)
4232
71
        ENDIF
4233
72
      ENDDO
4234
73
      ANS=ANS/DBLE(IDEN)
4235
74
      END
4236
75
4237
76
4238
77
      REAL*8 FUNCTION B_SF_001(P,HELL)
4239
78
C     
4240
79
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4241
80
C     By the MadGraph5_aMC@NLO Development Team
4242
81
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4243
82
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4244
83
C     FOR THE POINT WITH EXTERNAL LINES W(0:6,NEXTERNAL-1)
4245
84
4246
85
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4247
86
C     spectators: 1 2 
4248
87
4249
88
C     
4250
89
      IMPLICIT NONE
4251
90
C     
4252
91
C     CONSTANTS
4253
92
C     
4254
93
      INTEGER     NGRAPHS
4255
94
      PARAMETER ( NGRAPHS = 3 )
4256
95
      INTEGER NCOLOR1, NCOLOR2
4257
96
      PARAMETER (NCOLOR1=2, NCOLOR2=3)
4258
97
      REAL*8     ZERO
4259
98
      PARAMETER (ZERO=0D0)
4260
99
      COMPLEX*16 IMAG1
4261
100
      PARAMETER (IMAG1 = (0D0,1D0))
4262
101
      INCLUDE 'nexternal.inc'
4263
102
      INCLUDE 'born_nhel.inc'
4264
103
C     
4265
104
C     ARGUMENTS 
4266
105
C     
4267
106
      REAL*8 P(0:3,NEXTERNAL-1)
4268
107
      INTEGER HELL
4269
108
C     
4270
109
C     LOCAL VARIABLES 
4271
110
C     
4272
111
      INTEGER I,J
4273
112
      REAL*8 DENOM(NCOLOR1), CF(NCOLOR2,NCOLOR1)
4274
113
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP1(NCOLOR1), JAMP2(NCOLOR2)
4275
114
C     
4276
115
C     GLOBAL VARIABLES
4277
116
C     
4278
117
      DOUBLE COMPLEX SAVEAMP(NGRAPHS,MAX_BHEL)
4279
118
      COMMON/TO_SAVEAMP/SAVEAMP
4280
119
      LOGICAL CALCULATEDBORN
4281
120
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4282
121
C     
4283
122
C     COLOR DATA
4284
123
C     
4285
124
      DATA DENOM(1)/3/
4286
125
      DATA (CF(I,  1),I=  1,  3) /   16,   -2,    6/
4287
126
      DATA DENOM(2)/3/
4288
127
      DATA (CF(I,  2),I=  1,  3) /   -2,   16,    6/
4289
128
C     ----------
4290
129
C     BEGIN CODE
4291
130
C     ----------
4292
131
      IF (.NOT. CALCULATEDBORN) THEN
4293
132
        WRITE(*,*) 'Error in b_sf: color_linked borns should be called
4294
133
     $    only with calculatedborn = true'
4295
134
        STOP
4296
135
      ELSEIF (CALCULATEDBORN) THEN
4297
136
        DO I=1,NGRAPHS
4298
137
          AMP(I)=SAVEAMP(I,HELL)
4299
138
        ENDDO
4300
139
      ENDIF
4301
140
      JAMP1(1)=+IMAG1*AMP(1)-AMP(2)
4302
141
      JAMP1(2)=-IMAG1*AMP(1)-AMP(3)
4303
142
      JAMP2(1)=+1D0/2D0*(-3D0*IMAG1*AMP(1)+3D0*AMP(2))
4304
143
      JAMP2(2)=+1D0/2D0*(+3D0*IMAG1*AMP(1)+3D0*AMP(3))
4305
144
      JAMP2(3)=+1D0/2D0*(+AMP(2)+AMP(3))
4306
145
      B_SF_001 = 0.D0
4307
146
      DO I = 1, NCOLOR1
4308
147
        ZTEMP = (0.D0,0.D0)
4309
148
        DO J = 1, NCOLOR2
4310
149
          ZTEMP = ZTEMP + CF(J,I)*JAMP2(J)
4311
150
        ENDDO
4312
151
        B_SF_001 =B_SF_001+ZTEMP*DCONJG(JAMP1(I))/DENOM(I)
4313
152
      ENDDO
4314
153
      END
4315
154
4316
155
4317
156
4318
0
157
4319
=== added file 'tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_002.f'
4320
--- tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_002.f	1970-01-01 00:00:00 +0000
4321
+++ tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_002.f	2018-09-06 19:46:41 +0000
4322
@@ -0,0 +1,155 @@
4323
1
      SUBROUTINE SB_SF_002(P1,ANS)
4324
2
C     
4325
3
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4326
4
C     By the MadGraph5_aMC@NLO Development Team
4327
5
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4328
6
C     
4329
7
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4330
8
C     AND HELICITIES
4331
9
C     FOR THE POINT IN PHASE SPACE P(0:3,NEXTERNAL-1)
4332
10
C     
4333
11
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4334
12
C     spectators: 1 3 
4335
13
4336
14
C     
4337
15
      IMPLICIT NONE
4338
16
C     
4339
17
C     CONSTANTS
4340
18
C     
4341
19
      INCLUDE 'nexternal.inc'
4342
20
      INTEGER     NCOMB
4343
21
      PARAMETER ( NCOMB=  16 )
4344
22
      INTEGER    THEL
4345
23
      PARAMETER (THEL=NCOMB*8)
4346
24
      INTEGER NGRAPHS
4347
25
      PARAMETER (NGRAPHS=   3)
4348
26
C     
4349
27
C     ARGUMENTS 
4350
28
C     
4351
29
      REAL*8 P1(0:3,NEXTERNAL-1),ANS
4352
30
C     
4353
31
C     LOCAL VARIABLES 
4354
32
C     
4355
33
      INTEGER IHEL,IDEN,J
4356
34
      REAL*8 B_SF_002
4357
35
      INTEGER IDEN_VALUES(8)
4358
36
      DATA IDEN_VALUES /256, 256, 256, 256, 256, 256, 256, 256/
4359
37
C     
4360
38
C     GLOBAL VARIABLES
4361
39
C     
4362
40
      LOGICAL GOODHEL(NCOMB,8)
4363
41
      COMMON /C_GOODHEL/ GOODHEL
4364
42
      DOUBLE PRECISION SAVEMOM(NEXTERNAL-1,2)
4365
43
      COMMON/TO_SAVEMOM/SAVEMOM
4366
44
      LOGICAL CALCULATEDBORN
4367
45
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4368
46
      INTEGER NFKSPROCESS
4369
47
      COMMON/C_NFKSPROCESS/NFKSPROCESS
4370
48
C     ----------
4371
49
C     BEGIN CODE
4372
50
C     ----------
4373
51
      IDEN=IDEN_VALUES(NFKSPROCESS)
4374
52
      IF (CALCULATEDBORN) THEN
4375
53
        DO J=1,NEXTERNAL-1
4376
54
          IF (SAVEMOM(J,1).NE.P1(0,J) .OR. SAVEMOM(J,2).NE.P1(3
4377
55
     $     ,J)) THEN
4378
56
            CALCULATEDBORN=.FALSE.
4379
57
            WRITE(*,*) 'Error in sb_sf: momenta not the same in the
4380
58
     $        born'
4381
59
            STOP
4382
60
          ENDIF
4383
61
        ENDDO
4384
62
      ELSE
4385
63
        WRITE(*,*) 'Error in sb_sf: color_linked borns should be
4386
64
     $    called only with calculatedborn = true'
4387
65
        STOP
4388
66
      ENDIF
4389
67
      ANS = 0D0
4390
68
      DO IHEL=1,NCOMB
4391
69
        IF (GOODHEL(IHEL,NFKSPROCESS)) THEN
4392
70
          ANS=ANS+B_SF_002(P1,IHEL)
4393
71
        ENDIF
4394
72
      ENDDO
4395
73
      ANS=ANS/DBLE(IDEN)
4396
74
      END
4397
75
4398
76
4399
77
      REAL*8 FUNCTION B_SF_002(P,HELL)
4400
78
C     
4401
79
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4402
80
C     By the MadGraph5_aMC@NLO Development Team
4403
81
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4404
82
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4405
83
C     FOR THE POINT WITH EXTERNAL LINES W(0:6,NEXTERNAL-1)
4406
84
4407
85
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4408
86
C     spectators: 1 3 
4409
87
4410
88
C     
4411
89
      IMPLICIT NONE
4412
90
C     
4413
91
C     CONSTANTS
4414
92
C     
4415
93
      INTEGER     NGRAPHS
4416
94
      PARAMETER ( NGRAPHS = 3 )
4417
95
      INTEGER NCOLOR1, NCOLOR2
4418
96
      PARAMETER (NCOLOR1=2, NCOLOR2=2)
4419
97
      REAL*8     ZERO
4420
98
      PARAMETER (ZERO=0D0)
4421
99
      COMPLEX*16 IMAG1
4422
100
      PARAMETER (IMAG1 = (0D0,1D0))
4423
101
      INCLUDE 'nexternal.inc'
4424
102
      INCLUDE 'born_nhel.inc'
4425
103
C     
4426
104
C     ARGUMENTS 
4427
105
C     
4428
106
      REAL*8 P(0:3,NEXTERNAL-1)
4429
107
      INTEGER HELL
4430
108
C     
4431
109
C     LOCAL VARIABLES 
4432
110
C     
4433
111
      INTEGER I,J
4434
112
      REAL*8 DENOM(NCOLOR1), CF(NCOLOR2,NCOLOR1)
4435
113
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP1(NCOLOR1), JAMP2(NCOLOR2)
4436
114
C     
4437
115
C     GLOBAL VARIABLES
4438
116
C     
4439
117
      DOUBLE COMPLEX SAVEAMP(NGRAPHS,MAX_BHEL)
4440
118
      COMMON/TO_SAVEAMP/SAVEAMP
4441
119
      LOGICAL CALCULATEDBORN
4442
120
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4443
121
C     
4444
122
C     COLOR DATA
4445
123
C     
4446
124
      DATA DENOM(1)/3/
4447
125
      DATA (CF(I,  1),I=  1,  2) /   16,    6/
4448
126
      DATA DENOM(2)/3/
4449
127
      DATA (CF(I,  2),I=  1,  2) /   -2,    6/
4450
128
C     ----------
4451
129
C     BEGIN CODE
4452
130
C     ----------
4453
131
      IF (.NOT. CALCULATEDBORN) THEN
4454
132
        WRITE(*,*) 'Error in b_sf: color_linked borns should be called
4455
133
     $    only with calculatedborn = true'
4456
134
        STOP
4457
135
      ELSEIF (CALCULATEDBORN) THEN
4458
136
        DO I=1,NGRAPHS
4459
137
          AMP(I)=SAVEAMP(I,HELL)
4460
138
        ENDDO
4461
139
      ENDIF
4462
140
      JAMP1(1)=+IMAG1*AMP(1)-AMP(2)
4463
141
      JAMP1(2)=-IMAG1*AMP(1)-AMP(3)
4464
142
      JAMP2(1)=+1D0/2D0*(-3D0*IMAG1*AMP(1)+3D0*AMP(2))
4465
143
      JAMP2(2)=+1D0/2D0*(-IMAG1*AMP(1)-AMP(3))
4466
144
      B_SF_002 = 0.D0
4467
145
      DO I = 1, NCOLOR1
4468
146
        ZTEMP = (0.D0,0.D0)
4469
147
        DO J = 1, NCOLOR2
4470
148
          ZTEMP = ZTEMP + CF(J,I)*JAMP2(J)
4471
149
        ENDDO
4472
150
        B_SF_002 =B_SF_002+ZTEMP*DCONJG(JAMP1(I))/DENOM(I)
4473
151
      ENDDO
4474
152
      END
4475
153
4476
154
4477
155
4478
0
156
4479
=== added file 'tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_003.f'
4480
--- tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_003.f	1970-01-01 00:00:00 +0000
4481
+++ tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_003.f	2018-09-06 19:46:41 +0000
4482
@@ -0,0 +1,155 @@
4483
1
      SUBROUTINE SB_SF_003(P1,ANS)
4484
2
C     
4485
3
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4486
4
C     By the MadGraph5_aMC@NLO Development Team
4487
5
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4488
6
C     
4489
7
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4490
8
C     AND HELICITIES
4491
9
C     FOR THE POINT IN PHASE SPACE P(0:3,NEXTERNAL-1)
4492
10
C     
4493
11
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4494
12
C     spectators: 1 4 
4495
13
4496
14
C     
4497
15
      IMPLICIT NONE
4498
16
C     
4499
17
C     CONSTANTS
4500
18
C     
4501
19
      INCLUDE 'nexternal.inc'
4502
20
      INTEGER     NCOMB
4503
21
      PARAMETER ( NCOMB=  16 )
4504
22
      INTEGER    THEL
4505
23
      PARAMETER (THEL=NCOMB*8)
4506
24
      INTEGER NGRAPHS
4507
25
      PARAMETER (NGRAPHS=   3)
4508
26
C     
4509
27
C     ARGUMENTS 
4510
28
C     
4511
29
      REAL*8 P1(0:3,NEXTERNAL-1),ANS
4512
30
C     
4513
31
C     LOCAL VARIABLES 
4514
32
C     
4515
33
      INTEGER IHEL,IDEN,J
4516
34
      REAL*8 B_SF_003
4517
35
      INTEGER IDEN_VALUES(8)
4518
36
      DATA IDEN_VALUES /256, 256, 256, 256, 256, 256, 256, 256/
4519
37
C     
4520
38
C     GLOBAL VARIABLES
4521
39
C     
4522
40
      LOGICAL GOODHEL(NCOMB,8)
4523
41
      COMMON /C_GOODHEL/ GOODHEL
4524
42
      DOUBLE PRECISION SAVEMOM(NEXTERNAL-1,2)
4525
43
      COMMON/TO_SAVEMOM/SAVEMOM
4526
44
      LOGICAL CALCULATEDBORN
4527
45
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4528
46
      INTEGER NFKSPROCESS
4529
47
      COMMON/C_NFKSPROCESS/NFKSPROCESS
4530
48
C     ----------
4531
49
C     BEGIN CODE
4532
50
C     ----------
4533
51
      IDEN=IDEN_VALUES(NFKSPROCESS)
4534
52
      IF (CALCULATEDBORN) THEN
4535
53
        DO J=1,NEXTERNAL-1
4536
54
          IF (SAVEMOM(J,1).NE.P1(0,J) .OR. SAVEMOM(J,2).NE.P1(3
4537
55
     $     ,J)) THEN
4538
56
            CALCULATEDBORN=.FALSE.
4539
57
            WRITE(*,*) 'Error in sb_sf: momenta not the same in the
4540
58
     $        born'
4541
59
            STOP
4542
60
          ENDIF
4543
61
        ENDDO
4544
62
      ELSE
4545
63
        WRITE(*,*) 'Error in sb_sf: color_linked borns should be
4546
64
     $    called only with calculatedborn = true'
4547
65
        STOP
4548
66
      ENDIF
4549
67
      ANS = 0D0
4550
68
      DO IHEL=1,NCOMB
4551
69
        IF (GOODHEL(IHEL,NFKSPROCESS)) THEN
4552
70
          ANS=ANS+B_SF_003(P1,IHEL)
4553
71
        ENDIF
4554
72
      ENDDO
4555
73
      ANS=ANS/DBLE(IDEN)
4556
74
      END
4557
75
4558
76
4559
77
      REAL*8 FUNCTION B_SF_003(P,HELL)
4560
78
C     
4561
79
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4562
80
C     By the MadGraph5_aMC@NLO Development Team
4563
81
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4564
82
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4565
83
C     FOR THE POINT WITH EXTERNAL LINES W(0:6,NEXTERNAL-1)
4566
84
4567
85
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4568
86
C     spectators: 1 4 
4569
87
4570
88
C     
4571
89
      IMPLICIT NONE
4572
90
C     
4573
91
C     CONSTANTS
4574
92
C     
4575
93
      INTEGER     NGRAPHS
4576
94
      PARAMETER ( NGRAPHS = 3 )
4577
95
      INTEGER NCOLOR1, NCOLOR2
4578
96
      PARAMETER (NCOLOR1=2, NCOLOR2=2)
4579
97
      REAL*8     ZERO
4580
98
      PARAMETER (ZERO=0D0)
4581
99
      COMPLEX*16 IMAG1
4582
100
      PARAMETER (IMAG1 = (0D0,1D0))
4583
101
      INCLUDE 'nexternal.inc'
4584
102
      INCLUDE 'born_nhel.inc'
4585
103
C     
4586
104
C     ARGUMENTS 
4587
105
C     
4588
106
      REAL*8 P(0:3,NEXTERNAL-1)
4589
107
      INTEGER HELL
4590
108
C     
4591
109
C     LOCAL VARIABLES 
4592
110
C     
4593
111
      INTEGER I,J
4594
112
      REAL*8 DENOM(NCOLOR1), CF(NCOLOR2,NCOLOR1)
4595
113
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP1(NCOLOR1), JAMP2(NCOLOR2)
4596
114
C     
4597
115
C     GLOBAL VARIABLES
4598
116
C     
4599
117
      DOUBLE COMPLEX SAVEAMP(NGRAPHS,MAX_BHEL)
4600
118
      COMMON/TO_SAVEAMP/SAVEAMP
4601
119
      LOGICAL CALCULATEDBORN
4602
120
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4603
121
C     
4604
122
C     COLOR DATA
4605
123
C     
4606
124
      DATA DENOM(1)/3/
4607
125
      DATA (CF(I,  1),I=  1,  2) /   -2,    6/
4608
126
      DATA DENOM(2)/3/
4609
127
      DATA (CF(I,  2),I=  1,  2) /   16,    6/
4610
128
C     ----------
4611
129
C     BEGIN CODE
4612
130
C     ----------
4613
131
      IF (.NOT. CALCULATEDBORN) THEN
4614
132
        WRITE(*,*) 'Error in b_sf: color_linked borns should be called
4615
133
     $    only with calculatedborn = true'
4616
134
        STOP
4617
135
      ELSEIF (CALCULATEDBORN) THEN
4618
136
        DO I=1,NGRAPHS
4619
137
          AMP(I)=SAVEAMP(I,HELL)
4620
138
        ENDDO
4621
139
      ENDIF
4622
140
      JAMP1(1)=+IMAG1*AMP(1)-AMP(2)
4623
141
      JAMP1(2)=-IMAG1*AMP(1)-AMP(3)
4624
142
      JAMP2(1)=+1D0/2D0*(+3D0*IMAG1*AMP(1)+3D0*AMP(3))
4625
143
      JAMP2(2)=+1D0/2D0*(+IMAG1*AMP(1)-AMP(2))
4626
144
      B_SF_003 = 0.D0
4627
145
      DO I = 1, NCOLOR1
4628
146
        ZTEMP = (0.D0,0.D0)
4629
147
        DO J = 1, NCOLOR2
4630
148
          ZTEMP = ZTEMP + CF(J,I)*JAMP2(J)
4631
149
        ENDDO
4632
150
        B_SF_003 =B_SF_003+ZTEMP*DCONJG(JAMP1(I))/DENOM(I)
4633
151
      ENDDO
4634
152
      END
4635
153
4636
154
4637
155
4638
0
156
4639
=== added file 'tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_004.f'
4640
--- tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_004.f	1970-01-01 00:00:00 +0000
4641
+++ tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_004.f	2018-09-06 19:46:41 +0000
4642
@@ -0,0 +1,155 @@
4643
1
      SUBROUTINE SB_SF_004(P1,ANS)
4644
2
C     
4645
3
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4646
4
C     By the MadGraph5_aMC@NLO Development Team
4647
5
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4648
6
C     
4649
7
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4650
8
C     AND HELICITIES
4651
9
C     FOR THE POINT IN PHASE SPACE P(0:3,NEXTERNAL-1)
4652
10
C     
4653
11
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4654
12
C     spectators: 2 3 
4655
13
4656
14
C     
4657
15
      IMPLICIT NONE
4658
16
C     
4659
17
C     CONSTANTS
4660
18
C     
4661
19
      INCLUDE 'nexternal.inc'
4662
20
      INTEGER     NCOMB
4663
21
      PARAMETER ( NCOMB=  16 )
4664
22
      INTEGER    THEL
4665
23
      PARAMETER (THEL=NCOMB*8)
4666
24
      INTEGER NGRAPHS
4667
25
      PARAMETER (NGRAPHS=   3)
4668
26
C     
4669
27
C     ARGUMENTS 
4670
28
C     
4671
29
      REAL*8 P1(0:3,NEXTERNAL-1),ANS
4672
30
C     
4673
31
C     LOCAL VARIABLES 
4674
32
C     
4675
33
      INTEGER IHEL,IDEN,J
4676
34
      REAL*8 B_SF_004
4677
35
      INTEGER IDEN_VALUES(8)
4678
36
      DATA IDEN_VALUES /256, 256, 256, 256, 256, 256, 256, 256/
4679
37
C     
4680
38
C     GLOBAL VARIABLES
4681
39
C     
4682
40
      LOGICAL GOODHEL(NCOMB,8)
4683
41
      COMMON /C_GOODHEL/ GOODHEL
4684
42
      DOUBLE PRECISION SAVEMOM(NEXTERNAL-1,2)
4685
43
      COMMON/TO_SAVEMOM/SAVEMOM
4686
44
      LOGICAL CALCULATEDBORN
4687
45
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4688
46
      INTEGER NFKSPROCESS
4689
47
      COMMON/C_NFKSPROCESS/NFKSPROCESS
4690
48
C     ----------
4691
49
C     BEGIN CODE
4692
50
C     ----------
4693
51
      IDEN=IDEN_VALUES(NFKSPROCESS)
4694
52
      IF (CALCULATEDBORN) THEN
4695
53
        DO J=1,NEXTERNAL-1
4696
54
          IF (SAVEMOM(J,1).NE.P1(0,J) .OR. SAVEMOM(J,2).NE.P1(3
4697
55
     $     ,J)) THEN
4698
56
            CALCULATEDBORN=.FALSE.
4699
57
            WRITE(*,*) 'Error in sb_sf: momenta not the same in the
4700
58
     $        born'
4701
59
            STOP
4702
60
          ENDIF
4703
61
        ENDDO
4704
62
      ELSE
4705
63
        WRITE(*,*) 'Error in sb_sf: color_linked borns should be
4706
64
     $    called only with calculatedborn = true'
4707
65
        STOP
4708
66
      ENDIF
4709
67
      ANS = 0D0
4710
68
      DO IHEL=1,NCOMB
4711
69
        IF (GOODHEL(IHEL,NFKSPROCESS)) THEN
4712
70
          ANS=ANS+B_SF_004(P1,IHEL)
4713
71
        ENDIF
4714
72
      ENDDO
4715
73
      ANS=ANS/DBLE(IDEN)
4716
74
      END
4717
75
4718
76
4719
77
      REAL*8 FUNCTION B_SF_004(P,HELL)
4720
78
C     
4721
79
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4722
80
C     By the MadGraph5_aMC@NLO Development Team
4723
81
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4724
82
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4725
83
C     FOR THE POINT WITH EXTERNAL LINES W(0:6,NEXTERNAL-1)
4726
84
4727
85
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4728
86
C     spectators: 2 3 
4729
87
4730
88
C     
4731
89
      IMPLICIT NONE
4732
90
C     
4733
91
C     CONSTANTS
4734
92
C     
4735
93
      INTEGER     NGRAPHS
4736
94
      PARAMETER ( NGRAPHS = 3 )
4737
95
      INTEGER NCOLOR1, NCOLOR2
4738
96
      PARAMETER (NCOLOR1=2, NCOLOR2=2)
4739
97
      REAL*8     ZERO
4740
98
      PARAMETER (ZERO=0D0)
4741
99
      COMPLEX*16 IMAG1
4742
100
      PARAMETER (IMAG1 = (0D0,1D0))
4743
101
      INCLUDE 'nexternal.inc'
4744
102
      INCLUDE 'born_nhel.inc'
4745
103
C     
4746
104
C     ARGUMENTS 
4747
105
C     
4748
106
      REAL*8 P(0:3,NEXTERNAL-1)
4749
107
      INTEGER HELL
4750
108
C     
4751
109
C     LOCAL VARIABLES 
4752
110
C     
4753
111
      INTEGER I,J
4754
112
      REAL*8 DENOM(NCOLOR1), CF(NCOLOR2,NCOLOR1)
4755
113
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP1(NCOLOR1), JAMP2(NCOLOR2)
4756
114
C     
4757
115
C     GLOBAL VARIABLES
4758
116
C     
4759
117
      DOUBLE COMPLEX SAVEAMP(NGRAPHS,MAX_BHEL)
4760
118
      COMMON/TO_SAVEAMP/SAVEAMP
4761
119
      LOGICAL CALCULATEDBORN
4762
120
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4763
121
C     
4764
122
C     COLOR DATA
4765
123
C     
4766
124
      DATA DENOM(1)/3/
4767
125
      DATA (CF(I,  1),I=  1,  2) /   -2,    6/
4768
126
      DATA DENOM(2)/3/
4769
127
      DATA (CF(I,  2),I=  1,  2) /   16,    6/
4770
128
C     ----------
4771
129
C     BEGIN CODE
4772
130
C     ----------
4773
131
      IF (.NOT. CALCULATEDBORN) THEN
4774
132
        WRITE(*,*) 'Error in b_sf: color_linked borns should be called
4775
133
     $    only with calculatedborn = true'
4776
134
        STOP
4777
135
      ELSEIF (CALCULATEDBORN) THEN
4778
136
        DO I=1,NGRAPHS
4779
137
          AMP(I)=SAVEAMP(I,HELL)
4780
138
        ENDDO
4781
139
      ENDIF
4782
140
      JAMP1(1)=+IMAG1*AMP(1)-AMP(2)
4783
141
      JAMP1(2)=-IMAG1*AMP(1)-AMP(3)
4784
142
      JAMP2(1)=+1D0/2D0*(+3D0*IMAG1*AMP(1)+3D0*AMP(3))
4785
143
      JAMP2(2)=+1D0/2D0*(+IMAG1*AMP(1)-AMP(2))
4786
144
      B_SF_004 = 0.D0
4787
145
      DO I = 1, NCOLOR1
4788
146
        ZTEMP = (0.D0,0.D0)
4789
147
        DO J = 1, NCOLOR2
4790
148
          ZTEMP = ZTEMP + CF(J,I)*JAMP2(J)
4791
149
        ENDDO
4792
150
        B_SF_004 =B_SF_004+ZTEMP*DCONJG(JAMP1(I))/DENOM(I)
4793
151
      ENDDO
4794
152
      END
4795
153
4796
154
4797
155
4798
0
156
4799
=== added file 'tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_005.f'
4800
--- tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_005.f	1970-01-01 00:00:00 +0000
4801
+++ tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_005.f	2018-09-06 19:46:41 +0000
4802
@@ -0,0 +1,155 @@
4803
1
      SUBROUTINE SB_SF_005(P1,ANS)
4804
2
C     
4805
3
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4806
4
C     By the MadGraph5_aMC@NLO Development Team
4807
5
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4808
6
C     
4809
7
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4810
8
C     AND HELICITIES
4811
9
C     FOR THE POINT IN PHASE SPACE P(0:3,NEXTERNAL-1)
4812
10
C     
4813
11
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4814
12
C     spectators: 2 4 
4815
13
4816
14
C     
4817
15
      IMPLICIT NONE
4818
16
C     
4819
17
C     CONSTANTS
4820
18
C     
4821
19
      INCLUDE 'nexternal.inc'
4822
20
      INTEGER     NCOMB
4823
21
      PARAMETER ( NCOMB=  16 )
4824
22
      INTEGER    THEL
4825
23
      PARAMETER (THEL=NCOMB*8)
4826
24
      INTEGER NGRAPHS
4827
25
      PARAMETER (NGRAPHS=   3)
4828
26
C     
4829
27
C     ARGUMENTS 
4830
28
C     
4831
29
      REAL*8 P1(0:3,NEXTERNAL-1),ANS
4832
30
C     
4833
31
C     LOCAL VARIABLES 
4834
32
C     
4835
33
      INTEGER IHEL,IDEN,J
4836
34
      REAL*8 B_SF_005
4837
35
      INTEGER IDEN_VALUES(8)
4838
36
      DATA IDEN_VALUES /256, 256, 256, 256, 256, 256, 256, 256/
4839
37
C     
4840
38
C     GLOBAL VARIABLES
4841
39
C     
4842
40
      LOGICAL GOODHEL(NCOMB,8)
4843
41
      COMMON /C_GOODHEL/ GOODHEL
4844
42
      DOUBLE PRECISION SAVEMOM(NEXTERNAL-1,2)
4845
43
      COMMON/TO_SAVEMOM/SAVEMOM
4846
44
      LOGICAL CALCULATEDBORN
4847
45
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4848
46
      INTEGER NFKSPROCESS
4849
47
      COMMON/C_NFKSPROCESS/NFKSPROCESS
4850
48
C     ----------
4851
49
C     BEGIN CODE
4852
50
C     ----------
4853
51
      IDEN=IDEN_VALUES(NFKSPROCESS)
4854
52
      IF (CALCULATEDBORN) THEN
4855
53
        DO J=1,NEXTERNAL-1
4856
54
          IF (SAVEMOM(J,1).NE.P1(0,J) .OR. SAVEMOM(J,2).NE.P1(3
4857
55
     $     ,J)) THEN
4858
56
            CALCULATEDBORN=.FALSE.
4859
57
            WRITE(*,*) 'Error in sb_sf: momenta not the same in the
4860
58
     $        born'
4861
59
            STOP
4862
60
          ENDIF
4863
61
        ENDDO
4864
62
      ELSE
4865
63
        WRITE(*,*) 'Error in sb_sf: color_linked borns should be
4866
64
     $    called only with calculatedborn = true'
4867
65
        STOP
4868
66
      ENDIF
4869
67
      ANS = 0D0
4870
68
      DO IHEL=1,NCOMB
4871
69
        IF (GOODHEL(IHEL,NFKSPROCESS)) THEN
4872
70
          ANS=ANS+B_SF_005(P1,IHEL)
4873
71
        ENDIF
4874
72
      ENDDO
4875
73
      ANS=ANS/DBLE(IDEN)
4876
74
      END
4877
75
4878
76
4879
77
      REAL*8 FUNCTION B_SF_005(P,HELL)
4880
78
C     
4881
79
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4882
80
C     By the MadGraph5_aMC@NLO Development Team
4883
81
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4884
82
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4885
83
C     FOR THE POINT WITH EXTERNAL LINES W(0:6,NEXTERNAL-1)
4886
84
4887
85
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4888
86
C     spectators: 2 4 
4889
87
4890
88
C     
4891
89
      IMPLICIT NONE
4892
90
C     
4893
91
C     CONSTANTS
4894
92
C     
4895
93
      INTEGER     NGRAPHS
4896
94
      PARAMETER ( NGRAPHS = 3 )
4897
95
      INTEGER NCOLOR1, NCOLOR2
4898
96
      PARAMETER (NCOLOR1=2, NCOLOR2=2)
4899
97
      REAL*8     ZERO
4900
98
      PARAMETER (ZERO=0D0)
4901
99
      COMPLEX*16 IMAG1
4902
100
      PARAMETER (IMAG1 = (0D0,1D0))
4903
101
      INCLUDE 'nexternal.inc'
4904
102
      INCLUDE 'born_nhel.inc'
4905
103
C     
4906
104
C     ARGUMENTS 
4907
105
C     
4908
106
      REAL*8 P(0:3,NEXTERNAL-1)
4909
107
      INTEGER HELL
4910
108
C     
4911
109
C     LOCAL VARIABLES 
4912
110
C     
4913
111
      INTEGER I,J
4914
112
      REAL*8 DENOM(NCOLOR1), CF(NCOLOR2,NCOLOR1)
4915
113
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP1(NCOLOR1), JAMP2(NCOLOR2)
4916
114
C     
4917
115
C     GLOBAL VARIABLES
4918
116
C     
4919
117
      DOUBLE COMPLEX SAVEAMP(NGRAPHS,MAX_BHEL)
4920
118
      COMMON/TO_SAVEAMP/SAVEAMP
4921
119
      LOGICAL CALCULATEDBORN
4922
120
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
4923
121
C     
4924
122
C     COLOR DATA
4925
123
C     
4926
124
      DATA DENOM(1)/3/
4927
125
      DATA (CF(I,  1),I=  1,  2) /   16,    6/
4928
126
      DATA DENOM(2)/3/
4929
127
      DATA (CF(I,  2),I=  1,  2) /   -2,    6/
4930
128
C     ----------
4931
129
C     BEGIN CODE
4932
130
C     ----------
4933
131
      IF (.NOT. CALCULATEDBORN) THEN
4934
132
        WRITE(*,*) 'Error in b_sf: color_linked borns should be called
4935
133
     $    only with calculatedborn = true'
4936
134
        STOP
4937
135
      ELSEIF (CALCULATEDBORN) THEN
4938
136
        DO I=1,NGRAPHS
4939
137
          AMP(I)=SAVEAMP(I,HELL)
4940
138
        ENDDO
4941
139
      ENDIF
4942
140
      JAMP1(1)=+IMAG1*AMP(1)-AMP(2)
4943
141
      JAMP1(2)=-IMAG1*AMP(1)-AMP(3)
4944
142
      JAMP2(1)=+1D0/2D0*(-3D0*IMAG1*AMP(1)+3D0*AMP(2))
4945
143
      JAMP2(2)=+1D0/2D0*(-IMAG1*AMP(1)-AMP(3))
4946
144
      B_SF_005 = 0.D0
4947
145
      DO I = 1, NCOLOR1
4948
146
        ZTEMP = (0.D0,0.D0)
4949
147
        DO J = 1, NCOLOR2
4950
148
          ZTEMP = ZTEMP + CF(J,I)*JAMP2(J)
4951
149
        ENDDO
4952
150
        B_SF_005 =B_SF_005+ZTEMP*DCONJG(JAMP1(I))/DENOM(I)
4953
151
      ENDDO
4954
152
      END
4955
153
4956
154
4957
155
4958
0
156
4959
=== added file 'tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_006.f'
4960
--- tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_006.f	1970-01-01 00:00:00 +0000
4961
+++ tests/input_files/IOTestsComparison.moved/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%b_sf_006.f	2018-09-06 19:46:41 +0000
4962
@@ -0,0 +1,157 @@
4963
1
      SUBROUTINE SB_SF_006(P1,ANS)
4964
2
C     
4965
3
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4966
4
C     By the MadGraph5_aMC@NLO Development Team
4967
5
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
4968
6
C     
4969
7
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
4970
8
C     AND HELICITIES
4971
9
C     FOR THE POINT IN PHASE SPACE P(0:3,NEXTERNAL-1)
4972
10
C     
4973
11
C     Process: g g > t t~ WEIGHTED=2 [ QCD ]
4974
12
C     spectators: 3 3 
4975
13
4976
14
C     
4977
15
      IMPLICIT NONE
4978
16
C     
4979
17
C     CONSTANTS
4980
18
C     
4981
19
      INCLUDE 'nexternal.inc'
4982
20
      INTEGER     NCOMB
4983
21
      PARAMETER ( NCOMB=  16 )
4984
22
      INTEGER    THEL
4985
23
      PARAMETER (THEL=NCOMB*8)
4986
24
      INTEGER NGRAPHS
4987
25
      PARAMETER (NGRAPHS=   3)
4988
26
C     
4989
27
C     ARGUMENTS 
4990
28
C     
4991
29
      REAL*8 P1(0:3,NEXTERNAL-1),ANS
4992
30
C     
4993
31
C     LOCAL VARIABLES 
4994
32
C     
4995
33
      INTEGER IHEL,IDEN,J
4996
34
      REAL*8 B_SF_006
4997
35
      INTEGER IDEN_VALUES(8)
4998
36
      DATA IDEN_VALUES /256, 256, 256, 256, 256, 256, 256, 256/
4999
37
C     
5000
38
C     GLOBAL VARIABLES
Reviewer	Review Type	Date Requested	Status
mg5core1		2018-09-06	Pending
Review via email: mp+354418@code.launchpad.net