← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

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

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/symmetry_fks_test_ME.f
 
1
      program symmetry
 
2
c*****************************************************************************
 
3
c     Given identical particles, and the configurations. This program identifies
 
4
c     identical configurations and specifies which ones can be skipped
 
5
c*****************************************************************************
 
6
      implicit none
 
7
c
 
8
c     Constants
 
9
c
 
10
      include 'genps.inc'      
 
11
      include 'nexternal.inc'
 
12
      include '../../Source/run_config.inc'
 
13
      include 'nFKSconfigs.inc'
 
14
      
 
15
      double precision ZERO,one
 
16
      parameter       (ZERO = 0d0)
 
17
      parameter       (one = 1d0)
 
18
      integer   maxswitch
 
19
      parameter(maxswitch=99)
 
20
c
 
21
c     Local
 
22
c
 
23
      integer iforest(2,-max_branch:-1,lmaxconfigs)
 
24
      integer mapconfig(0:lmaxconfigs)
 
25
      integer sprop(-max_branch:-1,lmaxconfigs)
 
26
      integer itree(2,-max_branch:-1)
 
27
      integer imatch
 
28
      integer i,j, k, n, nsym,l,ii,jj
 
29
      double precision diff,xi_i_fks
 
30
c$$$      double precision pmass(-max_branch:-1,lmaxconfigs)   !Propagotor mass
 
31
      double precision pmass(nexternal)
 
32
      double precision pwidth(-max_branch:-1,lmaxconfigs)  !Propagotor width
 
33
      integer pow(-max_branch:-1,lmaxconfigs)
 
34
 
 
35
      integer biforest(2,-max_branch:-1,lmaxconfigs)
 
36
      integer fksmother,fksgrandmother,fksaunt,compare
 
37
      integer fksconfiguration,mapbconf(0:lmaxconfigs)
 
38
      integer r2b(lmaxconfigs),b2r(lmaxconfigs)
 
39
      logical searchforgranny,is_beta_cms,is_granny_sch,topdown,non_prop
 
40
      integer nbranch,ns_channel,nt_channel
 
41
c      include "fks.inc"
 
42
      integer fks_j_from_i(nexternal,0:nexternal)
 
43
     &     ,particle_type(nexternal),pdg_type(nexternal)
 
44
      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
 
45
      double precision fxl,limit(15),wlimit(15)
 
46
      double precision lxp(0:3,nexternal+1),xp(15,0:3,nexternal+1)
 
47
      double precision fks_Sij
 
48
      double precision check,tolerance,zh,h_damp
 
49
      parameter (tolerance=1.d-4)
 
50
      integer kk,ll,bs,bs_min,bs_max,iconfig_in
 
51
 
 
52
      integer nsofttests,ncolltests,nerr,imax,iflag,iret
 
53
c
 
54
c     Local for generating amps
 
55
c
 
56
      double precision p(0:3,99), wgt, x(99), fx
 
57
      double complex wgt1(2)
 
58
      double precision p1(0:3,99),xx(maxinvar)
 
59
      integer ninvar, ndim, iconfig, minconfig, maxconfig
 
60
      integer ncall,itmax,nconfigs,ntry, ngraphs
 
61
      integer ic(nexternal,maxswitch), jc(12),nswitch
 
62
      double precision saveamp(maxamps)
 
63
      integer nmatch, ibase
 
64
      logical mtc, even
 
65
 
 
66
 
 
67
      double precision xi_i_fks_fix_save,y_ij_fks_fix_save
 
68
      double precision xi_i_fks_fix,y_ij_fks_fix
 
69
      common/cxiyfix/xi_i_fks_fix,y_ij_fks_fix
 
70
c
 
71
c     Global
 
72
c
 
73
      Double Precision amp2(maxamps), jamp2(0:maxamps)
 
74
      common/to_amps/  amp2,       jamp2
 
75
      include 'coupl.inc'
 
76
 
 
77
      logical calculatedBorn
 
78
      common/ccalculatedBorn/calculatedBorn
 
79
 
 
80
      integer i_fks,j_fks
 
81
      common/fks_indices/i_fks,j_fks
 
82
 
 
83
      double precision p1_cnt(0:3,nexternal,-2:2)
 
84
      double precision wgt_cnt(-2:2)
 
85
      double precision pswgt_cnt(-2:2)
 
86
      double precision jac_cnt(-2:2)
 
87
      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
 
88
 
 
89
      double precision p_born(0:3,nexternal-1)
 
90
      common/pborn/p_born
 
91
 
 
92
      double precision xi_i_fks_ev,y_ij_fks_ev
 
93
      double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
 
94
      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
 
95
 
 
96
      double precision xi_i_fks_cnt(-2:2)
 
97
      common /cxiifkscnt/xi_i_fks_cnt
 
98
 
 
99
      logical rotategranny
 
100
      common/crotategranny/rotategranny
 
101
 
 
102
      logical softtest,colltest
 
103
      common/sctests/softtest,colltest
 
104
      
 
105
      logical xexternal
 
106
      common /toxexternal/ xexternal
 
107
 
 
108
c Particle types (=color) of i_fks, j_fks and fks_mother
 
109
      integer i_type,j_type,m_type
 
110
      common/cparticle_types/i_type,j_type,m_type
 
111
 
 
112
c
 
113
c     External
 
114
c
 
115
      logical pass_point
 
116
      logical check_swap
 
117
      double precision dsig,ran2
 
118
      external pass_point, dsig,ran2
 
119
      external check_swap, fks_Sij
 
120
 
 
121
c define here the maximum fraction of failures to consider the test
 
122
c   passed
 
123
      double precision max_fail, fail_frac
 
124
      parameter (max_fail=0.3d0)
 
125
 
 
126
c helicity stuff
 
127
      integer          isum_hel
 
128
      logical                    multi_channel
 
129
      common/to_matrix/isum_hel, multi_channel
 
130
 
 
131
      integer fks_conf_number,fks_loop_min,fks_loop_max,fks_loop
 
132
      INTEGER NFKSPROCESS
 
133
      COMMON/C_NFKSPROCESS/NFKSPROCESS
 
134
      
 
135
c      integer icomp
 
136
c
 
137
c     DATA
 
138
c
 
139
      integer tprid(-max_branch:-1,lmaxconfigs)
 
140
      include 'born_conf.inc'
 
141
c-----
 
142
c  Begin Code
 
143
c-----
 
144
      write(*,*)'Enter xi_i, y_ij to be used in coll/soft tests'
 
145
      write(*,*)' Enter -2 to generate them randomly'
 
146
      read(*,*)xi_i_fks_fix_save,y_ij_fks_fix_save
 
147
 
 
148
      write(*,*)'Enter number of tests for soft and collinear limits'
 
149
      read(*,*)nsofttests,ncolltests
 
150
 
 
151
      write(*,*)'Sum over helicity (0), or random helicity (1)'
 
152
      read(*,*) isum_hel
 
153
 
 
154
      call setrun                !Sets up run parameters
 
155
      call setpara('param_card.dat')   !Sets up couplings and masses
 
156
      call setcuts               !Sets up cuts 
 
157
c
 
158
 
 
159
      write (*,*) 'Give FKS configuration number ("0" loops over all)'
 
160
      read (*,*) fks_conf_number
 
161
 
 
162
      if (fks_conf_number.eq.0) then
 
163
         fks_loop_min=1
 
164
         fks_loop_max=fks_configs
 
165
      else
 
166
         fks_loop_min=fks_conf_number
 
167
         fks_loop_max=fks_conf_number
 
168
      endif
 
169
 
 
170
      do fks_loop=fks_loop_min,fks_loop_max
 
171
         nFKSprocess=fks_loop
 
172
         write (*,*) ''
 
173
         write (*,*) '================================================='
 
174
         write (*,*) ''
 
175
         write (*,*) 'NEW FKS CONFIGURATION:'
 
176
 
 
177
         call fks_inc_chooser()
 
178
         call leshouche_inc_chooser()
 
179
         write (*,*) 'FKS configuration number is ',nFKSprocess
 
180
         write (*,*) 'FKS partons are: i=',i_fks,'  j=',j_fks
 
181
         write (*,*) 'with PDGs:       i=',PDG_type(i_fks),'  j='
 
182
     $        ,PDG_type(j_fks)
 
183
 
 
184
 
 
185
c
 
186
      ndim = 22
 
187
      ncall = 10000
 
188
      itmax = 10
 
189
      ninvar = 35
 
190
      nconfigs = 1
 
191
 
 
192
c Set color types of i_fks, j_fks and fks_mother.
 
193
      i_type=particle_type(i_fks)
 
194
      j_type=particle_type(j_fks)
 
195
      if (abs(i_type).eq.abs(j_type)) then
 
196
         m_type=8
 
197
         if ( (j_fks.le.nincoming .and.
 
198
     &        abs(i_type).eq.3 .and. j_type.ne.i_type) .or.
 
199
     &        (j_fks.gt.nincoming .and.
 
200
     &        abs(i_type).eq.3 .and. j_type.ne.-i_type)) then
 
201
            write(*,*)'Flavour mismatch #1 in setfksfactor',
 
202
     &           i_fks,j_fks,i_type,j_type
 
203
            stop
 
204
         endif
 
205
      elseif(abs(i_type).eq.3 .and. j_type.eq.8)then
 
206
         if(j_fks.le.nincoming)then
 
207
            m_type=-i_type
 
208
         else
 
209
            write (*,*) 'Error in setfksfactor: (i,j)=(q,g)'
 
210
            stop
 
211
         endif
 
212
      elseif(i_type.eq.8 .and. abs(j_type).eq.3)then
 
213
         if (j_fks.le.nincoming) then
 
214
            m_type=j_type
 
215
         else
 
216
            m_type=j_type
 
217
         endif
 
218
      else
 
219
         write(*,*)'Flavour mismatch #2 in setfksfactor',
 
220
     &        i_type,j_type,m_type
 
221
         stop
 
222
      endif
 
223
 
 
224
 
 
225
c     
 
226
c     Get momentum configuration
 
227
c
 
228
 
 
229
c Set xexternal to true to use the x's from external vegas in the
 
230
c x_to_f_arg subroutine
 
231
      xexternal=.true.
 
232
      
 
233
      write(*,*)'  '
 
234
      write(*,*)'  '
 
235
      write(*,*)'Enter graph number (iconfig), '
 
236
     &     //"'0' loops over all graphs"
 
237
      read(*,*)iconfig_in
 
238
      
 
239
      if (iconfig_in.eq.0) then
 
240
         bs_min=1
 
241
         bs_max=mapconfig(0)
 
242
      elseif (iconfig_in.eq.-1) then
 
243
         bs_min=1
 
244
         bs_max=1
 
245
      else
 
246
         bs_min=iconfig_in
 
247
         bs_max=iconfig_in
 
248
      endif
 
249
 
 
250
      do iconfig=bs_min,bs_max       ! Born configurations
 
251
      call setcuts
 
252
      call setfksfactor(iconfig)
 
253
      wgt=1d0
 
254
      ntry=1
 
255
 
 
256
      softtest=.false.
 
257
      colltest=.false.
 
258
 
 
259
      do jj=1,ndim
 
260
         x(jj)=ran2()
 
261
      enddo
 
262
      call generate_momenta(ndim,iconfig,wgt,x,p)
 
263
      calculatedBorn=.false.
 
264
      do while (( wgt.lt.0 .or. p(0,1).le.0d0 .or. p_born(0,1).le.0d0
 
265
     &           ) .and. ntry .lt. 1000)
 
266
         do jj=1,ndim
 
267
            x(jj)=ran2()
 
268
         enddo
 
269
         wgt=1d0
 
270
         call generate_momenta(ndim,iconfig,wgt,x,p)
 
271
         calculatedBorn=.false.
 
272
         ntry=ntry+1
 
273
      enddo
 
274
 
 
275
      if (ntry.ge.1000) then
 
276
         write (*,*) 'No points passed cuts...'
 
277
         write (12,*) 'ERROR: no points passed cuts...'
 
278
     &        //' Cannot perform ME tests properly for config',iconfig
 
279
         exit
 
280
      endif
 
281
 
 
282
      call sborn(p_born,wgt1)
 
283
      
 
284
      write (*,*) ''
 
285
      write (*,*) ''
 
286
      write (*,*) ''
 
287
 
 
288
      softtest=.true.
 
289
      colltest=.false.
 
290
      nerr=0
 
291
      imax=10
 
292
      do j=1,nsofttests
 
293
      call get_helicity(i_fks,j_fks)
 
294
 
 
295
         if(nsofttests.le.10)then
 
296
           write (*,*) ' '
 
297
           write (*,*) ' '
 
298
         endif
 
299
 
 
300
         y_ij_fks_fix=y_ij_fks_fix_save
 
301
         xi_i_fks_fix=0.1d0
 
302
         ntry=1
 
303
         wgt=1d0
 
304
         do jj=1,ndim
 
305
            x(jj)=ran2()
 
306
         enddo
 
307
         call generate_momenta(ndim,iconfig,wgt,x,p)
 
308
         do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
 
309
            wgt=1d0
 
310
            do jj=1,ndim
 
311
               x(jj)=ran2()
 
312
            enddo
 
313
            call generate_momenta(ndim,iconfig,wgt,x,p)
 
314
            ntry=ntry+1
 
315
         enddo
 
316
         if(nsofttests.le.10)write (*,*) 'ntry',ntry
 
317
         calculatedBorn=.false.
 
318
         call set_cms_stuff(0)
 
319
         call sreal(p1_cnt(0,1,0),zero,y_ij_fks_ev,fxl) 
 
320
         fxl=fxl*jac_cnt(0)
 
321
 
 
322
         call set_cms_stuff(-100)
 
323
         call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx)
 
324
         limit(1)=fx*wgt
 
325
         wlimit(1)=wgt
 
326
 
 
327
         do k=1,nexternal
 
328
            do l=0,3
 
329
               lxp(l,k)=p1_cnt(l,k,0)
 
330
               xp(1,l,k)=p(l,k)
 
331
            enddo
 
332
         enddo
 
333
         do l=0,3
 
334
            lxp(l,nexternal+1)=p_i_fks_cnt(l,0)
 
335
            xp(1,l,nexternal+1)=p_i_fks_ev(l)
 
336
         enddo
 
337
 
 
338
         do i=2,imax
 
339
            xi_i_fks_fix=xi_i_fks_fix/10d0
 
340
            wgt=1d0
 
341
            call generate_momenta(ndim,iconfig,wgt,x,p)
 
342
            calculatedBorn=.false.
 
343
            call set_cms_stuff(-100)
 
344
            call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx)
 
345
            limit(i)=fx*wgt
 
346
            wlimit(i)=wgt
 
347
            do k=1,nexternal
 
348
               do l=0,3
 
349
                  xp(i,l,k)=p(l,k)
 
350
               enddo
 
351
            enddo
 
352
            do l=0,3
 
353
               xp(i,l,nexternal+1)=p_i_fks_ev(l)
 
354
            enddo
 
355
         enddo
 
356
 
 
357
         if(nsofttests.le.10)then
 
358
           write (*,*) 'Soft limit:'
 
359
           do i=1,imax
 
360
              call xprintout(6,limit(i),fxl)
 
361
           enddo
 
362
c
 
363
           write(80,*)'  '
 
364
           write(80,*)'****************************'
 
365
           write(80,*)'  '
 
366
           do k=1,nexternal+1
 
367
              write(80,*)''
 
368
              write(80,*)'part:',k
 
369
              do l=0,3
 
370
                 write(80,*)'comp:',l
 
371
                 do i=1,10
 
372
                    call xprintout(80,xp(i,l,k),lxp(l,k))
 
373
                 enddo
 
374
              enddo
 
375
           enddo
 
376
        else
 
377
           iflag=0
 
378
           call checkres(limit,fxl,wlimit,jac_cnt(0),xp,lxp,
 
379
     &                   iflag,imax,j,nexternal,i_fks,j_fks,iret)
 
380
           nerr=nerr+iret
 
381
        endif
 
382
 
 
383
      enddo
 
384
      if(nsofttests.gt.10)then
 
385
         write(*,*)'Soft tests done for (Born) config',iconfig
 
386
         write(*,*)'Failures:',nerr
 
387
         fail_frac= nerr/dble(nsofttests)
 
388
         if (fail_frac.lt.max_fail) then
 
389
             write(*,401) nFKSprocess, fail_frac
 
390
         else
 
391
             write(*,402) nFKSprocess, fail_frac
 
392
         endif
 
393
      endif
 
394
 
 
395
      write (*,*) ''
 
396
      write (*,*) ''
 
397
      write (*,*) ''
 
398
      
 
399
      include 'pmass.inc'
 
400
 
 
401
      if (pmass(j_fks).ne.0d0) then
 
402
         write (*,*) 'No collinear test for massive j_fks'
 
403
         goto 123
 
404
      endif
 
405
 
 
406
      softtest=.false.
 
407
      colltest=.true.
 
408
 
 
409
c Set rotategranny=.true. to align grandmother along the z axis, when 
 
410
c grandmother is not the c.m. system (if granny=cms, this rotation coincides
 
411
c with the identity, and the following is harmless).
 
412
c WARNING: the setting of rotategranny changes the definition of xij_aor
 
413
c in genps_fks_test.f
 
414
      rotategranny=.false.
 
415
 
 
416
      nerr=0
 
417
      imax=10
 
418
      do j=1,ncolltests
 
419
         call get_helicity(i_fks,j_fks)
 
420
 
 
421
         if(ncolltests.le.10)then
 
422
            write (*,*) ' '
 
423
            write (*,*) ' '
 
424
         endif
 
425
 
 
426
         y_ij_fks_fix=0.9d0
 
427
         xi_i_fks_fix=xi_i_fks_fix_save
 
428
         ntry=1
 
429
         wgt=1d0
 
430
         do jj=1,ndim
 
431
            x(jj)=ran2()
 
432
         enddo
 
433
         call generate_momenta(ndim,iconfig,wgt,x,p)
 
434
         do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
 
435
            wgt=1d0
 
436
            do jj=1,ndim
 
437
               x(jj)=ran2()
 
438
            enddo
 
439
            call generate_momenta(ndim,iconfig,wgt,x,p)
 
440
            ntry=ntry+1
 
441
         enddo
 
442
         if(ncolltests.le.10)write (*,*) 'ntry',ntry
 
443
         calculatedBorn=.false.
 
444
         call set_cms_stuff(1)
 
445
         call sreal(p1_cnt(0,1,1),xi_i_fks_cnt(1),one,fxl) 
 
446
         fxl=fxl*jac_cnt(1)
 
447
 
 
448
         call set_cms_stuff(-100)
 
449
         call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx)
 
450
         limit(1)=fx*wgt
 
451
         wlimit(1)=wgt
 
452
 
 
453
         do k=1,nexternal
 
454
            do l=0,3
 
455
               lxp(l,k)=p1_cnt(l,k,1)
 
456
               xp(1,l,k)=p(l,k)
 
457
            enddo
 
458
         enddo
 
459
         do l=0,3
 
460
            lxp(l,nexternal+1)=p_i_fks_cnt(l,1)
 
461
            xp(1,l,nexternal+1)=p_i_fks_ev(l)
 
462
         enddo
 
463
 
 
464
         do i=2,imax
 
465
            y_ij_fks_fix=1-0.1d0**i
 
466
            wgt=1d0
 
467
            call generate_momenta(ndim,iconfig,wgt,x,p)
 
468
            calculatedBorn=.false.
 
469
            call set_cms_stuff(-100)
 
470
            call sreal(p,xi_i_fks_ev,y_ij_fks_ev,fx)
 
471
            limit(i)=fx*wgt
 
472
            wlimit(i)=wgt
 
473
            do k=1,nexternal
 
474
               do l=0,3
 
475
                  xp(i,l,k)=p(l,k)
 
476
               enddo
 
477
            enddo
 
478
            do l=0,3
 
479
               xp(i,l,nexternal+1)=p_i_fks_ev(l)
 
480
            enddo
 
481
         enddo
 
482
         if(ncolltests.le.10)then
 
483
            write (*,*) 'Collinear limit:'
 
484
            do i=1,imax
 
485
               call xprintout(6,limit(i),fxl)
 
486
            enddo
 
487
c
 
488
            write(80,*)'  '
 
489
            write(80,*)'****************************'
 
490
            write(80,*)'  '
 
491
            do k=1,nexternal+1
 
492
               write(80,*)''
 
493
               write(80,*)'part:',k
 
494
               do l=0,3
 
495
                  write(80,*)'comp:',l
 
496
                  do i=1,10
 
497
                     call xprintout(80,xp(i,l,k),lxp(l,k))
 
498
                  enddo
 
499
               enddo
 
500
            enddo
 
501
         else
 
502
            iflag=1
 
503
            call checkres(limit,fxl,wlimit,jac_cnt(1),xp,lxp,
 
504
     &                    iflag,imax,j,nexternal,i_fks,j_fks,iret)
 
505
            nerr=nerr+iret
 
506
         endif
 
507
      enddo
 
508
      if(ncolltests.gt.10)then
 
509
         write(*,*)'Collinear tests done for (Born) config', iconfig
 
510
         write(*,*)'Failures:',nerr
 
511
         fail_frac= nerr/dble(ncolltests)
 
512
         if (fail_frac.lt.max_fail) then
 
513
             write(*,501) nFKSprocess, fail_frac
 
514
         else
 
515
             write(*,502) nFKSprocess, fail_frac
 
516
         endif
 
517
      endif
 
518
 
 
519
 123  continue
 
520
 
 
521
      enddo                     ! Loop over Born configurations
 
522
      enddo                     ! Loop over nFKSprocess
 
523
 
 
524
 
 
525
      return
 
526
 401  format('     Soft test ',i2,' PASSED. Fraction of failures: ',
 
527
     & f4.2) 
 
528
 402  format('     Soft test ',I2,' FAILED. Fraction of failures: ',
 
529
     & f4.2) 
 
530
 501  format('Collinear test ',i2,' PASSED. Fraction of failures: ',
 
531
     & f4.2) 
 
532
 502  format('Collinear test ',I2,' FAILED. Fraction of failures: ',
 
533
     & f4.2) 
 
534
      end
 
535
 
 
536
c
 
537
c
 
538
c Dummy routines
 
539
c
 
540
c
 
541
      subroutine clear_events()
 
542
      end
 
543
      subroutine store_events()
 
544
      end
 
545
      integer function n_unwgted()
 
546
      n_unwgted = 1
 
547
      end
 
548
 
 
549
      subroutine outfun(pp,www)
 
550
      implicit none
 
551
      include 'nexternal.inc'
 
552
      real*8 pp(0:3,nexternal),www
 
553
c
 
554
      write(*,*)'This routine should not be called here'
 
555
      stop
 
556
      end