Merge lp:~maddevelopers/mg5amcnlo/version_1_0_1 into lp:~madteam/mg5amcnlo/trunk

[Johan Alwall]

- Removed hard limit on number of external particles in
MadEvent, allowing for unlimited length decay chains there
(up to 14 final state particles successfully integrated)
- Improved helicity selection and automatic full helicity
sum if needed. Optimization of run parameters.
- Fixed MadEvent bugs for multiparton processes with
conflicting decays and fixed some faulty array limits.
- Suppressed scary but irrelevant warnings for compiling
2->1 and 1->2 processes in MadEvent
- Optimized Pythia 8 output further
- Several minor fixes relating to user interface etc.

[lp:~maddevelopers/mg5amcnlo/version_1_0_1]

132. By Johan Alwall on 2011-04-19

Include possibility to choose in run_card whether or not apply cuts on decay products (pt/E/eta/dr/mij)

133. By Johan Alwall on 2011-04-19

Added restrict files setting tau and b masses to 0 for mssm

134. By Olivier Mattelaer on 2011-04-20

      make output of launch to have correct unit when computing width
      correct bug 761743 linked to lha_read.f

135. By Johan Alwall on 2011-04-20

Merged with reduce_aloha_call branch rev. 141: Combine multiple lorentz with the same color in the same vertex into a single aloha call, speeding up things like SM Z production and decay by a factor 2 or more.

136. By Johan Alwall on 2011-04-20

Ensure to flip color info etc for mirror processes after addmothers, since this might change there

137. By Johan Alwall on 2011-04-21

Fixed number of digits for Gxx.xx directories to work for all cases, fixed grid presetting also for forced BW without cuts

138. By Johan Alwall on 2011-04-21

Woops, forgot symmetry_v4.f

139. By Johan Alwall on 2011-04-21

Updated version number and UpdateNotes

[Olivier Mattelaer (olivier-mattelaer)]

Nice Job,

here is the first part of my review (the diff is too long to be on the web:)

1) At a point you use this:

2749+ start_time = time.time()
2750 for perm in itertools.permutations(range(ninitial, nexternal)):
2751 if [equivalent_process.get('legs')[i].get('id') for i in perm] != \
2752 final_states:
2753 # Non-identical particles permutated
2754 continue
2755+ # Check if we need to cancel the check due to time limit
2756+ if not max_time or time.time() - start_time > max_time:
2757+ logger.warning("Cancel diagram symmetry - time exceeded")
2758+ break

I would suggest to rather use a clock alarm.

for example:
def handle_alarm(signum, frame):
    raise TimeOutError

signal.signal(signal.SIGALRM, handle_alarm)
signal.alarm(timeout)
try:
 WHAT YOU HAVE TO DO
except TimeOutError:
 STOPPED IT NICELY
finally:
 signal.alarm(0) # stop the alarm since the job is done

*********************************************************************

2) models/RS/model.pkl need to be suppresed
It seemed to be include in the repository (this should be my fault. sorry for that)

So I'm at 50% of the review, and so far so good.

[Johan Alwall (johan-alwall)]

Hello Olivier,

> I would suggest to rather use a clock alarm.
>
> for example:
> def handle_alarm(signum, frame):
>    raise TimeOutError
>
> signal.signal(signal.SIGALRM, handle_alarm)
> signal.alarm(timeout)
> try:
>        WHAT YOU HAVE TO DO
> except TimeOutError:
>        STOPPED IT NICELY
> finally:
>        signal.alarm(0) # stop the alarm since the job is done

Excellent! This is much, much better! I've made the change.

> 2) models/RS/model.pkl need to be suppresed
> It seemed to be include in the repository (this should be my fault. sorry for that)

Already removed!

Thanks! Looking forward to the other half!

Johan

[lp:~maddevelopers/mg5amcnlo/version_1_0_1]

140. By Johan Alwall on 2011-04-21

Replaced stupid time.time in diagram symmetry with elegant Python alarm

141. By Johan Alwall on 2011-04-21

Cutoff when fact scale too low instead of ungraceful exit. Set minimum for grid even if no cuts (at 1/10000 x sqrt(s)). Changed refresh in crossx.html to 10 s instead of 30 s.

142. By Johan Alwall on 2011-04-21

Added floating point trap to commented line in makefiles

143. By Johan Alwall on 2011-04-21

Merged with updated reduce_aloha_call branch, including open command for relevant file types, and automatic opening of crossx.html when running launch for MadEvent.

144. By Johan Alwall on 2011-04-21

Set all applications directly in configuration step, no selection done in open_file

145. By Johan Alwall on 2011-04-21

Added .ps as valid extention for eps_viewer

146. By Johan Alwall on 2011-04-21

Updated UpdateNotes.

147. By Johan Alwall on 2011-04-21

Finally fixed the way default directories are set. Now, always check for present directory if none is given.

148. By Johan Alwall on 2011-04-22

Finally removed superfluous and confusing check if export has already happened to this directory.

149. By Johan Alwall on 2011-04-22

Remove scary output about no valid eps viewer or web browser

=== modified file 'Template/Cards/run_card.dat'
--- Template/Cards/run_card.dat	2011-03-02 18:01:07 +0000
+++ Template/Cards/run_card.dat	2011-04-22 01:24:30 +0000
@@ -70,6 +70,11 @@
 # BW cutoff (M+/-bwcutoff*Gamma)
 #**********************************
   15  = bwcutoff
+#**********************************************************
+# Apply pt/E/eta/dr/mij cuts on decay products or not
+# (note that etmiss/ptll/ptheavy/sorted cuts always apply)
+#**********************************************************
+   F  = cut_decays
 #*******************                                                 
 # Standard Cuts
 #*******************                                                 

=== modified file 'Template/Source/combine_events.f'
--- Template/Source/combine_events.f	2011-03-01 07:02:57 +0000
+++ Template/Source/combine_events.f	2011-04-22 01:24:30 +0000
@@ -13,15 +13,16 @@
       parameter (cmax_events=500000)
       integer    sfnum
       parameter (sfnum=17)   !Unit number for scratch file
+      include 'maxparticles.inc'
       integer    maxexternal
-      parameter (maxexternal=17)
+      parameter (maxexternal=2*max_particles-3)
       integer maxpara
       parameter (maxpara=1000)
 c     
 c     Local
 c
-      character*100 subname(maxsubprocesses)
-      character*110 pathsubname(maxsubprocesses)         !needed for MadWeight
+      character*300 subname(maxsubprocesses)
+      character*310 pathsubname(maxsubprocesses)         !needed for MadWeight
       character*80 down_path                           !needed for MadWeight
       character*40 filename                         !needed for MadWeight
       character*4  card_number                         !needed for MadWeight
@@ -442,6 +443,7 @@
 c
       character*(*) symfile
       parameter (symfile='symfact.dat')
+      include 'maxparticles.inc'
 c
 c     Arguments
 c
@@ -453,37 +455,31 @@
 c
       integer i,j, k, ip
       double precision xi
-      character*150 dirname,dname,channame
+      character*300 dirname,dname,channame
+      integer ncode,npos
+      character*20 formstr
 c-----
 c  Begin Code
 c-----
-       i = index(dir," ")
+      i = index(dir," ")
+c     ncode is number of digits needed for the bw coding
+      ncode=int(dlog10(3d0)*(max_particles-3))+1
       dname = dir(1:i-1)// "/" // symfile
       open(unit=35, file=dname ,status='old',err=59)
       do while (.true.)
          read(35,*,err=99,end=99) xi,j
          if (j .gt. 0) then
-            if ( (xi-int(xi+.01)) .lt. 1d-5) then
-               k = int(xi+.01)
-               if (k .lt. 10) then
-                  write(dirname,'(a,i1,a)') 'G',k,'/'
-               else if (k .lt. 100) then
-                  write(dirname,'(a,i2,a)') 'G',k,'/'
-               else if (k .lt. 1000) then
-                  write(dirname,'(a,i3,a)') 'G',k,'/'
-               else if (k .lt. 10000) then
-                  write(dirname,'(a,i4,a)') 'G',k,'/'
-               endif
+            k = int(xi*(1+10**-ncode))
+            npos=int(dlog10(dble(k)))+1
+            if ( (xi-k) .eq. 0) then
+c              Write with correct number of digits
+               write(formstr,'(a,i1,a)') '(a,i',npos,',a)'
+               write(dirname, formstr) 'G',k,'/'
             else               !Handle B.W.
-               if (xi .lt. 10) then
-                  write(dirname,'(a,f5.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 100) then
-                  write(dirname,'(a,f6.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 1000) then
-                  write(dirname,'(a,f7.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 10000) then
-                  write(dirname,'(a,f8.3,a,a)') 'G',xi,'/'
-               endif
+c              Write with correct number of digits
+               write(formstr,'(a,i1,a,i1,a)') '(a,f',npos+ncode+1,
+     $                 '.',ncode,',a)'
+               write(dirname,formstr)  'G',xi,'/'
             endif     
             ip = index(dirname,'/')
             channame = dname(1:i-1)// "/" //dirname(1:ip)
@@ -517,8 +513,9 @@
       parameter (sfnum=17)   !Unit number for scratch file
       character*(*) scaled_file
       parameter (scaled_file='events.lhe')
-      integer maxexternal
-      parameter (maxexternal=17)
+      include 'maxparticles.inc'
+      integer    maxexternal
+      parameter (maxexternal=2*max_particles-3)
       include 'run_config.inc'
       integer    max_read
       parameter (max_read = 2000000)
@@ -540,7 +537,7 @@
       integer ievent,iseed
       logical done,found
       character*140 buff
-      character*150 fullname
+      character*300 fullname
 c
 c     Les Houches init block (for the <init> info)
 c
@@ -654,7 +651,7 @@
 c
 c     Arguments
 c
-      character*100 subname(maxsubprocesses)
+      character*300 subname(maxsubprocesses)
       integer ns
 c-----
 c  Begin Code

=== modified file 'Template/Source/combine_runs.f'
--- Template/Source/combine_runs.f	2010-12-20 20:08:13 +0000
+++ Template/Source/combine_runs.f	2011-04-22 01:24:30 +0000
@@ -14,9 +14,9 @@
 c     
 c     Local
 c
-      character*150 procname(maxsubprocesses)
-      character*150 channame(maxsubprocesses)
-      character*150 pathname
+      character*300 procname(maxsubprocesses)
+      character*300 channame(maxsubprocesses)
+      character*300 pathname
       integer iproc,nproc
       integer ichan,nchan
       integer jc
@@ -88,6 +88,7 @@
 c
       character*(*) symfile
       parameter    (symfile='symfact.dat')
+      include 'maxparticles.inc'
       integer    maxsubprocesses
       parameter (maxsubprocesses=999)
 c
@@ -99,42 +100,36 @@
 c
 c     Local
 c
-      character*150 fname
+      character*300 fname
       character*50 dirname
       integer jc,ip
       double precision xi
       integer j,k
+      integer ncode,npos
+      character*20 formstr
 c-----
 c  Begin Code
 c-----
       jc = index(pathname," ")
+c     ncode is number of digits needed for the bw coding
+      ncode=int(dlog10(3d0)*(max_particles-3))+1
       fname = pathname(1:jc-1) // "/" // symfile
       nchan = 0
       open(unit=35, file=fname,status='old',err=99)
       do while (.true.)
          read(35,*,err=99,end=99) xi,j
          if (j .gt. 0) then
-            if ( (xi-int(xi+.01)) .lt. 1d-5) then
-               k = int(xi+.01)
-               if (k .lt. 10) then
-                  write(dirname,'(a,i1,a)') 'G',k,'/'
-               else if (k .lt. 100) then
-                  write(dirname,'(a,i2,a)') 'G',k,'/'
-               else if (k .lt. 1000) then
-                  write(dirname,'(a,i3,a)') 'G',k,'/'
-               else if (k .lt. 10000) then
-                  write(dirname,'(a,i4,a)') 'G',k,'/'
-               endif
+            k = int(xi*(1+10**-ncode))
+            npos=int(dlog10(dble(k)))+1
+            if ( (xi-k) .eq. 0) then
+c              Write with correct number of digits
+               write(formstr,'(a,i1,a)') '(a,i',npos,',a)'
+               write(dirname, formstr) 'G',k,'/'
             else               !Handle B.W.
-               if (xi .lt. 10) then
-                  write(dirname,'(a,f5.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 100) then
-                  write(dirname,'(a,f6.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 1000) then
-                  write(dirname,'(a,f7.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 10000) then
-                  write(dirname,'(a,f8.3,a,a)') 'G',xi,'/'
-               endif
+c              Write with correct number of digits
+               write(formstr,'(a,i1,a,i1,a)') '(a,f',npos+ncode+1,
+     $                 '.',ncode,',a)'
+               write(dirname,formstr)  'G',xi,'/'
             endif     
             ip = index(dirname,'/')
             nchan=nchan+1            
@@ -165,8 +160,8 @@
 c
 c     Local
 c
-      character*150 fname
-      character*150 pathname
+      character*300 fname
+      character*300 pathname
       character*26 alphabet
       integer jc,i
       integer njobs
@@ -263,8 +258,9 @@
 c
       character*(*) eventfname
       parameter    (eventfname="events.lhe")
+      include 'maxparticles.inc'
       integer    maxexternal     !Max number external momenta
-      parameter (maxexternal=17)
+      parameter (maxexternal=2*max_particles-3)
 c     
 c     Arguments
 c
@@ -274,9 +270,9 @@
 c     
 c     Local
 c
-      character*150 fname
+      character*300 fname
       integer jc
-      double precision P(0:3,maxexternal),xwgt
+      double precision P(0:4,maxexternal),xwgt
       integer n,ic(7,maxexternal),ievent
       double precision scale,aqcd,aqed
       character*140 buff
@@ -328,7 +324,7 @@
 c     
 c     Local
 c
-      character*150 fname
+      character*300 fname
       integer jc,i,j
       double precision x1,x2,x3,x4
 c-----
@@ -385,7 +381,7 @@
 c     
 c     Local
 c
-      character*150 fname
+      character*300 fname
       character*26 alphabet
       integer jc,i,j
 c-----
@@ -429,7 +425,7 @@
 c     
 c     Local
 c
-      character*150 fname
+      character*300 fname
       integer jc,i,j
       double precision x1,x2
       integer i3,i4

=== modified file 'Template/Source/cuts.inc'
--- Template/Source/cuts.inc	2010-12-19 17:10:07 +0000
+++ Template/Source/cuts.inc	2011-04-22 01:24:30 +0000
@@ -65,3 +65,9 @@
       real*8 xptj,xptb,xpta,xptl,xmtc
       real*8 xetamin,xqcut,deltaeta
       common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
+
+C
+C    FLAG FOR CUTTING DECAY PRODUCTS
+C
+      logical cut_decays
+      common /to_cut_decays/cut_decays

=== modified file 'Template/Source/dsample.f'
--- Template/Source/dsample.f	2011-04-04 06:18:33 +0000
+++ Template/Source/dsample.f	2011-04-22 01:24:30 +0000
@@ -1755,7 +1755,7 @@
 c
 c              special routines to deal with xmin cutoff
 c
-               do while(grid(1,k,j) .le. 0d0 .and. k .lt. ng)
+               do while(grid(1,k,j) .le. 0d0 .and. k+1 .lt. ng)
                   k=k+1
                enddo
 
@@ -1940,7 +1940,6 @@
                tmeant = 0d0
                tsigmat = 0d0
                do i=cur_it-3,cur_it-1
-                  print *,'mean, sigma: ',i,ymean(i),ysigma(i)
                   tmeant = tmeant+ymean(i)*ymean(i)**2/ysigma(i)**2
                   tsigmat = tsigmat + ymean(i)**2/ ysigma(i)**2
                enddo

=== modified file 'Template/Source/gen_ximprove.f'
--- Template/Source/gen_ximprove.f	2011-03-15 23:28:24 +0000
+++ Template/Source/gen_ximprove.f	2011-04-22 01:24:30 +0000
@@ -17,6 +17,7 @@
       parameter (maxpara=1000)
 
       include 'run_config.inc'
+      include 'maxparticles.inc'
 c
 c     global
 c
@@ -46,6 +47,8 @@
       logical Gridpack,gridrun
       logical split_channels
       common /to_split/split_channels
+      integer ncode,npos
+      character*20 formstr
 
 c-----
 c  Begin Code
@@ -106,33 +109,25 @@
       errtotc=0d0
       xtot = 0d0
       i = 0
+c     ncode is number of digits needed for the bw coding
+      ncode=int(dlog10(3d0)*(max_particles-3))+1
       do while (.true.)
          read(15,*,err=99,end=99) xi,j
          if (j .gt. 0) then
             i = i+1
-            if ( (xi-int(xi+.01)) .lt. 1d-5) then
-               n = int(xi+.01)
-               if (n .lt. 10) then
-                  write(fname,'(a,i1,a,a)') 'G',n,'/',rfile
-               else if (n .lt. 100) then
-                  write(fname,'(a,i2,a,a)') 'G',n,'/',rfile
-               else if (n .lt. 1000) then
-                  write(fname,'(a,i3,a,a)') 'G',n,'/',rfile
-               else if (n .lt. 10000) then
-                  write(fname,'(a,i4,a,a)') 'G',n,'/',rfile
-               endif
-            else
-               if (xi .lt. 10) then
-                  write(fname,'(a,f5.3,a,a)') 'G',xi,'/',rfile
-               else if (xi .lt. 100) then
-                  write(fname,'(a,f6.3,a,a)') 'G',xi,'/',rfile
-               else if (xi .lt. 1000) then
-                  write(fname,'(a,f7.3,a,a)') 'G',xi,'/',rfile
-               else if (xi .lt. 10000) then
-                  write(fname,'(a,f8.3,a,a)') 'G',xi,'/',rfile
-               endif
+            k = int(xi*(1+10**-ncode))
+            npos=int(dlog10(dble(k)))+1
+            if ( (xi-k) .eq. 0) then
+c              Write with correct number of digits
+               write(formstr,'(a,i1,a)') '(a,i',npos,',a,a)'
+               write(fname, formstr) 'G',k,'/',rfile
+            else               !Handle B.W.
+c              Write with correct number of digits
+               write(formstr,'(a,i1,a,i1,a)') '(a,f',npos+ncode+1,
+     $                 '.',ncode,',a,a)'
+               write(fname, formstr) 'G',xi,'/',rfile
+            endif
 c              write(*,*) 'log name ',fname
-            endif
          endif
          if (j .gt. 0) then
             gname(i)=fname
@@ -314,17 +309,6 @@
 
          ip = index(gn(io(i)),'/')
          write(26,'(2a)') 'j=',gn(io(i))(1:ip-1)
-         if (.false.) then
-         if (io(i) .lt. 10) then
-            write(26,'(a,i1)') 'j=G',io(i)
-         elseif (io(i) .lt. 100) then
-            write(26,'(a,i2)') 'j=G',io(i)
-         elseif (io(i) .lt. 1000) then
-            write(26,'(a,i3)') 'j=G',io(i)
-         elseif (io(i) .lt. 10000) then
-            write(26,'(a,i4)') 'j=G',io(i)
-         endif
-         endif
 c
 c     Determine estimates for getting the desired accuracy
 c
@@ -386,7 +370,8 @@
 c     local
 c
       character*30 fname
-      integer ic
+      integer ic, npos
+      character*10 formstr
 
       data ic/0/
 c-----
@@ -394,20 +379,12 @@
 c-----
       ic=ic+1
       fname='ajob'
-      if (ic .lt. 10) then
-         write(fname(5:5),'(i1)') ic
-      elseif (ic .lt. 100) then
-         write(fname(5:6),'(i2)') ic
-      elseif (ic .lt. 1000) then
-         write(fname(5:7),'(i3)') ic
-      elseif (ic .lt. 10000) then
-         write(fname(5:8),'(i4)') ic
-      else
-         write(*,*) 'Error too many jobs in gen_ximprove',ic
-         stop
-      endif
+c     Write ic with correct number of digits
+      npos=int(dlog10(dble(ic)))+1
+      write(formstr,'(a,i1,a)') '(I',npos,')'
+      write(fname(5:(5+npos-1)),formstr) ic
+
       write(*,*) 'Opening file ',fname
-
       open (unit=26, file = fname, status='unknown')
       write(26,15) '#!/bin/bash'
       write(26,15) '#PBS -q ' // PBS_QUE
@@ -559,9 +536,10 @@
          k=k+1
       enddo
       kk=k
-c     Check error for the rest of the channels
+c     Check error for the rest of the channels - rerun if 
+c     bigger than channel xsec and bigger than 1% of largest channel
       do while( kk .le. ng)
-         if (xerr(io(kk)).gt.xsec(io(kk))) then
+         if (xerr(io(kk)).gt.max(xsec(io(kk)),0.01*xsec(io(1)))) then
             write(*,*) 'Improving for error ',kk,gn(io(kk)),xt(kk),xsec(io(kk)),xerr(io(kk))
             io(k)=io(kk)
             xt(k)=xt(kk)

=== modified file 'Template/Source/genps.inc'
--- Template/Source/genps.inc	2011-04-03 22:06:00 +0000
+++ Template/Source/genps.inc	2011-04-22 01:24:30 +0000
@@ -2,12 +2,11 @@
 c       Parameters used by genps and dsample, you must recompile
 c	dsample if you change anything below
 c*************************************************************************
-      integer    max_particles  , max_branch
-      parameter (max_particles=10, max_branch=max_particles-2)
-      integer    max_tree    , max_dim
-      parameter (max_tree=200, max_dim=max_tree*(max_branch+1))
+      include 'maxparticles.inc'
+      integer    max_branch
+      parameter (max_branch=max_particles-1)
       integer    ng     , maxdim     , maxinvar     , maxconfigs
-      parameter (ng = 96, maxdim = 25, maxinvar= 50 , maxconfigs=10)
+      parameter (ng = 96, maxdim = 3*(max_particles-2)-2, maxinvar= 4*max_particles, maxconfigs=10)
 c      parameter (ng = 43, maxdim = 25, maxinvar= 57 , maxconfigs=20)
 
       double precision  xgmin,     xgmax

=== modified file 'Template/Source/invarients.f'
--- Template/Source/invarients.f	2011-04-03 22:06:00 +0000
+++ Template/Source/invarients.f	2011-04-22 01:24:30 +0000
@@ -226,7 +226,8 @@
 c
 
       if (nconfigs .eq. 1) then
-         do j=1,3*nbranch-4+2
+c         do j=1,3*nbranch-4+2
+         do j=1,maxdim
             minvar(j,mincfig)=j
          enddo
          jgrid=j-1

=== modified file 'Template/Source/makefile'
--- Template/Source/makefile	2010-12-19 04:12:16 +0000
+++ Template/Source/makefile	2011-04-22 01:24:30 +0000
@@ -1,6 +1,6 @@
 g77   = g77
 FFLAGS= -O -ffixed-line-length-132
-#FFLAGS= -O 
+#FFLAGS= -g -ffixed-line-length-132 -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal
 LIBDIR= ../lib/
 BINDIR= ../bin/
 PDFDIR= ./PDF/

=== added file 'Template/Source/maxparticles.inc'
--- Template/Source/maxparticles.inc	1970-01-01 00:00:00 +0000
+++ Template/Source/maxparticles.inc	2011-04-22 01:24:30 +0000
@@ -0,0 +1,2 @@
+      integer    max_particles
+      parameter (max_particles=10)

=== modified file 'Template/Source/run_config.inc'
--- Template/Source/run_config.inc	2010-12-19 04:12:16 +0000
+++ Template/Source/run_config.inc	2011-04-22 01:24:30 +0000
@@ -9,18 +9,18 @@
       integer    icomp
       parameter (icomp = 3)   !BW + Symmetry compression  0 == none	
       integer    min_events_subprocess          !Minimum number of events
-      parameter (min_events_subprocess = 2000) !per iteration in each subprocess
+      parameter (min_events_subprocess = 1000) !per iteration in each subprocess
       integer    min_events_channel         !Minimum number of events
       parameter (min_events_channel = 1000) !per iteration in each channel
       integer    iter_survey                !Number of iterations for survey
-      parameter (iter_survey=4)
+      parameter (iter_survey=3)
       integer    nhel_survey                !Number of helicities to sum for each
       parameter (nhel_survey=1)             !phase space point. (0 = sum_all)
 c
 c     The following parameters are used by gen_ximprove.f in running refine
 c
       integer    min_events                 !Minimum number of events/iteration
-      parameter (min_events = 2000)         !to refine a channel
+      parameter (min_events = 1000)         !to refine a channel
       integer    max_events                 !Maximum number of events/iteration
       parameter (max_events = 2000)        !to refine a channel
       integer    max_iter                   !Maximum number of iterations
@@ -51,6 +51,6 @@
       double precision acc_wu
       parameter       (acc_wu = 0.01)   !Desired accuracy from warmup run
       integer      npoints_wu,        itmax_wu !warmup # points/iterations
-      parameter   (npoints_wu = 4000, itmax_wu = 8)
+      parameter   (npoints_wu = 2000, itmax_wu = 8)
       integer      min_gevents_wu
       parameter   (min_gevents_wu=1000)   !Minumum # unweighted events to generate from channel

=== modified file 'Template/Source/setrun.f'
--- Template/Source/setrun.f	2010-12-21 00:56:27 +0000
+++ Template/Source/setrun.f	2011-04-22 01:24:30 +0000
@@ -87,6 +87,12 @@
       call  get_integer (npara,param,value,"maxjetflavor",maxjetflavor,4)
 
 c*********************************************************************
+c cut on decay products (pt/e/eta/dr/mij) or not
+c*********************************************************************
+
+      call  get_logical (npara,param,value,"cut_decays",cut_decays,.true.)
+
+c*********************************************************************
 c     Minimum pt's                                                   *
 c*********************************************************************
 

=== modified file 'Template/Source/sum_html.f'
--- Template/Source/sum_html.f	2011-03-23 21:27:45 +0000
+++ Template/Source/sum_html.f	2011-04-22 01:24:30 +0000
@@ -7,6 +7,7 @@
 c
 c     Constants
 c
+      include 'maxparticles.inc'
       character*(*) subfile
       parameter (subfile='subproc.mg')
       character*(*) symfile
@@ -14,7 +15,7 @@
       character*(*) rfile
       parameter (rfile='results.dat')
       integer   max_amps      , max_iter
-      parameter (max_amps=9999, max_iter=50)
+      parameter (max_amps=9999, max_iter=50)      
 c
 c     local
 c
@@ -34,8 +35,8 @@
       integer icor(max_amps)
       integer nunwgt(max_amps)
       integer minit
-      character*100 fname, pname
-      character*120 linkname(max_amps)
+      character*300 fname, pname
+      character*320 linkname(max_amps)
       integer i,j,k,l
       double precision xtot,errtot,err_goal, xi
       double precision errtotc, errtotu
@@ -44,6 +45,8 @@
       integer ntevents, ntw 
       integer ilen
       logical errex
+      integer ncode,npos
+      character*20 formstr
 
       logical sumproc
       common/to_sumproc/sumproc
@@ -70,6 +73,8 @@
       ntevents = 0
       ntw = 0
       i = 0
+c     ncode is number of digits needed for the bw coding
+      ncode=int(dlog10(3d0)*(max_particles-3))+1
       do while (.true.)
          if (sumproc) then
             j = 1
@@ -84,38 +89,22 @@
             read(15,*,err=99,end=99) xi,j
             if (j .gt. 0) then
             i=i+1
-            if ( (xi-int(xi+.01)) .lt. 1d-5) then
-               k = int(xi+.01)
-            if (k .lt. 10) then
-               write(fname,'(a,i1,a,a)') 'G',k,'/',rfile
-               write(linkname(i),'(a,i1,a,a)') 'G',k,'/','log.txt'
-            else if (k .lt. 100) then
-               write(fname,'(a,i2,a,a)') 'G',k,'/',rfile
-               write(linkname(i),'(a,i2,a,a)') 'G',k,'/','log.txt'
-            else if (k .lt. 1000) then
-               write(fname,'(a,i3,a,a)') 'G',k,'/',rfile
-               write(linkname(i),'(a,i3,a,a)') 'G',k,'/','log.txt'
-            else if (k .lt. 10000) then
-               write(fname,'(a,i4,a,a)') 'G',k,'/',rfile
-               write(linkname(i),'(a,i4,a,a)') 'G',k,'/','log.txt'
-            endif
-            else
-            if (xi .lt. 10) then
-               write(fname,'(a,f5.3,a,a)') 'G',xi,'/',rfile
-               write(linkname(i),'(a,f5.3,a,a)') 'G',xi,'/','log.txt'
-            else if (xi .lt. 100) then
-               write(fname,'(a,f6.3,a,a)') 'G',xi,'/',rfile
-               write(linkname(i),'(a,f6.3,a,a)') 'G',xi,'/','log.txt'
-            else if (xi .lt. 1000) then
-               write(fname,'(a,f7.3,a,a)') 'G',xi,'/',rfile
-               write(linkname(i),'(a,f7.3,a,a)') 'G',xi,'/','log.txt'
-            else if (xi .lt. 10000) then
-               write(fname,'(a,f8.3,a,a)') 'G',xi,'/',rfile
-               write(linkname(i),'(a,f8.3,a,a)') 'G',xi,'/','log.txt'
+            k = int(xi*(1+10**-ncode))
+            npos=int(dlog10(dble(k)))+1
+            if ( (xi-k) .eq. 0) then
+c              Write with correct number of digits
+               write(formstr,'(a,i1,a)') '(a,i',npos,',a,a)'
+               write(fname, formstr) 'G',k,'/',rfile
+               write(linkname(i),formstr) 'G',k,'/','log.txt'
+            else               !Handle B.W.
+c              Write with correct number of digits
+               write(formstr,'(a,i1,a,i1,a)') '(a,f',npos+ncode+1,
+     $                 '.',ncode,',a,a)'
+               write(fname, formstr) 'G',xi,'/',rfile
+               write(linkname(i),formstr) 'G',xi,'/','log.txt'
             endif
 c            write(*,*) 'log name ',fname
             endif
-            endif
          endif
 
 
@@ -139,7 +128,7 @@
             ntw=ntw+nw(i)
 c            maxit = min(maxit,2)
             if (sumproc) then
-               write(*,'(a20,e15.5)') pname(2:ilen), xsec(i)
+               write(*,'(a,e15.5)') pname(2:ilen), xsec(i)
             else
                write(*,*) fname,i,xsec(i),mfact(i)
             endif
@@ -320,7 +309,7 @@
       integer nsubproc          !Number of specific processes requested
       logical found
       integer ig
-      character*120 linkname(max_amps)
+      character*320 linkname(max_amps)
       integer sname(256)
       integer gname
 c
@@ -328,7 +317,7 @@
 c
       integer i,j,k, io(max_amps), ik
       integer ntot, ip,jp
-      character*100 procname
+      character*300 procname
       character*4 cpref
       character*20 fnamel, fnamee
       double precision scale,xt(max_amps), teff
@@ -445,7 +434,7 @@
       if (sumproc) then
          nsubproc=0
          do i=1,ng
-            procname = linkname(io(i))(:100)
+            procname = linkname(io(i))(:300)
             gname=0
             read(procname(2:index(procname,'_')-1),*,err=20) gname
  20         found = .false.
@@ -488,22 +477,6 @@
 c     Create directory names using the linkname
 c
 c
-            if (.false.) then
-            if (io(i) .lt. 10) then
-               write(fnamel,'(a,i1,a,a)') 'G',io(i),'/',logfile
-               write(fnamee,'(a,i1,a,a)') 'G',io(i),'/',eventfile
-            else if (io(i) .lt. 100) then
-               write(fnamel,'(a,i2,a,a)') 'G',io(i),'/',logfile
-               write(fnamee,'(a,i2,a,a)') 'G',io(i),'/',eventfile
-            else if (io(i) .lt. 1000) then
-               write(fnamel,'(a,i3,a,a)') 'G',io(i),'/',logfile
-               write(fnamee,'(a,i3,a,a)') 'G',io(i),'/',eventfile
-            else if (io(i) .lt. 10000) then
-               write(fnamel,'(a,i4,a,a)') 'G',io(i),'/',logfile
-               write(fnamee,'(a,i4,a,a)') 'G',io(i),'/',eventfile
-            endif
-            endif
-
             ik = index(linkname(io(i)),'log.txt')-1
             fnamel = linkname(io(i))(1:ik) // logfile
             fnamee = linkname(io(i))(1:ik) // eventfile
@@ -523,7 +496,7 @@
      $           xlum(io(i))/scale,'</td></tr>'
 c            write(*,*) io(i),xmax(io(i))
             else
-                procname = linkname(io(i))(:100)
+                procname = linkname(io(i))(:300)
 cxx   tjs 3-20-2006  + cfax 12.05.2006
 c                ip = index(procname,'P')+2 !Strip off first P_
                 ip = index(procname,'P')+1 !Strip off first P

=== modified file 'Template/SubProcesses/cluster.f'
--- Template/SubProcesses/cluster.f	2011-04-03 22:06:00 +0000
+++ Template/SubProcesses/cluster.f	2011-04-22 01:24:30 +0000
@@ -376,8 +376,8 @@
 
       logical             OnBW(-nexternal:0)     !Set if event is on B.W.
       common/to_BWEvents/ OnBW
-      integer           mincfig, maxcfig
-      common/to_configs/mincfig, maxcfig
+      integer            mapconfig(0:lmaxconfigs), this_config
+      common/to_mconfigs/mapconfig, this_config
 C $B$ IFOREST $B$ !this is a tag for MadWeight
       integer i, iforest(2,-max_branch:-1,lmaxconfigs)
       common/to_forest/ iforest
@@ -405,8 +405,8 @@
       enddo
       ibw=0
       do i=-1,-(nexternal-3),-1
-        icl(i)=icl(iforest(1,i,mincfig))+
-     $     icl(iforest(2,i,mincfig))
+        icl(i)=icl(iforest(1,i,this_config))+
+     $     icl(iforest(2,i,this_config))
         isbw(icl(i))=.false.
 C $B$ ONBW $B$ !this is a tag for MadWeight
         if(OnBW(i))then
@@ -414,7 +414,7 @@
           ibw=ibw+1
           ibwlist(ibw)=icl(i)
           isbw(icl(i))=.true.
-c          print *,'Added BW for resonance ',i,icl(i),mincfig
+c          print *,'Added BW for resonance ',i,icl(i),this_config
           if(ibw.eq.nbw) return
         endif
       enddo

=== modified file 'Template/SubProcesses/cuts.f'
--- Template/SubProcesses/cuts.f	2011-04-13 01:05:26 +0000
+++ Template/SubProcesses/cuts.f	2011-04-22 01:24:30 +0000
@@ -54,6 +54,7 @@
       integer i,j,njets,hardj1,hardj2
       REAL*8 XVAR,ptmax1,ptmax2,htj,tmp
       real*8 ptemp(0:3)
+      character*20 formstr
 C
 C     PARAMETERS
 C
@@ -135,12 +136,16 @@
          call initcluster()
 c
 c
-         write(*,'(a10,10i8)') 'Particle',(i,i=nincoming+1,nexternal)
-         write(*,'(a10,10f8.1)') 'Et >',(etmin(i),i=nincoming+1,nexternal)
-         write(*,'(a10,10f8.1)') 'E >',(emin(i),i=nincoming+1,nexternal)
-         write(*,'(a10,10f8.1)') 'Eta <',(etamax(i),i=nincoming+1,nexternal)
+         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'i8)'
+         write(*,formstr) 'Particle',(i,i=nincoming+1,nexternal)
+         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'f8.1)'
+         write(*,formstr) 'Et >',(etmin(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'E >',(emin(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'Eta <',(etamax(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
+         write(formstr,'(a,i2.2,a)')'(a,i2,a,',nexternal,'f8.1)'
          do j=nincoming+1,nexternal-1
-            write(*,'(a,i2,a,10f8.1)') 'd R #',j,'  >',(-0.0,i=nincoming+1,j),
+            write(*,formstr) 'd R #',j,'  >',(-0.0,i=nincoming+1,j),
      &           (r2min(i,j),i=j+1,nexternal)
             do i=j+1,nexternal
                r2min(i,j)=r2min(i,j)*dabs(r2min(i,j))    !Since r2 returns distance squared
@@ -148,12 +153,11 @@
             enddo
          enddo
          do j=nincoming+1,nexternal-1
-            write(*,'(a,i2,a,10f8.1)') 's min #',j,'>',
+            write(*,formstr) 's min #',j,'>',
      &           (s_min(i,j),i=nincoming+1,nexternal)
          enddo
-         write(*,'(a10,10f8.1)') 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
          do j=nincoming+1,nexternal-1
-            write(*,'(a,i2,a,10f8.1)') 'xqcutij #',j,'>',
+            write(*,formstr) 'xqcutij #',j,'>',
      &           (xqcutij(i,j),i=nincoming+1,nexternal)
          enddo
 

=== modified file 'Template/SubProcesses/driver.f'
--- Template/SubProcesses/driver.f	2011-04-03 22:06:00 +0000
+++ Template/SubProcesses/driver.f	2011-04-22 01:24:30 +0000
@@ -31,8 +31,6 @@
       common /sample_status/result_file,where_file,nsteps
       integer           Minvar(maxdim,lmaxconfigs)
       common /to_invar/ Minvar
-      real*8          dsigtot(10)
-      common/to_dsig/ dsigtot
       integer ngroup
       common/to_group/ngroup
       data ngroup/0/
@@ -177,14 +175,15 @@
 c     Constants
 c
       include 'nexternal.inc'
+      include 'maxparticles.inc'
 c
 c     Arguments
 c
-      integer ncall,itmax,iconfig, jconfig
+      integer ncall,itmax,iconfig
 c
 c     Local
 c
-      integer i, j
+      integer i, j, jconfig, ncode
       double precision dconfig
 c
 c     Global
@@ -239,7 +238,9 @@
 
       write(*,10) 'Enter Configuration Number: '
       read(*,*) dconfig
-      iconfig = int(dconfig)
+c     ncode is number of digits needed for the BW code
+      ncode=int(dlog10(3d0)*(max_particles-3))+1
+      iconfig = int(dconfig*(1+10**-ncode))
       write(*,12) 'Running Configuration Number: ',iconfig
 c
 c     Here I want to set up with B.W. we map and which we don't
@@ -250,7 +251,7 @@
          lbw(0)=0
       else
          lbw(0)=1
-         jconfig=dconfig*1000.1
+         jconfig=dconfig*(10**ncode + 0.1)
          write(*,*) 'Using dconfig=',jconfig
          call DeCode(jconfig,lbw(1),3,nexternal)
          write(*,*) 'BW Setting ', (lbw(j),j=1,nexternal-2)

=== modified file 'Template/SubProcesses/makefile'
--- Template/SubProcesses/makefile	2010-11-25 00:43:00 +0000
+++ Template/SubProcesses/makefile	2011-04-22 01:24:30 +0000
@@ -1,6 +1,6 @@
 g77   = g77
-FFLAGS= -O -ffixed-line-length-132
-#FFLAGS = -g -ffixed-line-length-132
+FFLAGS= -O -w -ffixed-line-length-132
+#FFLAGS = -g -ffixed-line-length-132 -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal
 LIBDIR = ../../lib/
 BINDIR = ../../bin/
 # $B$ obj_for_MW $B$ !this is a tag for MadWeight. Don't edit this line

=== modified file 'Template/SubProcesses/myamp.f'
--- Template/SubProcesses/myamp.f	2011-04-03 22:06:00 +0000
+++ Template/SubProcesses/myamp.f	2011-04-22 01:24:30 +0000
@@ -215,32 +215,38 @@
               idenpart=0
               do j=1,2
                 ida(j)=iforest(j,i,iconfig)
-                if((ida(j).lt.0.and.sprop(i,iconfig).eq.sprop(ida(j),iconfig))
-     $             .or.(ida(j).gt.0.and.sprop(i,iconfig).eq.IDUP(ida(j),1,1)))
-     $             idenpart=ida(j)
+                if(ida(j).lt.0) then
+                   if(sprop(i,iconfig).eq.sprop(ida(j),iconfig))
+     $                  idenpart=ida(j)
+                elseif (ida(j).gt.0) then
+                   if(sprop(i,iconfig).eq.IDUP(ida(j),1,1))
+     $                  idenpart=ida(j)
+                endif
               enddo
 c           Always remove if daughter final-state
               if(idenpart.gt.0) then
                 OnBW(i)=.false.
 c           Else remove if daughter forced to be onshell
-              elseif(idenpart.lt.0.and.gForceBW(idenpart, iconfig)) then
-                OnBW(i)=.false.
+              elseif(idenpart.lt.0)then
+                 if(gForceBW(idenpart, iconfig)) then
+                    OnBW(i)=.false.
 c           Else remove daughter if forced to be onshell
-              elseif(idenpart.lt.0.and.gForceBW(i, iconfig)) then
-                OnBW(idenpart)=.false.
+                 elseif(gForceBW(i, iconfig)) then
+                    OnBW(idenpart)=.false.
 c           Else remove either this resonance or daughter, which is closer to mass shell
-              elseif(idenpart.lt.0.and.abs(xmass-pmass(i,iconfig)).gt.
-     $             abs(sqrt(dot(xp(0,idenpart),xp(0,idenpart)))-
-     $             pmass(i,iconfig))) then
-                OnBW(i)=.false.
+                 elseif(abs(xmass-pmass(i,iconfig)).gt.
+     $                   abs(sqrt(dot(xp(0,idenpart),xp(0,idenpart)))-
+     $                   pmass(i,iconfig))) then
+                    OnBW(i)=.false.
 c           Else remove OnBW for daughter
-              else if(idenpart.lt.0) then
-                OnBW(idenpart)=.false.
+                 else
+                    OnBW(idenpart)=.false.
+                 endif
               endif
             endif
             if (onshell .and. (lbw(nbw).eq. 2) ) cut_bw=.true.
             if (.not. onshell .and. (lbw(nbw).eq. 1)) cut_bw=.true.
-c            write(*,*) nbw,xmass,onshell,lbw(nbw),cut_bw
+c            write(*,*) 'cut_bw: ',nbw,xmass,onshell,lbw(nbw),cut_bw
          endif
 
       enddo
@@ -286,6 +292,9 @@
 
       logical gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
       include 'decayBW.inc'
+
+      double precision forced_mass
+      data forced_mass/0d0/
 c
 c     Global
 c
@@ -318,7 +327,7 @@
 
       integer        lbw(0:nexternal)  !Use of B.W.
       common /to_BW/ lbw
-      
+
       include 'coupl.inc'
 
 c
@@ -351,7 +360,7 @@
       tsgn    = +1d0
 c     Reset variables
       nbw = 0
-      do i=1,nexternal
+      do i=1,nexternal-2
          iden_part(-i)=0
          spole(i)=0
          swidth(i)=0
@@ -359,22 +368,28 @@
       do i=-1,-(nexternal-3),-1              !Find all the propagotors
 c     JA 3/31/11 Keep track of identical particles (i.e., radiation vertices)
 c     by tracing the particle identity from the external particle.
-         if(iforest(1,i,iconfig).gt.0.and.
-     $        sprop(i,iconfig).eq.idup(iforest(1,i,iconfig),1,1).or.
-     $      iforest(2,i,iconfig).gt.0.and.
-     $        sprop(i,iconfig).eq.idup(iforest(2,i,iconfig),1,1).or.
-     $      iforest(1,i,iconfig).lt.0.and.
-     $        (iden_part(iforest(1,i,iconfig)).ne.0.and.
+         if(iforest(1,i,iconfig).gt.0) then
+             if (sprop(i,iconfig).eq.idup(iforest(1,i,iconfig),1,1))
+     $        iden_part(i) = sprop(i,iconfig)
+          endif
+         if(iforest(2,i,iconfig).gt.0) then
+            if(sprop(i,iconfig).eq.idup(iforest(2,i,iconfig),1,1))
+     $           iden_part(i) = sprop(i,iconfig)
+         endif
+         if(iforest(1,i,iconfig).lt.0) then
+            if((iden_part(iforest(1,i,iconfig)).ne.0.and.
      $        sprop(i,iconfig).eq.iden_part(iforest(1,i,iconfig)) .or.
      $        gforcebw(iforest(1,i,iconfig),iconfig).and.
-     $        sprop(i,iconfig).eq.sprop(iforest(1,i,iconfig),iconfig)).or.
-     $      iforest(2,i,iconfig).lt.0.and.
-     $        (iden_part(iforest(2,i,iconfig)).ne.0.and.
+     $        sprop(i,iconfig).eq.sprop(iforest(1,i,iconfig),iconfig)))
+     $       iden_part(i) = sprop(i,iconfig)
+         endif
+         if(iforest(2,i,iconfig).lt.0) then
+            if((iden_part(iforest(2,i,iconfig)).ne.0.and.
      $        sprop(i,iconfig).eq.iden_part(iforest(2,i,iconfig)).or.
      $        gforcebw(iforest(2,i,iconfig),iconfig).and.
-     $        sprop(i,iconfig).eq.sprop(iforest(2,i,iconfig),iconfig)))then
-            iden_part(i) = sprop(i,iconfig)
-            endif
+     $        sprop(i,iconfig).eq.sprop(iforest(2,i,iconfig),iconfig)))
+     $           iden_part(i) = sprop(i,iconfig)
+         endif
          if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
          if (tsgn .eq. 1d0) then                         !s channel
             xm(i) = xm(iforest(1,i,iconfig))+xm(iforest(2,i,iconfig))
@@ -418,6 +433,9 @@
                   spole(-i)=pmass(i,iconfig)*pmass(i,iconfig)/stot
                   swidth(-i) = pwidth(i,iconfig)*pmass(i,iconfig)/stot
                   xm(i) = pmass(i,iconfig)
+c                 Remember largest BW mass for better grid setting
+                  forced_mass = max(forced_mass,
+     $                 pmass(i,iconfig)-bwcutoff*pwidth(i,iconfig))
 c RF & TJS, should start from final state particle masses, not only at resonance.
 c Therefore remove the next line.
 c                  xe(i) = max(xe(i),xm(i))
@@ -474,7 +492,8 @@
 
             xo = min(x1,x2)
 
-            xo = xo*xo/stot
+c           Use 1/10000 of sqrt(s) as minimum, to always get integration
+            xo = max(xo*xo/stot,1e-8)
             a=-pmass(i,iconfig)**2/stot
 c            call setgrid(-i,xo,a,pow(i,iconfig))
 
@@ -485,8 +504,8 @@
          endif
       enddo
       if (abs(lpp(1)) .eq. 1 .or. abs(lpp(2)) .eq. 1) then
-         write(*,*) 'etot',etot,nexternal
-         xo = etot**2/stot
+c     Set minimum based on: 1) required energy 2) resonances 3) 1/10000 of sqrt(s)
+         xo = max(max(etot**2, forced_mass**2)/stot,1e-8)
          i = 3*(nexternal-2) - 4 + 1
 c-----------------------
 c     tjs  4/29/2008 use analytic transform for s-hat

=== modified file 'Template/SubProcesses/reweight.f'
--- Template/SubProcesses/reweight.f	2011-04-03 22:06:00 +0000
+++ Template/SubProcesses/reweight.f	2011-04-22 01:24:30 +0000
@@ -373,10 +373,11 @@
       integer idfl, idmap(-nexternal:nexternal)
       integer ipart(2,n_max_cl)
       double precision xnow(2)
-      integer jlast(2),jfirst(2)
+      integer jlast(2),jfirst(2),nwarning
       logical qcdline(2),qcdrad(2)
       logical failed,first
       data first/.true./
+      data nwarning/0/
 
       logical isqcd,isjet,isparton,isjetvx,cluster
       double precision alphas
@@ -585,6 +586,22 @@
      $        q2fact(2)=pt2ijcl(nexternal-2)
       endif
 
+c     Check that factorization scale is >= 2 GeV
+      if(lpp(1).ne.0.and.q2fact(1).lt.4d0.or.
+     $   lpp(2).ne.0.and.q2fact(2).lt.4d0)then
+         if(nwarning.le.10) then
+             nwarning=nwarning+1
+             write(*,*) 'Warning: Too low fact scales: ',
+     $            sqrt(q2fact(1)), sqrt(q2fact(2))
+          endif
+         if(nwarning.eq.11) then
+             nwarning=nwarning+1
+             write(*,*) 'No more warnings written out this run.'
+          endif
+         setclscales=.false.
+         return
+      endif
+
       if (btest(mlevel,3))
      $     write(*,*) 'Set fact scales to ',sqrt(q2fact(1)),sqrt(q2fact(2))
       return

=== modified file 'Template/SubProcesses/setcuts.f'
--- Template/SubProcesses/setcuts.f	2011-04-03 22:06:00 +0000
+++ Template/SubProcesses/setcuts.f	2011-04-22 01:24:30 +0000
@@ -28,6 +28,7 @@
       logical  do_cuts(nexternal)
       integer ncheck
       logical done,fopened
+      logical from_decay(-(nexternal+3):nexternal)
 C     
 C     GLOBAL
 C
@@ -109,7 +110,7 @@
          endif
 c
 c     set ptj and s_min if xqcut and ktscheme = 1, to improve
-c     integration speed, and set drjj to 0.
+c     integration speed, and set drjj and drjl to 0.
 c
         if(xqcut.gt.0.and.ktscheme.eq.1) then
            if(ptj.ge.0d0) ptj=max(ptj, xqcut)
@@ -120,14 +121,23 @@
               write(*,*) 'Warning! drjj > 0 with xqcut > 0, set to 0'
               drjj = 0d0
            endif
+           if(drjl.gt.0d0) then
+              write(*,*) 'Warning! drjl > 0 with xqcut > 0, set to 0'
+              drjl = 0d0
+           endif
         endif
 
+c     Check which particles come from decays
+      if(.not.cut_decays)
+     $       call check_decay(from_decay)
+
 c
 c     check if I have to apply cuts on the particles
 c
       do i=nincoming+1,nexternal
          do_cuts(i)=.true.
          if(nincoming.eq.1) do_cuts(i)=.false.
+         if(.not.cut_decays.and.from_decay(i)) do_cuts(i)=.false.
          is_a_j(i)=.false.
          is_a_l(i)=.false.
          is_a_b(i)=.false.
@@ -518,3 +528,50 @@
 
       end
 
+
+c************************************************************************
+c     Subroutine to check which external particles are from decays
+c************************************************************************
+      subroutine check_decay(from_decay)
+      implicit none
+
+      include 'nexternal.inc'
+      include 'maxconfigs.inc'
+      include 'genps.inc'
+c
+c    Arguments
+c
+      logical from_decay(-(nexternal+3):nexternal)
+c
+c     Global
+c
+      integer iforest(2,-max_branch:-1,lmaxconfigs)
+      common/to_forest/ iforest
+      integer sprop(-max_branch:-1,lmaxconfigs)
+      integer tprid(-max_branch:-1,lmaxconfigs)
+      common/to_sprop/sprop,tprid
+      logical gForceBW(-max_branch:-1,lmaxconfigs)  ! Forced BW
+      include 'decayBW.inc'
+
+c
+c     Local
+c
+      integer i,j
+
+      do i=-(nexternal-3),nexternal
+         from_decay(i)=.false.
+      enddo
+
+c     Set who comes from decay based on forced BW
+      do i=-(nexternal-3),-1
+         if(tprid(i,1).eq.0.and.gForceBW(i,1).or.
+     $        from_decay(i)) then
+            do j=i,-1
+               from_decay(j)=.true.
+               from_decay(iforest(1,j,1))=.true.
+               from_decay(iforest(2,j,1))=.true.
+            enddo
+         endif
+      enddo
+
+      end

=== modified file 'Template/SubProcesses/symmetry.f'
--- Template/SubProcesses/symmetry.f	2011-04-08 00:33:38 +0000
+++ Template/SubProcesses/symmetry.f	2011-04-22 01:24:30 +0000
@@ -164,7 +164,8 @@
 c     Local
 c
       character*30 fname
-      integer ic
+      integer ic, npos
+      character*10 formstr
 
       data ic/0/
 c-----
@@ -172,13 +173,10 @@
 c-----
       ic=ic+1
       fname='ajob'
-      if (ic .lt. 10) then
-         write(fname(5:5),'(i1)') ic
-      elseif (ic .lt. 100) then
-         write(fname(5:6),'(i2)') ic
-      elseif (ic .lt. 1000) then
-         write(fname(5:7),'(i3)') ic
-      endif
+c     Write ic with correct number of digits
+      npos=int(dlog10(dble(ic)))+1
+      write(formstr,'(a,i1,a)') '(I',npos,')'
+      write(fname(5:(5+npos-1)),formstr) ic
       open (unit=lun, file = fname, status='unknown')
       write(lun,15) '#!/bin/bash'
       write(lun,15) '#PBS -q ' // PBS_QUE
@@ -226,6 +224,9 @@
 c     local
 c
       integer i, j, nbw, ic, icode
+      integer ncode, nconf
+      double precision dconfig
+      character*10 formstr
       integer iarray(imax)
       logical lconflict(-max_branch:nexternal)
       logical done
@@ -237,6 +238,8 @@
 c-----
       call open_bash_file(26)
       ic = 0      
+c     ncode is number of digits needed for the code
+      ncode=int(dlog10(3d0)*(max_particles-3))+1
       do i=1,mapconfig(0)
          if (use_config(i) .gt. 0) then
             call bw_conflict(i,iforest(1,-max_branch,i),lconflict,
@@ -273,26 +276,17 @@
                   ic = 1
                endif
 c               write(*,*) 'mapping',ic,mapconfig(i)
-               if (mapconfig(i) .lt. 10) then
-                  write(26,'(i1$)') mapconfig(i)
-               elseif (mapconfig(i) .lt. 100) then
-                  write(26,'(i2$)') mapconfig(i)
-               elseif (mapconfig(i) .lt. 1000) then
-                  write(26,'(i3$)') mapconfig(i)
-               elseif (mapconfig(i) .lt. 10000) then
-                  write(26,'(i4$)') mapconfig(i)
-               endif
+               nconf=int(dlog10(dble(mapconfig(i))))+1
                if (icode .eq. 0) then
-c                 write(26,'($a)') '.000'
-               elseif (icode .lt. 10) then
-                  write(26,'(a,i1$)') '.00', icode
-               elseif (icode .lt. 100) then
-                  write(26,'(a,i2$)') '.0', icode
-               elseif (icode .lt. 1000) then
-                  write(26,'(a,i3$)') '.', icode
+c                 Create format string based on number of digits
+                  write(formstr,'(a,i1,a)') '(I',nconf,'$)'
+                  write(26,formstr) mapconfig(i)
                else
-                  write(*,*) 'Error too many B.W. in symmetry.f',icode
-                  stop
+c                 Create format string based on number of digits
+                  dconfig=mapconfig(i)+icode*1d0/10**ncode
+                  write(formstr,'(a,i1,a,i1,a)') '(F',nconf+ncode+1,
+     $                 '.',ncode,'$)'
+                  write(26,formstr) dconfig
                endif
                write(26,'(a$)') ' '
                call bw_increment_array(iarray,imax,ibase,gForceBW(-imax,i),done)
@@ -300,8 +294,8 @@
          endif
       enddo
       call close_bash_file(26)
-      if (mapconfig(0) .gt. 9999) then
-         write(*,*) 'Only writing first 9999 jobs',mapconfig(0)
+      if (mapconfig(0) .gt. 99999) then
+         write(*,*) 'Only writing first 99999 jobs',mapconfig(0)
       endif
 c
 c     Now write out the symmetry factors for each graph
@@ -337,8 +331,11 @@
             do while (.not. done)
                call enCode(icode,iarray,ibase,imax)
                if (icode .gt. 0) then
-                  write(26,'(f9.3,i6)') mapconfig(i)+real(icode)/1000.,
-     $                 use_config(i)
+                  nconf=int(dlog10(dble(mapconfig(i))))+1
+                  write(formstr,'(a,i1,a,i1,a)') '(F',nconf+ncode+1,
+     $                 '.',ncode,',i6)'
+                  dconfig=mapconfig(i)+icode*1d0/10**ncode
+                  write(26,formstr) dconfig,use_config(i)
                else
                   write(26,'(2i6)') mapconfig(i),use_config(i)
                endif
@@ -401,6 +398,8 @@
 c      
       do i=1,nexternal
          xmass(i) = 0d0
+      enddo
+      do i=1,nexternal-1
          lconflict(-i) = .false.
          iden_part(-i)=0
       enddo
@@ -414,19 +413,27 @@
          if (pwidth(-i,iconfig) .gt. 0d0) then
 c     JA 3/31/11 Keep track of identical particles (i.e., radiation vertices)
 c     by tracing the particle identity from the external particle.
-            if(itree(1,-i).gt.0.and.
-     $           sprop(-i).eq.idup(itree(1,-i),1,1).or.
-     $         itree(2,-i).gt.0.and.
-     $           sprop(-i).eq.idup(itree(2,-i),1,1).or.
-     $         itree(1,-i).lt.0.and.(iden_part(itree(1,-i)).ne.0.and.
+            if(itree(1,-i).gt.0) then
+               if(sprop(-i).eq.idup(itree(1,-i),1,1))
+     $              iden_part(-i) = sprop(-i)
+            endif
+            if(itree(2,-i).gt.0) then
+               if(sprop(-i).eq.idup(itree(2,-i),1,1))
+     $              iden_part(-i) = sprop(-i)
+            endif
+            if(itree(1,-i).lt.0) then
+               if(iden_part(itree(1,-i)).ne.0.and.
      $           sprop(-i).eq.iden_part(itree(1,-i)) .or.
      $         forcebw(itree(1,-i)).and.
-     $           sprop(-i).eq.sprop(itree(1,-i))).or.
-     $         itree(2,-i).lt.0.and.(iden_part(itree(2,-i)).ne.0.and.
+     $           sprop(-i).eq.sprop(itree(1,-i)))
+     $              iden_part(-i) = sprop(-i)
+            endif
+            if(itree(2,-i).lt.0) then
+               if(iden_part(itree(2,-i)).ne.0.and.
      $           sprop(-i).eq.iden_part(itree(2,-i)).or.
      $         forcebw(itree(2,-i)).and.
-     $           sprop(-i).eq.sprop(itree(2,-i))))then
-               iden_part(-i) = sprop(-i)
+     $           sprop(-i).eq.sprop(itree(2,-i)))
+     $              iden_part(-i) = sprop(-i)
             endif
             if (xmass(-i) .gt. pmass(-i,iconfig) .and.
      $           iden_part(-i).eq.0) then !Can't be on shell, and not radiation

=== modified file 'Template/SubProcesses/symmetry_v4.f'
--- Template/SubProcesses/symmetry_v4.f	2011-04-03 22:06:00 +0000
+++ Template/SubProcesses/symmetry_v4.f	2011-04-22 01:24:30 +0000
@@ -15,8 +15,6 @@
       
       double precision ZERO
       parameter       (ZERO = 0d0)
-      integer   maxswitch
-      parameter(maxswitch=99)
       integer lun
       parameter (lun=28)
 c
@@ -26,7 +24,9 @@
       integer mapconfig(0:lmaxconfigs)
       integer sprop(-max_branch:-1,lmaxconfigs)
       integer itree(2,-max_branch:-1)
-      integer imatch
+      integer imatch, nperm
+      integer   maxswitch
+      parameter (maxswitch=10000000)
       integer use_config(0:lmaxconfigs)
       integer i,j, k, n, nsym
       double precision diff
@@ -36,11 +36,11 @@
 c
 c     Local for generating amps
 c
-      double precision p(0:3,99), wgt, x(99), fx
-      double precision p1(0:3,99)
+      double precision p(0:3,nexternal), wgt, x(maxinvar), fx
+      double precision p1(0:3,nexternal)
       integer ninvar, ndim, iconfig, minconfig, maxconfig
       integer ncall,itmax,nconfigs,ntry, ngraphs
-      integer ic(nexternal,maxswitch), jc(12),nswitch
+      integer ic(nexternal), jc(12),nswitch
       double precision saveamp(maxamps)
       integer nmatch, ibase
       logical mtc, even
@@ -99,10 +99,10 @@
       call printout
       call run_printout
 c
-      ndim = 22
+      ndim = maxdim
       ncall = 10000
       itmax = 10
-      ninvar = 35
+      ninvar = maxinvar
       nconfigs = 1
       iconfig=1
       include 'props.inc'
@@ -149,46 +149,52 @@
 c
       nswitch = 1
       do k=1,nexternal
-         ic(k,1)=k
+         ic(k)=k
       enddo
       nmatch = 0
       mtc=.false.
       nsym = 1
       write(lun,'(a,i2,a$)') '       data (isym(i,',nsym,'),i=1,nexternal) /1 '
       do i=2,nexternal
-         write(lun,'(a,i2$)') ",",ic(i,1)
-c         write(*,'(a,i2$)') ",",ic(i,1)
+         write(lun,'(a,i2$)') ",",ic(i)
+c         write(*,'(a,i2$)') ",",ic(i)
       enddo
 c      write(*,*)
       write(lun,'(a)') "/"
-      call nexper(nexternal-2,ic(3,1),mtc,even)
-c      write(*,*) 'mtc',mtc, (ic(i,1),i=1,nexternal)
+      call nexper(nexternal-2,ic(3),mtc,even)
+c      write(*,*) 'mtc',mtc, (ic(i),i=1,nexternal)
+      nperm=0
       do while(mtc)
-         call nexper(nexternal-2,ic(3,1),mtc,even)
-c         write(*,*) 'mtc',mtc, (ic(i,1),i=1,nexternal)
+         nperm=nperm+1
+         if (nperm.gt.maxswitch) goto 10
+         call nexper(nexternal-2,ic(3),mtc,even)
+c         write(*,*) 'mtc',mtc, (ic(i),i=1,nexternal)
          do j=3,nexternal
-            ic(j,1)=ic(j,1)+2
+            ic(j)=ic(j)+2
          enddo
-         if (mtc) then
-             CALL SWITCHMOM(P,P1,IC(1,1),JC,NEXTERNAL)
+         if (.not.check_swap(ic(1))) then
+            do j=3,nexternal
+               ic(j)=ic(j)-2
+            enddo
+            cycle
+         endif
+         CALL SWITCHMOM(P,P1,IC(1),JC,NEXTERNAL)
 c
 c     Now check if it is a valid swap to make
 c
-
-         if (check_swap(ic(1,1))) then
-            write(*,*) 'Good swap', (ic(i,1),i=1,nexternal)
-            nsym=nsym+1
-            write(lun,'(a,i2,a$)') '       data (isym(i,',nsym,'),i=1,nexternal) /1 '
-            do i=2,nexternal
-               write(lun,'(a,i2$)') ",",ic(i,1)
-            enddo
-            write(lun,'(a)') "/"
-
-            call smatrix(p1,fx) 
+         write(*,*) 'Good swap', (ic(i),i=1,nexternal)
+         nsym=nsym+1
+         write(lun,'(a,i2,a$)') '       data (isym(i,',nsym,'),i=1,nexternal) /1 '
+         do i=2,nexternal
+            write(lun,'(a,i2$)') ",",ic(i)
+         enddo
+         write(lun,'(a)') "/"
+
+         call smatrix(p1,fx) 
 c            write(*,*) 'dsig=',fx
-c         write(*,'(7i6)') i,(ic(k,i),k=1,nexternal)
+c         write(*,'(7i6)') i,(ic(k),k=1,nexternal)
 c         do n = i+1,nswitch
-c            CALL SWITCHMOM(P,P1,IC(1,n),JC,NEXTERNAL)
+c            CALL SWITCHMOM(P,P1,IC(1),JC,NEXTERNAL)
 c            call smatrix(p1,fx)
 c
 c        Look for matches, but only for diagrams < current diagram
@@ -213,15 +219,11 @@
                endif
             enddo
          enddo
-         else
-            write(*,*) 'Bad swap', (ic(i,1),i=1,nexternal)
-         endif   !Good Swap
-         endif   !Real Swap
          do j=3,nexternal
-            ic(j,1)=ic(j,1)-2
+            ic(j)=ic(j)-2
          enddo
       enddo
-      write(lun,*) '      data nsym /',nsym,'/'
+ 10   write(lun,*) '      data nsym /',nsym,'/'
       close(lun)
       write(*,*) 'Found ',nmatch, ' matches. ',mapconfig(0)-nmatch,
      $     ' channels remain for integration.'
@@ -396,7 +398,8 @@
 c     Local
 c
       character*30 fname
-      integer ic
+      integer ic, npos
+      character*10 formstr
 
       data ic/0/
 c-----
@@ -404,13 +407,10 @@
 c-----
       ic=ic+1
       fname='ajob'
-      if (ic .lt. 10) then
-         write(fname(5:5),'(i1)') ic
-      elseif (ic .lt. 100) then
-         write(fname(5:6),'(i2)') ic
-      elseif (ic .lt. 1000) then
-         write(fname(5:7),'(i3)') ic
-      endif
+c     Write ic with correct number of digits
+      npos=int(dlog10(dble(ic)))+1
+      write(formstr,'(a,i1,a)') '(I',npos,')'
+      write(fname(5:(5+npos-1)),formstr) ic
       open (unit=lun, file = fname, status='unknown')
       write(lun,15) '#!/bin/bash'
       write(lun,15) '#PBS -q ' // PBS_QUE
@@ -458,6 +458,9 @@
 c     local
 c
       integer i, j, nbw, ic, icode
+      integer ncode, nconf
+      double precision dconfig
+      character*10 formstr
       integer iarray(imax)
       logical lconflict(-max_branch:nexternal)
       logical done
@@ -469,6 +472,8 @@
 c-----
       call open_bash_file(26)
       ic = 0      
+c     ncode is number of digits needed for the code
+      ncode=int(dlog10(3d0)*(max_particles-3))+1
       do i=1,mapconfig(0)
          if (use_config(i) .gt. 0) then
             call bw_conflict(i,iforest(1,-max_branch,i),lconflict,
@@ -505,26 +510,17 @@
                   ic = 1
                endif
 c               write(*,*) 'mapping',ic,mapconfig(i)
-               if (mapconfig(i) .lt. 10) then
-                  write(26,'(i1$)') mapconfig(i)
-               elseif (mapconfig(i) .lt. 100) then
-                  write(26,'(i2$)') mapconfig(i)
-               elseif (mapconfig(i) .lt. 1000) then
-                  write(26,'(i3$)') mapconfig(i)
-               elseif (mapconfig(i) .lt. 10000) then
-                  write(26,'(i4$)') mapconfig(i)
-               endif
-               if (icode .eq. 0) then
-c                 write(26,'($a)') '.000'
-               elseif (icode .lt. 10) then
-                  write(26,'(a,i1$)') '.00', icode
-               elseif (icode .lt. 100) then
-                  write(26,'(a,i2$)') '.0', icode
-               elseif (icode .lt. 1000) then
-                  write(26,'(a,i3$)') '.', icode
+               nconf=int(dlog10(dble(mapconfig(i))))+1
+               if(icode.eq.0)then
+c                 Create format string based on number of digits
+                  write(formstr,'(a,i1,a)') '(I',nconf,'$)'
+                  write(26,formstr) mapconfig(i)
                else
-                  write(*,*) 'Error too many B.W. in symmetry.f',icode
-                  stop
+c                 Create format string based on number of digits
+                  dconfig=mapconfig(i)+icode*1d0/10**ncode
+                  write(formstr,'(a,i1,a,i1,a)') '(F',nconf+ncode+1,
+     $                 '.',ncode,'$)'
+                  write(26,formstr) dconfig
                endif
                write(26,'(a$)') ' '
                call bw_increment_array(iarray,imax,ibase,gForceBW(-imax,i),done)
@@ -532,9 +528,6 @@
          endif
       enddo
       call close_bash_file(26)
-      if (mapconfig(0) .gt. 9999) then
-         write(*,*) 'Only writing first 9999 jobs',mapconfig(0)
-      endif
 c
 c     Now write out the symmetry factors for each graph
 c
@@ -569,8 +562,11 @@
             do while (.not. done)
                call enCode(icode,iarray,ibase,imax)
                if (icode .gt. 0) then
-                  write(26,'(f9.3,i6)') mapconfig(i)+real(icode)/1000.,
-     $                 use_config(i)
+                  nconf=int(dlog10(dble(mapconfig(i))))+1
+                  write(formstr,'(a,i1,a,i1,a)') '(F',nconf+ncode+1,
+     $                 '.',ncode,',i6)'
+                  dconfig=mapconfig(i)+icode*1d0/10**ncode
+                  write(26,formstr) dconfig,use_config(i)
                else
                   write(26,'(2i6)') mapconfig(i),use_config(i)
                endif
@@ -632,6 +628,8 @@
 c      
       do i=1,nexternal
          xmass(i) = 0d0
+      enddo
+      do i=1,nexternal-1
          lconflict(-i) = .false.
       enddo
 c
@@ -644,19 +642,27 @@
          if (pwidth(-i,iconfig) .gt. 0d0) then
 c     JA 3/31/11 Keep track of identical particles (i.e., radiation vertices)
 c     by tracing the particle identity from the external particle.
-            if(itree(1,-i).gt.0.and.
-     $           sprop(-i).eq.idup(itree(1,-i),1,1).or.
-     $         itree(2,-i).gt.0.and.
-     $           sprop(-i).eq.idup(itree(2,-i),1,1).or.
-     $         itree(1,-i).lt.0.and.(iden_part(itree(1,-i)).ne.0.and.
+            if(itree(1,-i).gt.0) then
+               if(sprop(-i).eq.idup(itree(1,-i),1,1))
+     $              iden_part(-i) = sprop(-i)
+            endif
+            if(itree(2,-i).gt.0) then
+               if(sprop(-i).eq.idup(itree(2,-i),1,1))
+     $              iden_part(-i) = sprop(-i)
+            endif
+            if(itree(1,-i).lt.0) then
+               if(iden_part(itree(1,-i)).ne.0.and.
      $           sprop(-i).eq.iden_part(itree(1,-i)) .or.
      $         forcebw(itree(1,-i)).and.
-     $           sprop(-i).eq.sprop(itree(1,-i))).or.
-     $         itree(2,-i).lt.0.and.(iden_part(itree(2,-i)).ne.0.and.
+     $           sprop(-i).eq.sprop(itree(1,-i)))
+     $              iden_part(-i) = sprop(-i)
+            endif
+            if(itree(2,-i).lt.0) then
+               if(iden_part(itree(2,-i)).ne.0.and.
      $           sprop(-i).eq.iden_part(itree(2,-i)).or.
      $         forcebw(itree(2,-i)).and.
-     $           sprop(-i).eq.sprop(itree(2,-i))))then
-               iden_part(-i) = sprop(-i)
+     $           sprop(-i).eq.sprop(itree(2,-i)))
+     $              iden_part(-i) = sprop(-i)
             endif
             if (xmass(-i) .gt. pmass(-i,iconfig) .and.
      $           iden_part(-i).eq.0) then !Can't be on shell, and not radiation

=== modified file 'Template/SubProcesses/unwgt.f'
--- Template/SubProcesses/unwgt.f	2010-12-03 19:28:13 +0000
+++ Template/SubProcesses/unwgt.f	2011-04-22 01:24:30 +0000
@@ -391,6 +391,7 @@
 
       double precision sscale,aaqcd,aaqed
       integer ievent,npart
+      logical flip
 
       real ran1
       external ran1
@@ -411,7 +412,7 @@
       Double Precision amp2(maxamps), jamp2(0:maxflow)
       common/to_amps/  amp2,       jamp2
 
-      character*79       hel_buf
+      character*101       hel_buf
       common/to_helicity/hel_buf
 
       integer           mincfig, maxcfig
@@ -501,7 +502,7 @@
 
 c   Set helicities
 c      write(*,*) 'Getting helicity',hel_buf(1:50)
-      read(hel_buf,'(15i5)') (jpart(7,isym(i, jsym)),i=1,nexternal)
+      read(hel_buf,'(20i5)') (jpart(7,isym(i, jsym)),i=1,nexternal)
 c      write(*,*) 'ihel',jpart(7,1),jpart(7,2)
 
 c   Fix ordering of ptclus
@@ -513,20 +514,14 @@
       enddo
 
 c     Check if we have flipped particle 1 and 2, and flip back
+      flip = .false.
       if (p(3,1).lt.0) then
          do j=0,3
             pdum(j)=p(j,1)
             p(j,1)=p(j,2)
             p(j,2)=pdum(j)
          enddo
-         do i=1,7
-            j=jpart(i,1)
-            jpart(i,1)=jpart(i,2)
-            jpart(i,2)=j
-         enddo
-         ptcltmp(1)=ptclus(1)
-         ptclus(1)=ptclus(2)
-         ptclus(2)=ptcltmp(1)
+         flip = .true.
       endif
 
 c
@@ -556,6 +551,18 @@
       call addmothers(ip,jpart,pb,isym,jsym,sscale,aaqcd,aaqed,buff,
      $                npart,numproc)
 
+c     Need to flip after addmothers, since color might get overwritten
+      if (flip) then
+         do i=1,7
+            j=jpart(i,1)
+            jpart(i,1)=jpart(i,2)
+            jpart(i,2)=j
+         enddo
+         ptcltmp(1)=ptclus(1)
+         ptclus(1)=ptclus(2)
+         ptclus(2)=ptcltmp(1)
+      endif
+
 c
 c     Write events to lun
 c

=== modified file 'Template/bin/gen_crossxhtml-pl'
--- Template/bin/gen_crossxhtml-pl	2011-04-04 01:05:14 +0000
+++ Template/bin/gen_crossxhtml-pl	2011-04-22 01:24:30 +0000
@@ -98,7 +98,7 @@
   print PAGE  "<HTML> <HEAD> \n";
 
   if (-e "RunWeb") {
-      print PAGE "<META HTTP-EQUIV=\"Refresh\" CONTENT=\"30\" > ";
+      print PAGE "<META HTTP-EQUIV=\"Refresh\" CONTENT=\"10\" > ";
   }
   print PAGE "<META HTTP-EQUIV=\"EXPIRES\" CONTENT=\"20\" > ";
   print PAGE "\n <TITLE>Online Event Generation</TITLE></HEAD> \n <BODY>";

=== modified file 'Template/bin/generate_events'
--- Template/bin/generate_events	2011-04-04 01:05:14 +0000
+++ Template/bin/generate_events	2011-04-22 01:24:30 +0000
@@ -201,7 +201,6 @@
    mv survey refine
    $dirbin/gen_crossxhtml-pl $t
    $dirbin/refine $a $mode $n 5 $t
-   #$dirbin/refine 0.01 $mode $n
    mv refine refine2
    $dirbin/gen_crossxhtml-pl $t
    $dirbin/refine $a $mode $n 1 $t

=== added file 'UpdateNotes.txt'
--- UpdateNotes.txt	1970-01-01 00:00:00 +0000
+++ UpdateNotes.txt	2011-04-22 01:24:30 +0000
@@ -0,0 +1,36 @@
+Update notes for MadGraph 5 (in reverse time order)
+
+1.1.0 (21/04/11): JA: Removed hard limit on number of external particles in 
+                      MadEvent, allowing for unlimited length decay chains there
+                      (up to 14 final state particles successfully integrated).
+                  JA: Improved helicity selection and automatic full helicity 
+                      sum if needed. Optimization of run parameters.
+                  JA: New flag in run_card.dat to decide whether basic cuts
+                      are applied to decay products or not.
+                  OM: Merged ALOHA calls for different lorentz structures 
+                      with the same color structures, increasing the speed and 
+                      efficiency of matrix element evaluations.
+                  OM: Added new "open" command in command line interface,
+                      allowing to open standard file types directly.
+                      Automatically open crossx.html at launch.
+                  JA: Fixed MadEvent bugs for multiparton processes with 
+                      conflicting decays and some faulty array limits.
+                  JA: Suppressed scary but irrelevant warnings for compiling 
+                      2->1 and 1->2 processes in MadEvent.
+                  JA: Pythia 8 output further optimized.
+                  JA: Several minor fixes relating to user interface etc.
+
+1.0.0 (12/04/11): Official release of MadGraph 5. Some of the features:
+                  - Complete FeynRules compatibility through the UFO interface
+                  - Automatic writing of HELAS routines for any model in
+                    Fortran, C++ or Python through ALOHA
+                  - Matrix element output in Fortran, C++ or Python
+                  - Output formats: MadEvent, Pythia 8, Standalone (Fortran/C++)
+                  - Support for spin 0, 1/2, 1, 2 particles
+                  - Support for color 1, 3, 6, 8
+                  - Revamped MadEvent with improved subprocess directory 
+                    organization and vastly increased speed and stability
+                  - Unlimited length decay chains (up to 12 final state
+                    particles tested with MadEvent, see v. 1.0.1)
+                  - Process checks for new model implementations
+                  - ...and much more (see paper "MadGraph 5: Going Beyond")

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2011-04-01 22:38:18 +0000
+++ aloha/aloha_writers.py	2011-04-22 01:24:30 +0000
@@ -19,10 +19,16 @@
     
     def __init__(self, abstract_routine, dirpath):
 
+
+        name = get_routine_name(abstract_routine.name, abstract_routine.outgoing)
+        if dirpath:
+            self.dir_out = dirpath
+            self.out_path = os.path.join(dirpath, name + self.extension)
+        else:
+            self.out_path = None
+            self.dir_out = None
+
         self.obj = abstract_routine.expr
-        name = get_routine_name(abstract_routine.name, abstract_routine.outgoing)
-        self.out_path = os.path.join(dirpath, name + self.extension)
-        self.dir_out = dirpath
         self.particles =  [self.type_to_variable[spin] for spin in \
                           abstract_routine.spins]
         self.namestring = name
@@ -36,6 +42,9 @@
         self.make_all_lists()    # Compute the expression for the call ordering
                                  #the definition of objects,..
                                  
+                                 
+                                 
+                                 
     def pass_to_HELAS(self, indices, start=0):
         """find the Fortran HELAS position for the list of index""" 
         
@@ -170,18 +179,19 @@
                 
         return call_arg
 
-    def reorder_call_list(self, call_list, old, new):
-        """ restore the correct order for symmetries """
-        spins = self.particles
-        assert(0 < old < new)
-        old, new = old -1, new -1 # pass in real position in particles list
-        assert(spins[old] == spins[new])
-        spin =spins[old]
-        
-        new_call = call_list[:]
-        val = new_call.pop(old)
-        new_call.insert(new - 1, val)
-        return new_call
+#    def reorder_call_list(self, call_list, old, new):
+#        """ restore the correct order for symmetries """
+#        raise
+#        #spins = self.particles
+#        #assert(0 < old < new)
+#        #old, new = old -1, new -1 # pass in real position in particles list
+#        #assert(spins[old] == spins[new])
+#        #spin =spins[old]
+#        
+#        new_call = call_list[:]
+#        #val = new_call.pop(old)
+#        #new_call.insert(new - 1, val)
+#        return new_call
     
         
     def make_momentum_conservation(self):
@@ -249,42 +259,45 @@
                     'V':'double complex V%d(6)',
                     'T':'double complex T%s(18)'}
     
-    def define_header(self):
+    def define_header(self, name=None):
         """Define the Header of the fortran file. This include
             - function tag
             - definition of variable
         """
-            
+        
+        if not name:
+            name = self.namestring
+             
         Momenta = self.collected['momenta']
         OverM = self.collected['om']
         
         CallList = self.calllist['CallList']
         declare_list = self.calllist['DeclareList']
-        if 'double complex C' in declare_list:
+        if 'double complex COUP' in declare_list:
             alredy_update = True
         else:
             alredy_update = False
         
         if not alredy_update:    
-            declare_list.append('double complex C')
+            declare_list.append('double complex COUP')
 
         # define the type of function and argument        
         if not self.offshell:
             str_out = 'subroutine %(name)s(%(args)s,vertex)\n' % \
-               {'name': self.namestring,
-                'args': ','.join(CallList+ ['C']) }
+               {'name': name,
+                'args': ','.join(CallList+ ['COUP']) }
             if not alredy_update: 
-                declare_list.append('double complex vertex\n') 
+                declare_list.append('double complex vertex') 
         else:
             if not alredy_update:
                 declare_list.append('double complex denom')
                 declare_list.append('double precision M%(id)d, W%(id)d' % 
                                                           {'id': self.offshell})
-            call_arg = '%(args)s, C, M%(id)d, W%(id)d, %(spin)s%(id)d' % \
+            call_arg = '%(args)s, COUP, M%(id)d, W%(id)d, %(spin)s%(id)d' % \
                     {'args': ', '.join(CallList), 
                      'spin': self.particles[self.offshell -1],
                      'id': self.offshell}
-            str_out = ' subroutine %s(%s)\n' % (self.namestring, call_arg) 
+            str_out = ' subroutine %s(%s)\n' % (name, call_arg) 
 
         # Forcing implicit None
         str_out += 'implicit none \n'
@@ -402,7 +415,7 @@
         OutString = ''
         if not self.offshell:
             for ind in self.obj.listindices():
-                string = 'Vertex = C*' + self.write_obj(self.obj.get_rep(ind))
+                string = 'Vertex = COUP*' + self.write_obj(self.obj.get_rep(ind))
                 string = string.replace('+-', '-')
                 string = re.sub('\((?P<num>[+-]*[0-9])(?P<num2>[+-][0-9])[Jj]\)\.', '(\g<num>d0,\g<num2>d0)', string)
                 string = re.sub('(?P<num>[0-9])[Jj]\.', '\g<num>.*(0d0,1d0)', string)
@@ -420,7 +433,7 @@
             string = re.sub('(?P<num>[0-9])[Jj]\.', '\g<num>*(0d0,1d0)', string)
             OutString = OutString + string + '\n'
             for ind in numerator.listindices():
-                string = '%s(%d)= C*denom*' % (OffShellParticle, self.pass_to_HELAS(ind, start=1))
+                string = '%s(%d)= COUP*denom*' % (OffShellParticle, self.pass_to_HELAS(ind, start=1))
                 string += self.write_obj(numerator.get_rep(ind))
                 string = string.replace('+-', '-')
                 string = re.sub('\((?P<num>[+-][0-9])\+(?P<num2>[+-][0-9])[Jj]\)\.', '(\g<num>d0,\g<num2>d0)', string)
@@ -430,45 +443,178 @@
     
     def define_symmetry(self, new_nb):
         number = self.offshell 
-        calls = self.reorder_call_list(self.calllist['CallList'], self.offshell,
-                                                                        new_nb)
-        Outstring = 'call '+self.namestring+'('+','.join(calls)+',C,M%s,W%s,%s%s)' \
+        calls = self.calllist['CallList']
+                                                                
+        Outstring = 'call '+self.namestring+'('+','.join(calls)+',COUP,M%s,W%s,%s%s)\n' \
                          %(number,number,self.particles[number-1],number)
         return Outstring
     
     def define_foot(self):
-        return 'end' 
+        return 'end\n\n' 
 
     def write(self, mode='self'):
-                
-        writer = writers.FortranWriter(self.out_path)
-        writer.downcase = False 
-        commentstring = 'This File is Automatically generated by ALOHA \n'
-        commentstring += 'The process calculated in this file is: \n'
-        commentstring += self.comment + '\n'
-        writer.write_comments(commentstring)
-         
+                         
         # write head - momenta - body - foot
-        writer.writelines(self.define_header())
-        writer.writelines(self.define_momenta())
-        writer.writelines(self.define_expression())
-        writer.writelines(self.define_foot())
+        text = self.define_header()+'\n'
+        text += self.define_momenta()+'\n'
+        text += self.define_expression()
+        text += self.define_foot()
         
+        sym_text = []
         for elem in self.symmetries: 
             symmetryhead = self.define_header().replace( \
                              self.namestring,self.namestring[0:-1]+'%s' %(elem))
             symmetrybody = self.define_symmetry(elem)
-            writer.write_comments('\n%s\n' % ('#'*65))
-            writer.writelines(symmetryhead)
-            writer.writelines(symmetrybody)
-            writer.writelines(self.define_foot())
-        
-        
+            
+            sym_text.append(symmetryhead + symmetrybody + self.define_foot())
+            
+            
+        if self.out_path:
+            writer = writers.FortranWriter(self.out_path)
+            writer.downcase = False 
+            commentstring = 'This File is Automatically generated by ALOHA \n'
+            commentstring += 'The process calculated in this file is: \n'
+            commentstring += self.comment + '\n'
+            writer.write_comments(commentstring)
+            writer.writelines(text)
+            for text in sym_text:
+                writer.write_comments('\n%s\n' % ('#'*65))
+                writer.writelines(text)
+        else:
+            for stext in sym_text:
+                text += '\n\n' + stext
+
+        return text + '\n'
+            
+    def write_combined(self, lor_names, mode='self', offshell=None):
+        """Write routine for combine ALOHA call (more than one coupling)"""
+        
+        # Set some usefull command
+        if offshell is None:
+            sym = 1 
+            offshell = self.offshell
+        else:
+            sym = None  # deactivate symetry
+            
+        name = combine_name(self.abstractname, lor_names, offshell)
+
+                 
+        # write header 
+        header = self.define_header(name=name)
+        new_couplings = ['COUP%s' % (i+1) for i in range(len(lor_names)+1)]
+        text = header.replace('COUP', ','.join(new_couplings))
+        # define the TMP for storing output        
+        if not offshell:
+            text += ' double complex TMP\n'
+        else:
+            spin = self.particles[offshell -1] 
+            text += ' double complex TMP(%s)\n integer i' % self.type_to_size[spin]         
+        
+        # Define which part of the routine should be called
+        addon = ''
+        if 'C' in self.namestring:
+            short_name, addon = name.split('C',1)
+            if addon.split('_')[0].isdigit():
+                addon = 'C' +self.namestring.split('C',1)[1]
+            else:
+                addon = '_%s' % self.offshell
+        else:
+            addon = '_%s' % self.offshell
+
+        # how to call the routine
+        if not offshell:
+            main = 'vertex'
+            call_arg = '%(args)s, %(COUP)s, %(LAST)s' % \
+                    {'args': ', '.join(self.calllist['CallList']), 
+                     'COUP':'COUP%d',
+                     'spin': self.particles[self.offshell -1],
+                     'LAST': '%s'}
+        else:
+            main = '%(spin)s%(id)d' % \
+                          {'spin': self.particles[offshell -1],
+                           'id': self.offshell}
+            call_arg = '%(args)s, %(COUP)s, M%(id)d, W%(id)d, %(LAST)s' % \
+                    {'args': ', '.join(self.calllist['CallList']), 
+                     'COUP':'COUP%d',
+                     'id': self.offshell,
+                     'LAST': '%s'}
+
+        # make the first call
+        line = " CALL %s%s("+call_arg+")\n"
+        text += '\n\n' + line % (self.namestring, '', 1, main)
+        
+        # make the other call
+        for i,lor in enumerate(lor_names):
+            text += line % (lor, addon, i+2, 'TMP')
+            if not offshell:
+                text += ' vertex = vertex + tmp\n'
+            else:
+                size = self.type_to_size[spin] -2
+                text += """ do i=1,%(id)d
+                %(main)s(i) = %(main)s(i) + tmp(i)
+                enddo\n""" %  {'id': size, 'main':main}
+                
+
+        text += self.define_foot()
+        
+        #ADD SYMETRY
+        if sym:
+            for elem in self.symmetries:
+                text += self.write_combined(lor_names, mode, elem)
+            
+            
+        if self.out_path:
+            # Prepare a specific file
+            path = os.path.join(os.path.dirname(self.out_path), name+'.f')
+            writer = writers.FortranWriter(path)
+            writer.downcase = False 
+            commentstring = 'This File is Automatically generated by ALOHA \n'
+            writer.write_comments(commentstring)
+            writer.writelines(text)
+        
+        return text
+    
 def get_routine_name(name,outgoing):
     """ build the name of the aloha function """
     
     return '%s_%s' % (name, outgoing) 
 
+def combine_name(name, other_names, outgoing):
+    """ build the name for combined aloha function """
+
+    # Two possible scheme FFV1C1_2_X or FFV1__FFV2C1_X
+    # If they are all in FFVX scheme then use the first
+    p=re.compile('^(?P<type>[FSVT]+)(?P<id>\d+)')
+    routine = ''
+    if p.search(name):
+        base, id = p.search(name).groups()
+        routine = name
+        for s in other_names:
+            try:
+                base2,id2 = p.search(s).groups()
+            except:
+                routine = ''
+                break # one matching not good -> other scheme
+            if base != base2:
+                routine = ''
+                break  # one matching not good -> other scheme
+            else:
+                routine += '_%s' % id2
+    if routine:
+        return routine +'_%s' % outgoing
+
+    addon = ''
+    if 'C' in name:
+        short_name, addon = name.split('C',1)
+        try:
+            addon = 'C' + str(int(addon))
+        except:
+            addon = ''
+        else:
+            name = short_name
+
+    return '_'.join((name,) + tuple(other_names)) + addon + '_%s' % outgoing
+ 
 class ALOHAWriterForCPP(WriteALOHA): 
     """Routines for writing out helicity amplitudes as C++ .h and .cc files."""
     
@@ -477,17 +623,20 @@
                     'V':'double complex V%d[6]',
                     'T':'double complex T%s[18]'}
     
-    def define_header(self):
+    def define_header(self, name=None):
         """Define the headers for the C++ .h and .cc files. This include
             - function tag
             - definition of variable
             - momentum conservation
         """
+        
+        if name is None:
+            name = self.namestring
             
         #Width = self.collected['width']
         #Mass = self.collected['mass']
         
-        CallList = self.calllist['CallList']
+        CallList = self.calllist['CallList'][:]
         
         local_declare = []
         OffShell = self.offshell
@@ -503,12 +652,12 @@
         # define the type of function and argument
         if not OffShell:
             str_out = 'void %(name)s(%(args)s, complex<double>& vertex)' % \
-               {'name': self.namestring,
-                'args': ','.join(CallList + ['complex<double> C'])}
+               {'name': name,
+                'args': ','.join(CallList + ['complex<double> COUP'])}
         else: 
             str_out = 'void %(name)s(%(args)s, double M%(number)d, double W%(number)d, complex<double>%(out)s%(number)d[])' % \
-              {'name': self.namestring,
-               'args': ','.join(CallList+ ['complex<double> C']),
+              {'name': name,
+               'args': ','.join(CallList+ ['complex<double> COUP']),
                'out': self.particles[OffShellParticle],
                'number': OffShellParticle + 1 
                }
@@ -613,7 +762,7 @@
         OutString = '' 
         if not self.offshell:
             for ind in self.obj.listindices():
-                string = 'vertex = C*' + self.write_obj(self.obj.get_rep(ind))
+                string = 'vertex = COUP*' + self.write_obj(self.obj.get_rep(ind))
                 string = string.replace('+-', '-')
                 OutString = OutString + string + ';\n'
         else:
@@ -626,7 +775,7 @@
             string = string.replace('+-', '-')
             OutString = OutString + string + ';\n'
             for ind in numerator.listindices():
-                string = '%s[%d]= C*denom*' % (OffShellParticle, self.pass_to_HELAS(ind))
+                string = '%s[%d]= COUP*denom*' % (OffShellParticle, self.pass_to_HELAS(ind))
                 string += self.write_obj(numerator.get_rep(ind))
                 string = string.replace('+-', '-')
                 OutString = OutString + string + ';\n' 
@@ -636,8 +785,7 @@
     remove_double = re.compile('complex<double> (?P<name>[\w]+)\[\]')
     def define_symmetry(self, new_nb):
         """Write the call for symmetric routines"""
-        calls = self.reorder_call_list(self.calllist['CallList'],
-                                       self.offshell, new_nb)
+        calls = self.calllist['CallList']
         
         for i, call in enumerate(calls):
             if self.remove_double.match(call):
@@ -647,22 +795,24 @@
         #calls = [self.remove_double.match(call).group('name') for call in \
         #         calls]
         number = self.offshell 
-        Outstring = self.namestring+'('+','.join(calls)+',C,M%s,W%s,%s%s);' \
+        Outstring = self.namestring+'('+','.join(calls)+',COUP,M%s,W%s,%s%s);\n' \
                          %(number,number,self.particles[self.offshell-1],number)
         return Outstring
     
     def define_foot(self):
         """Return the end of the function definition"""
 
-        return '}\n' 
+        return '}\n\n' 
 
-    def write_h(self, header):
+    def write_h(self, header, compiler_cmd=True):
         """Return the full contents of the .h file"""
 
-        h_string = '#ifndef '+ self.namestring + '_guard\n'
-        h_string += '#define ' + self.namestring + '_guard\n'
-        h_string += '#include <complex>\n'
-        h_string += 'using namespace std;\n\n'
+        h_string = ''
+        if compiler_cmd:
+            h_string = '#ifndef '+ self.namestring + '_guard\n'
+            h_string += '#define ' + self.namestring + '_guard\n'
+            h_string += '#include <complex>\n'
+            h_string += 'using namespace std;\n\n'
 
         h_header = header['h_header']
 
@@ -672,15 +822,44 @@
             symmetryhead = h_header.replace( \
                              self.namestring,self.namestring[0:-1]+'%s' %(elem))
             h_string += symmetryhead
-            
-        h_string += '#endif'
+
+        if compiler_cmd:
+            h_string += '#endif\n\n'
 
         return h_string
 
-    def write_cc(self, header):
+    def write_combined_h(self, lor_names, offshell=None, compiler_cmd=True):
+        """Return the content of the .h file linked to multiple lorentz call."""
+        
+        name = combine_name(self.abstractname, lor_names, offshell)
+        text= ''
+        if compiler_cmd:
+            text = '#ifndef '+ name + '_guard\n'
+            text += '#define ' + name + '_guard\n'
+            text += '#include <complex>\n'
+            text += 'using namespace std;\n\n'
+        
+        # write header 
+        header = self.define_header(name=name)
+        h_header = header['h_header']
+        new_couplings = ['COUP%s' % (i+1) for i in range(len(lor_names)+1)]
+        h_string = h_header.replace('COUP', ', complex <double>'.join(new_couplings))
+        text += h_string 
+        
+        for elem in self.symmetries: 
+            text += h_string.replace(name, name[0:-1]+str(elem))
+        
+        if compiler_cmd:
+            text += '#endif'
+        
+        return text
+
+    def write_cc(self, header, compiler_cmd=True):
         """Return the full contents of the .cc file"""
 
-        cc_string = '#include \"%s.h\"\n\n' % self.namestring
+        cc_string = ''
+        if compiler_cmd:
+            cc_string = '#include \"%s.h\"\n\n' % self.namestring
         cc_header = header['cc_header']
         cc_string += cc_header
         cc_string += self.define_momenta()
@@ -696,24 +875,154 @@
             cc_string += self.define_foot()
 
         return cc_string
+
+    def write_combined_cc(self, lor_names, offshell=None, compiler_cmd=True, sym=True):
+        "Return the content of the .cc file linked to multiple lorentz call."
+        
+        # Set some usefull command
+        if offshell is None:
+            offshell = self.offshell
+            
+        name = combine_name(self.abstractname, lor_names, offshell)
+
+        text = ''
+        if compiler_cmd:
+            text += '#include "%s.h"\n\n' % name
+           
+        # write header 
+        header = self.define_header(name=name)['cc_header']
+        new_couplings = ['COUP%s' % (i+1) for i in range(len(lor_names)+1)]
+        text += header.replace('COUP', ', complex<double>'.join(new_couplings))
+        # define the TMP for storing output        
+        if not offshell:
+            text += 'complex<double> tmp;\n'
+        else:
+            spin = self.particles[offshell -1] 
+            text += 'complex<double> tmp[%s];\n int i = 0;' % self.type_to_size[spin]         
+
+        # Define which part of the routine should be called
+        addon = ''
+        if 'C' in self.namestring:
+            short_name, addon = name.split('C',1)
+            if addon.split('_')[0].isdigit():
+                addon = 'C' +self.namestring.split('C',1)[1]
+            else:
+                addon = '_%s' % self.offshell
+        else:
+            addon = '_%s' % self.offshell
+
+        # how to call the routine
+        if not offshell:
+            main = 'vertex'
+            call_arg = '%(args)s, %(COUP)s, %(LAST)s' % \
+                    {'args': ', '.join(self.calllist['CallList']), 
+                     'COUP':'COUP%d',
+                     'spin': self.particles[self.offshell -1],
+                     'LAST': '%s'}
+        else:
+            main = '%(spin)s%(id)d' % \
+                          {'spin': self.particles[self.offshell -1],
+                           'id': self.offshell}
+            call_arg = '%(args)s, %(COUP)s, M%(id)d, W%(id)d, %(LAST)s' % \
+                    {'args': ', '.join(self.calllist['CallList']), 
+                     'COUP':'COUP%d',
+                     'id': self.offshell,
+                     'LAST': '%s'}
+
+        # make the first call
+        line = "%s%s("+call_arg+");\n"
+        text += '\n\n' + line % (self.namestring, '', 1, main)
+        
+        # make the other call
+        for i,lor in enumerate(lor_names):
+            text += line % (lor, addon, i+2, 'tmp')
+            if not offshell:
+                text += ' vertex = vertex + tmp;\n'
+            else:
+                size = self.type_to_size[spin] -2
+                text += """ while (i < %(id)d)
+                {
+                %(main)s[i] = %(main)s[i] + tmp[i];
+                i++;
+                }\n""" %  {'id': size, 'main':main}
+                
+
+        text += self.define_foot()
+        
+        if sym:
+            for elem in self.symmetries:
+                text += self.write_combined_cc(lor_names, elem, 
+                                                      compiler_cmd=compiler_cmd,
+                                                      sym=False)
+            
+            
+        if self.out_path:
+            # Prepare a specific file
+            path = os.path.join(os.path.dirname(self.out_path), name+'.f')
+            writer = writers.FortranWriter(path)
+            writer.downcase = False 
+            commentstring = 'This File is Automatically generated by ALOHA \n'
+            writer.write_comments(commentstring)
+            writer.writelines(text)
+        
+        return text
+
+
     
-    def write(self, mode='self'):
+    def write(self, mode='self', **opt):
         """Write the .h and .cc files"""
 
-        writer_h = writers.CPPWriter(self.out_path + ".h")
-        writer_cc = writers.CPPWriter(self.out_path + ".cc")
-        commentstring = 'This File is Automatically generated by ALOHA \n'
-        commentstring += 'The process calculated in this file is: \n'
-        commentstring += self.comment + '\n'
-        writer_h.write_comments(commentstring)
-        writer_cc.write_comments(commentstring)
-         
+
         # write head - momenta - body - foot
-        
         header = self.define_header()
-        writer_h.writelines(self.write_h(header))
-        writer_cc.writelines(self.write_cc(header))
- 
+        h_text = self.write_h(header, **opt)
+        cc_text = self.write_cc(header, **opt)
+        
+        # write in two file
+        if self.out_path:
+            writer_h = writers.CPPWriter(self.out_path + ".h")
+            writer_cc = writers.CPPWriter(self.out_path + ".cc")
+            commentstring = 'This File is Automatically generated by ALOHA \n'
+            commentstring += 'The process calculated in this file is: \n'
+            commentstring += self.comment + '\n'
+            writer_h.write_comments(commentstring)
+            writer_cc.write_comments(commentstring)
+            writer_h.writelines(h_text)
+            writer_cc.writelines(cc_text)
+            
+        return h_text, cc_text
+ 
+ 
+ 
+    def write_combined(self, lor_names, mode='self', offshell=None, **opt):
+        """Write the .h and .cc files associated to the combined file"""
+
+        # Set some usefull command
+        if offshell is None:
+            offshell = self.offshell
+        
+        name = combine_name(self.abstractname, lor_names, offshell)
+        
+        h_text = self.write_combined_h(lor_names, offshell, **opt)
+        cc_text = self.write_combined_cc(lor_names, offshell, sym=True, **opt)
+        
+        if self.out_path:
+            # Prepare a specific file
+            path = os.path.join(os.path.dirname(self.out_path), name)
+            commentstring = 'This File is Automatically generated by ALOHA \n'
+            
+            writer_h = writers.CPPWriter(path + ".h")
+            writer_h.write_comments(commentstring)
+            writer_h.writelines(h_text)
+            
+            writer_cc = writers.CPPWriter(path + ".cc")
+            writer_cc.write_comments(commentstring)
+            writer_cc.writelines(cc_text)
+        else:
+            return h_text, cc_text
+        
+        return h_text, cc_text
+        
 class ALOHAWriterForPython(WriteALOHA):
     """ A class for returning a file/a string for python evaluation """
     
@@ -723,18 +1032,11 @@
         """ standard init but if dirpath is None the write routine will
         return a string (which can be evaluated by an exec"""
         
-        if dirpath:
-            WriteALOHA.__init__(self, abstract_routine, dirpath)
-        else:
-            WriteALOHA.__init__(self, abstract_routine, '')
-            self.out_path = None
-            self.dir_out = None
+        WriteALOHA.__init__(self, abstract_routine, dirpath)
         self.outname = '%s%s' % (self.particles[self.offshell -1], \
                                                                self.offshell)
     
-    
-    
-    
+        
     def change_var_format(self, name): 
         """Formatting the variable name to Python format
         start to count at zero"""
@@ -759,7 +1061,7 @@
         OutString = ''
         if not self.offshell:
             for ind in self.obj.listindices():
-                string = 'vertex = C*' + self.write_obj(self.obj.get_rep(ind))
+                string = 'vertex = COUP*' + self.write_obj(self.obj.get_rep(ind))
                 string = string.replace('+-', '-')
                 OutString = OutString + string + '\n'
         else:
@@ -773,7 +1075,7 @@
             string = string.replace('+-', '-')
             OutString += string + '\n'
             for ind in numerator.listindices():
-                string = '%s[%d]= C*denom*' % (self.outname, self.pass_to_HELAS(ind))
+                string = '%s[%d]= COUP*denom*' % (self.outname, self.pass_to_HELAS(ind))
                 string += self.write_obj(numerator.get_rep(ind))
                 string = string.replace('+-', '-')
                 OutString += string + '\n' 
@@ -781,17 +1083,19 @@
     
     def define_foot(self):
         if not self.offshell:
-            return 'return vertex'
+            return 'return vertex\n\n'
         else:
-            return 'return %s' % (self.outname)
+            return 'return %s\n\n' % (self.outname)
             
     
-    def define_header(self):
+    def define_header(self, name=None):
         """Define the Header of the fortran file. This include
             - function tag
             - definition of variable
         """
-            
+        if name is None:
+            name = self.namestring
+           
         Momenta = self.collected['momenta']
         OverM = self.collected['om']
         
@@ -801,11 +1105,11 @@
         # define the type of function and argument        
         if not self.offshell:
             str_out += 'def %(name)s(%(args)s):\n' % \
-                {'name': self.namestring,
-                 'args': ','.join(CallList+ ['C']) }
+                {'name': name,
+                 'args': ','.join(CallList+ ['COUP']) }
         else:
-            str_out += 'def %(name)s(%(args)s, C, M%(id)d, W%(id)d):\n' % \
-                {'name': self.namestring,
+            str_out += 'def %(name)s(%(args)s, COUP, M%(id)d, W%(id)d):\n' % \
+                {'name': name,
                  'args': ', '.join(CallList), 
                      'id': self.offshell}            
         return str_out     
@@ -871,9 +1175,8 @@
 
     def define_symmetry(self, new_nb):
         number = self.offshell 
-        calls = self.reorder_call_list(self.calllist['CallList'], self.offshell,
-                                                                        new_nb)
-        Outstring = 'return '+self.namestring+'('+','.join(calls)+',C,M%s,W%s)' \
+        calls = self.calllist['CallList']
+        Outstring = 'return '+self.namestring+'('+','.join(calls)+',COUP,M%s,W%s)\n' \
                          %(number,number)
         return Outstring        
 
@@ -900,5 +1203,82 @@
             
         if self.out_path:
             ff = open(self.out_path,'w').write(text)
-        else:
-            return text
+
+        return text + '\n'
+
+    def write_combined(self, lor_names, mode='self', offshell=None):
+        """Write routine for combine ALOHA call (more than one coupling)"""
+        
+        # Set some usefull command
+        if offshell is None:
+            sym = 1
+            offshell = self.offshell  
+        else:
+            sym = None
+        name = combine_name(self.abstractname, lor_names, offshell)
+
+        # write head - momenta - body - foot
+        text = ''
+        #if mode == 'mg5':
+        #    text = 'import aloha.template_files.wavefunctions as wavefunctions\n'
+        #else:
+        #    text = 'import wavefunctions\n'
+                    
+                 
+        # write header 
+        header = self.define_header(name=name)
+        new_couplings = ['COUP%s' % (i+1) for i in range(len(lor_names)+1)]
+        header = header.replace('COUP', ','.join(new_couplings))
+        
+        text += header
+  
+        # Define which part of the routine should be called
+        addon = ''
+        if 'C' in self.namestring:
+            short_name, addon = name.split('C',1)
+            if addon.split('_')[0].isdigit():
+                addon = 'C' +self.namestring.split('C',1)[1]
+            else:
+                addon = '_%s' % self.offshell
+        else:
+            addon = '_%s' % self.offshell
+
+        # how to call the routine
+        if not offshell:
+            main = 'vertex'
+            call_arg = '%(args)s, %(COUP)s' % \
+                    {'args': ', '.join(self.calllist['CallList']), 
+                     'COUP':'COUP%d',
+                     'spin': self.particles[self.offshell -1]}
+        else:
+            main = '%(spin)s%(id)d' % \
+                          {'spin': self.particles[self.offshell -1],
+                           'id': self.offshell}
+            call_arg = '%(args)s, %(COUP)s, M%(id)d, W%(id)d' % \
+                    {'args': ', '.join(self.calllist['CallList']), 
+                     'COUP':'COUP%d',
+                     'id': self.offshell}
+
+        # make the first call
+        line = "    %s = %s%s("+call_arg+")\n"
+        text += '\n\n' + line % (main, self.namestring, '', 1)
+        
+        # make the other call
+        for i,name in enumerate(lor_names):
+            text += line % ('tmp',name, addon, i+2)
+            if not offshell:
+               text += '    vertex += tmp\n'
+            else:
+                size = self.type_to_size[self.particles[offshell -1]] -2
+                text += "    for i in range(%(id)d):\n" % {'id': size}
+                text += "        %(main)s[i] += tmp[i]\n" %{'main': main}
+        
+        text += '    '+self.define_foot()
+
+        #ADD SYMETRY
+        if sym:
+            for elem in self.symmetries:
+                text += self.write_combined(lor_names, mode, elem)
+            
+        return text
+        
\ No newline at end of file

=== modified file 'aloha/bin/aloha'
--- aloha/bin/aloha	2011-02-03 09:33:24 +0000
+++ aloha/bin/aloha	2011-04-22 01:24:30 +0000
@@ -31,7 +31,7 @@
 
 (options, args) = parser.parse_args()
 if len(args) != 1:
-    sys.exit('ALOHA requires as arguments the path to UFO directory')
+    sys.exit('ALOHA requires exactly one argument (the path to UFO directory)')
     
 # Set logging level according to the logging level given by options
 logging.basicConfig(level=vars(logging)[options.logging])

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2011-04-14 21:57:13 +0000
+++ aloha/create_aloha.py	2011-04-22 01:24:30 +0000
@@ -56,18 +56,32 @@
         self.outgoing = outgoing
         self.infostr = infostr
         self.symmetries = []
+        self.combined = []
         
     def add_symmetry(self, outgoing):
         """ add an outgoing """
         
         if not outgoing in self.symmetries:
             self.symmetries.append(outgoing)
-        
-    def write(self, output_dir, language='Fortran', mode='self'):
+    
+    def add_combine(self, lor_list):
+        """add a combine rule """
+        
+        if lor_list not in self.combined:
+            self.combined.append(lor_list)
+        
+    def write(self, output_dir, language='Fortran', mode='self', **opt):
         """ write the content of the object """
-        return getattr(aloha_writers, 'ALOHAWriterFor%s' % language)(self, output_dir).write(mode=mode)
-
-
+        
+        writer = getattr(aloha_writers, 'ALOHAWriterFor%s' % language)(self, output_dir)
+        text = writer.write(mode=mode, **opt)
+        for grouped in self.combined:
+            if isinstance(text, tuple):
+                text = tuple([old.__add__(new)  for old, new in zip(text, 
+                             writer.write_combined(grouped, mode=mode, **opt))])
+            else:
+                text += writer.write_combined(grouped, mode=mode, **opt)
+        return text
 
 class AbstractRoutineBuilder(object):
     """ Launch the creation of the Helicity Routine"""
@@ -346,6 +360,7 @@
         # init the dictionary
         dict.__init__(self)
         self.symmetries = {}
+        self.multiple_lor = {}
         
         # check the mass of spin2 (if any)
         self.massless_spin2 = self.has_massless_spin2()
@@ -394,6 +409,7 @@
         
     def load(self, filepos=None):
         """ reload the pickle file """
+        return False
         if not filepos:
             filepos = os.path.join(self.model_pos,'aloha.pkl') 
         if os.path.exists(filepos):
@@ -426,6 +442,7 @@
         #to compute identical contribution
         self.look_for_symmetries()
         conjugate_list = self.look_for_conjugate()
+        self.look_for_multiple_lorentz_interactions()
         if not wanted_lorentz:
             wanted_lorentz = [l.name for l in self.model.all_lorentz]
         for lorentz in self.model.all_lorentz:
@@ -445,6 +462,13 @@
             if 5 in lorentz.spins and self.massless_spin2 is not None:
                 builder.spin2_massless = self.massless_spin2 
             self.compute_aloha(builder)
+            
+            if lorentz.name in self.multiple_lor:
+                for m in self.multiple_lor[lorentz.name]:
+                    for outgoing in range(len(lorentz.spins)+1):
+                        self[(lorentz.name, outgoing)].add_combine(m)
+                        
+                    
             if lorentz.name in conjugate_list:
                 conjg_builder_list= builder.define_all_conjugate_builder(\
                                                    conjugate_list[lorentz.name])
@@ -452,6 +476,11 @@
                     # No duplication of conjugation:
                     assert conjg_builder_list.count(conjg_builder) == 1
                     self.compute_aloha(conjg_builder, lorentz.name)
+                    if lorentz.name in self.multiple_lor:
+                        for m in self.multiple_lor[lorentz.name]:
+                            for outgoing in range(len(lorentz.spins)+1):
+                                self[(conjg_builder.name, outgoing)].add_combine(m)
+                                                
                     
                     
         if save:
@@ -471,14 +500,16 @@
         # reorganize the data (in order to use optimization for a given lorentz
         #structure
         request = {}
-        for l_name, conjugate, outgoing in data:
-            try:
-                request[l_name][conjugate].append(outgoing)
-            except:
+        for list_l_name, conjugate, outgoing in data:
+            for l_name in list_l_name:
                 try:
-                    request[l_name][conjugate] = [outgoing]
+                    request[l_name][conjugate].append(outgoing)
                 except:
-                    request[l_name] = {conjugate: [outgoing]}
+                    try:
+                        request[l_name][conjugate] = [outgoing]
+                    except:
+                        request[l_name] = {conjugate: [outgoing]}
+                        
         # Loop on the structure to build exactly what is request
         for l_name in request:
             lorentz = eval('self.model.lorentz.%s' % l_name)
@@ -504,7 +535,14 @@
                     # Compute routines
                     self.compute_aloha(conjg_builder, symmetry=lorentz.name,
                                         routines=routines)
-            
+        
+        # Build mutiple lorentz call
+        for list_l_name, conjugate, outgoing in data:
+            if len(list_l_name) >1:
+                lorentzname = list_l_name[0]
+                for c in conjugate:
+                    lorentzname += 'C%s' % c
+                self[(lorentzname, outgoing)].add_combine(list_l_name[1:])
                         
     def compute_aloha(self, builder, symmetry=None, routines=None):
         """ define all the AbstractRoutine linked to a given lorentz structure
@@ -521,7 +559,7 @@
         for outgoing in routines:
             symmetric = self.has_symmetries(symmetry, outgoing, valid_output=routines)
             if symmetric:
-                self.get(symmetry, symmetric).add_symmetry(outgoing)
+                self.get(name, symmetric).add_symmetry(outgoing)
             else:
                 wavefunction = builder.compute_routine(outgoing)
                 #Store the information
@@ -590,9 +628,9 @@
             for i, part1 in enumerate(vertex.particles):
                 for j in range(i-1,-1,-1):
                     part2 = vertex.particles[j]
-                    if part1.name == part2.name and \
-                                        part1.color == part2.color == 1 and\
-                                        part1.spin != 2:
+                    if part1.pdg_code == part2.pdg_code:
+                        if part1.spin == 2 and (i % 2 != j % 2 ):
+                            continue 
                         for lorentz in vertex.lorentz:
                             if self.symmetries.has_key(lorentz.name):
                                 if self.symmetries[lorentz.name].has_key(i+1):
@@ -602,7 +640,46 @@
                             else:
                                 self.symmetries[lorentz.name] = {i+1:j+1}
                         break
+    
+    def look_for_multiple_lorentz_interactions(self):
+        """Search the interaction associate with more than one lorentz structure.
+        If those lorentz structure have the same order and the same color then
+        associate a multiple lorentz routines to ALOHA """
+        
+        orders = {}
+        for coup in self.model.all_couplings:
+            orders[coup.name] = str(coup.order)
+        
+        for vertex in self.model.all_vertices:
+            if len(vertex.lorentz) == 1:
+                continue
+            #remove ghost
+            if -1 in vertex.lorentz[0].spins:
+                continue
+            
+            # assign each order/color to a set of lorentz routine
+            combine = {}
+            for (id_col, id_lor), coup in vertex.couplings.items():
+                order = orders[coup.name]
+                key = (id_col, order)
+                if key in combine:
+                    combine[key].append(id_lor)
+                else:
+                    combine[key] = [id_lor]
                     
+            # Check if more than one routine are associated
+            for list_lor in combine.values():
+                if len(list_lor) == 1:
+                    continue
+                list_lor.sort() 
+                main = vertex.lorentz[list_lor[0]].name 
+                if main not in self.multiple_lor:
+                    self.multiple_lor[main] = []
+                
+                info = tuple([vertex.lorentz[id].name for id in list_lor[1:]])
+                if info not in self.multiple_lor[main]:
+                    self.multiple_lor[main].append(info)
+                
     def has_massless_spin2(self):
         """Search if the spin2 particles are massless or not"""
         
@@ -691,10 +768,10 @@
     aloha_files = []
     
     # Identify the valid files
-    alohafile_pattern = re.compile(r'''^[STFV]*[_C\d]*_\d%s''' % file_ext)
+    alohafile_pattern = re.compile(r'''_\d%s''' % file_ext)
     for filename in os.listdir(aloha_dir):
         if os.path.isfile(os.path.join(aloha_dir, filename)):
-            if alohafile_pattern.match(filename):
+            if alohafile_pattern.search(filename):
                 aloha_files.append(filename.replace(file_ext, comp_ext))
 
     text="ALOHARoutine = "

=== modified file 'input/mg5_configuration.txt'
--- input/mg5_configuration.txt	2011-04-01 09:24:18 +0000
+++ input/mg5_configuration.txt	2011-04-22 01:24:30 +0000
@@ -14,6 +14,7 @@
 ################################################################################
 #
 # This File contains some configuration variable for mg5
+# This File is use if the file ~/.mg5_config is not present.
 #
 ################################################################################
 
@@ -21,3 +22,20 @@
 # Define the path to the pythia8 directory.
 # relative path start from mg5 directory
 pythia8_path = ./pythia8
+# symmetry_max_time
+# Set maximum time for trying to find symmetric diagrams for a single
+# matrix element. 0 means no timeout.
+symmetry_max_time= 600
+
+# Prefered Text Editor
+#  Default: use the shell default Editor
+#           or try to find one available on the system
+text_editor = None
+
+# Prefered WebBrower
+# If None: try to find one available on the system
+web_browser = None
+
+# Prefered PS viewer
+# If None: try to find one available on the system
+eps_viewer = None

=== modified file 'madgraph/VERSION'
--- madgraph/VERSION	2011-04-12 16:49:18 +0000
+++ madgraph/VERSION	2011-04-22 01:24:30 +0000
@@ -1,2 +1,2 @@
-version = 1.0.0
-date = 2011-04-12
+version = 1.1.0
+date = 2011-04-22

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2011-04-03 22:06:00 +0000
+++ madgraph/core/helas_objects.py	2011-04-22 01:24:30 +0000
@@ -72,16 +72,16 @@
         # interaction_id = the id of the interaction in the model
         # pdg_codes = the pdg_codes property of the interaction, [11, -11, 22]
         # inter_color = the 'color' property of the interaction: []
-        # lorentz = the 'lorentz' property of the interaction: ['']
+        # lorentz = the 'lorentz' property of the interaction: ('')
         # couplings = the coupling names from the interaction: {(0,0):'MGVX12'}
         self['interaction_id'] = 0
         self['pdg_codes'] = []
         self['orders'] = {}
         self['inter_color'] = None
-        self['lorentz'] = ''
-        self['coupling'] = 'none'
-        # The Lorentz and color index used in this wavefunction
-        self['coupl_key'] = (0, 0)
+        self['lorentz'] = []
+        self['coupling'] = ['none']
+        # The color index used in this wavefunction
+        self['color_key'] = 0
         # Properties relating to the leg/vertex
         # state = initial/final (for external bosons),
         #         intermediate (for intermediate bosons),
@@ -206,28 +206,32 @@
                             "%s is not a valid Color String" % str(value)
 
         if name == 'lorentz':
-            #Should be a string
-            if not isinstance(value, str):
-                    raise self.PhysicsObjectError, \
-                        "%s is not a valid string" % str(value)
+            #Should be a list of string
+            if not isinstance(value, list):
+                    raise self.PhysicsObjectError, \
+                        "%s is not a valid tuple" % str(value)
+            for name in value:
+                if not isinstance(name, str):
+                    raise self.PhysicsObjectError, \
+                        "%s doesn't contain only string" % str(value)
 
         if name == 'coupling':
-            #Should be a string
-            if not isinstance(value, str):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid coupling string" % \
-                                                                str(value)
+            #Should be a list of string
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid coupling string" % str(value)
+            for name in value:
+                if not isinstance(name, str):
+                    raise self.PhysicsObjectError, \
+                        "%s doesn't contain only string" % str(value)
+            if len(value) == 0:
+                raise self.PhysicsObjectError, \
+                        "%s should have at least one value" % str(value)
 
-        if name == 'coupl_key':
-            if value and not isinstance(value, tuple):
-                raise self.PhysicsObjectError, \
-                      "%s is not a valid tuple" % str(value)
-            if len(value) != 2:
-                raise self.PhysicsObjectError, \
-                      "%s is not a valid tuple with 2 elements" % str(value)
-            if not isinstance(value[0], int) or not isinstance(value[1], int):
-                raise self.PhysicsObjectError, \
-                      "%s is not a valid tuple of integer" % str(value)
+        if name == 'color_key':
+            if value and not isinstance(value, int):
+                raise self.PhysicsObjectError, \
+                      "%s is not a valid integer" % str(value)
 
         if name == 'state':
             if not isinstance(value, str):
@@ -320,9 +324,9 @@
                     if inter.get('color'):
                         self.set('inter_color', inter.get('color')[0])
                     if inter.get('lorentz'):
-                        self.set('lorentz', inter.get('lorentz')[0])
+                        self.set('lorentz', [inter.get('lorentz')[0]])
                     if inter.get('couplings'):
-                        self.set('coupling', inter.get('couplings').values()[0])
+                        self.set('coupling', [inter.get('couplings').values()[0]])
                 return True
             elif name == 'particle':
                 self.set('particle', model.get('particle_dict')[value])
@@ -343,7 +347,7 @@
 
         return ['particle', 'antiparticle', 'is_part',
                 'interaction_id', 'pdg_codes', 'orders', 'inter_color', 
-                'lorentz', 'coupling', 'coupl_key', 'state', 'number_external',
+                'lorentz', 'coupling', 'color_key', 'state', 'number_external',
                 'number', 'fermionflow', 'mothers']
 
     # Helper functions
@@ -369,8 +373,8 @@
         array_rep = array.array('i', [self['interaction_id']])
         # Need the coupling key, to distinguish between
         # wavefunctions from the same interaction but different
-        # Lorentz or color structures
-        array_rep.extend(list(self['coupl_key']))
+        # color structures
+        array_rep.append(self['color_key'])
         # Finally, the mother numbers
         array_rep.extend([mother['number'] for \
                           mother in self['mothers']])
@@ -406,7 +410,7 @@
                    particles[1].get_anti_pdg_code()\
                    and self.get('coupling')[0] != '-':
                 # We need a minus sign in front of the coupling
-                self.set('coupling', '-' + self.get('coupling'))
+                self.set('coupling', ['-' + c for c in self.get('coupling')])
 
     def set_octet_majorana_coupling_sign(self):
         """For octet Majorana fermions, need an extra minus sign in
@@ -418,7 +422,7 @@
                self.get_spin_state_number() == -2 and \
                self.get('self_antipart') and \
                [m.get('color') for m in self.get('mothers')] == [8, 8]:
-            self.set('coupling', '-' + self.get('coupling'))
+            self.set('coupling', ['-' + c for c in self.get('coupling')])
         
     def set_state_and_particle(self, model):
         """Set incoming/outgoing state according to mother states and
@@ -794,7 +798,7 @@
         if self.needs_hermitian_conjugate():
             res.append(self.get_conjugate_index())
 
-        return (tuple(res), self.get('lorentz'))
+        return (tuple(res), tuple(self.get('lorentz')))
 
     def get_base_vertices(self, wf_dict = {}, vx_list = [], optimization = 1):
         """Recursive method to get a base_objects.VertexList
@@ -876,7 +880,7 @@
             # This is where recursion happens
             color_indices.extend(mother.get_color_indices())
         # Add this wf's color index
-        color_indices.append(self.get('coupl_key')[0])
+        color_indices.append(self.get('color_key'))
 
         return color_indices
 
@@ -1058,7 +1062,7 @@
         # Check relevant directly defined properties
         if self['number_external'] != other['number_external'] or \
            self['fermionflow'] != other['fermionflow'] or \
-           self['coupl_key'] != other['coupl_key'] or \
+           self['color_key'] != other['color_key'] or \
            self['lorentz'] != other['lorentz'] or \
            self['coupling'] != other['coupling'] or \
            self['state'] != other['state'] or \
@@ -1293,10 +1297,10 @@
         self['pdg_codes'] = []
         self['orders'] = {}
         self['inter_color'] = None
-        self['lorentz'] = ''
-        self['coupling'] = 'none'
+        self['lorentz'] = []
+        self['coupling'] = ['none']
         # The Lorentz and color index used in this amplitude
-        self['coupl_key'] = (0, 0)
+        self['color_key'] = 0
         # Properties relating to the vertex
         self['number'] = 0
         self['fermionfactor'] = 0
@@ -1359,28 +1363,33 @@
                             "%s is not a valid Color String" % str(value)
 
         if name == 'lorentz':
-            #Should be a string
-            if not isinstance(value, str):
-                    raise self.PhysicsObjectError, \
-                        "%s is not a valid string" % str(value)
-
+            #Should be a list of string
+            if not isinstance(value, list):
+                    raise self.PhysicsObjectError, \
+                        "%s is not a valid list of string" % str(value)
+            for name in value:
+                if not isinstance(name, str):
+                    raise self.PhysicsObjectError, \
+                        "%s doesn't contain only string" % str(value)
+                        
         if name == 'coupling':
-            #Should be a string
-            if not isinstance(value, str):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid coupling string" % \
-                                                                str(value)
+            #Should be a list of string
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError, \
+                      "%s is not a valid coupling (list of string)" % str(value)
+            
+            for name in value:
+                if not isinstance(name, str):
+                    raise self.PhysicsObjectError, \
+                        "%s doesn't contain only string" % str(value)
+            if not len(value):
+                raise self.PhysicsObjectError, \
+                                      'coupling should have at least one value'
 
-        if name == 'coupl_key':
-            if value and not isinstance(value, tuple):
-                raise self.PhysicsObjectError, \
-                      "%s is not a valid tuple" % str(value)
-            if len(value) != 2:
-                raise self.PhysicsObjectError, \
-                      "%s is not a valid tuple with 2 elements" % str(value)
-            if not isinstance(value[0], int) or not isinstance(value[1], int):
-                raise self.PhysicsObjectError, \
-                      "%s is not a valid tuple of integer" % str(value)
+        if name == 'color_key':
+            if value and not isinstance(value, int):
+                raise self.PhysicsObjectError, \
+                      "%s is not a valid integer" % str(value)
 
         if name == 'number':
             if not isinstance(value, int):
@@ -1416,6 +1425,30 @@
 
         return True
 
+    def __str__(self):
+        """ practicle way to represent an HelasAmplitude"""
+        
+        mystr = '{\n'
+        for prop in self.get_sorted_keys():
+            if isinstance(self[prop], str):
+                mystr = mystr + '    \'' + prop + '\': \'' + \
+                        self[prop] + '\',\n'
+            elif isinstance(self[prop], float):
+                mystr = mystr + '    \'' + prop + '\': %.2f,\n' % self[prop]
+            elif isinstance(self[prop], int):
+                mystr = mystr + '    \'' + prop + '\': %s,\n' % self[prop]                
+            elif prop != 'mothers':
+                mystr = mystr + '    \'' + prop + '\': ' + \
+                       str(self[prop]) + ',\n'
+            else:
+                info = [m.get('pdg_code') for m in self['mothers']]
+                mystr += '    \'%s\': %s,\n' % (prop, info) 
+                
+        mystr = mystr.rstrip(',\n')
+        mystr = mystr + '\n}'
+
+        return mystr
+
     # Enhanced get function
     def get(self, name):
         """Get the value of the property name."""
@@ -1453,9 +1486,9 @@
                     if inter.get('color'):
                         self.set('inter_color', inter.get('color')[0])
                     if inter.get('lorentz'):
-                        self.set('lorentz', inter.get('lorentz')[0])
+                        self.set('lorentz', [inter.get('lorentz')[0]])
                     if inter.get('couplings'):
-                        self.set('coupling', inter.get('couplings').values()[0])
+                        self.set('coupling', [inter.get('couplings').values()[0]])
                 return True
             else:
                 raise self.PhysicsObjectError, \
@@ -1467,7 +1500,7 @@
         """Return particle property names as a nicely sorted list."""
 
         return ['interaction_id', 'pdg_codes', 'orders', 'inter_color', 
-                'lorentz', 'coupling', 'coupl_key', 'number', 'color_indices',
+                'lorentz', 'coupling', 'color_key', 'number', 'color_indices',
                 'fermionfactor', 'mothers']
 
 
@@ -1516,7 +1549,7 @@
         if self.needs_hermitian_conjugate():
             res.append(self.get_conjugate_index())
 
-        return (tuple(res), self.get('lorentz'))
+        return (tuple(res), tuple(self.get('lorentz')))
 
     def calculate_fermionfactor(self):
         """Calculate the fermion factor for the diagram corresponding
@@ -1770,7 +1803,7 @@
 
         # Add this amp's color index
         if self.get('interaction_id'):
-            color_indices.append(self.get('coupl_key')[0])
+            color_indices.append(self.get('color_key'))
 
         return color_indices
 
@@ -2218,15 +2251,22 @@
 
                     # Now generate new wavefunction for the last leg
 
-                    # Need one amplitude for each Lorentz/color structure,
-                    # i.e. for each coupling
+                    # Need one amplitude for each color structure,
+                    done_color = {} # store link to color
                     for coupl_key in sorted(inter.get('couplings').keys()):
+                        color = coupl_key[0]
+                        if color in done_color:
+                            wf = done_color[color]
+                            wf.get('coupling').append(inter.get('couplings')[coupl_key])
+                            wf.get('lorentz').append(inter.get('lorentz')[coupl_key[1]])
+                            continue
                         wf = HelasWavefunction(last_leg, vertex.get('id'), model)
-                        wf.set('coupling', inter.get('couplings')[coupl_key])
+                        wf.set('coupling', [inter.get('couplings')[coupl_key]])
                         if inter.get('color'):
                             wf.set('inter_color', inter.get('color')[coupl_key[0]])
-                        wf.set('lorentz', inter.get('lorentz')[coupl_key[1]])
-                        wf.set('coupl_key', coupl_key)
+                        done_color[color] = wf
+                        wf.set('lorentz', [inter.get('lorentz')[coupl_key[1]]])
+                        wf.set('color_key', color)
                         wf.set('mothers', mothers)
                         # Need to set incoming/outgoing and
                         # particle/antiparticle according to the fermion flow
@@ -2312,24 +2352,29 @@
                                               wf_number,
                                               False,
                                               number_to_wavefunctions)
-
+                done_color = {}
                 for i, coupl_key in enumerate(keys):
+                    color = coupl_key[0]
+                    if inter and color in done_color.keys():
+                        amp = done_color[color]
+                        amp.get('coupling').append(inter.get('couplings')[coupl_key])
+                        amp.get('lorentz').append(inter.get('lorentz')[coupl_key[1]])
+                        continue
                     amp = HelasAmplitude(lastvx, model)
                     if inter:
-                        amp.set('coupling', inter.get('couplings')[coupl_key])
-                        amp.set('lorentz', inter.get('lorentz')[\
-                                coupl_key[1]])
+                        amp.set('coupling', [inter.get('couplings')[coupl_key]])
+                        amp.set('lorentz', [inter.get('lorentz')[coupl_key[1]]])
                         if inter.get('color'):
-                            amp.set('inter_color', inter.get('color')[\
-                                coupl_key[0]])
-                        amp.set('coupl_key', coupl_key)
+                            amp.set('inter_color', inter.get('color')[color])
+                        amp.set('color_key', color)
+                        done_color[color] = amp
                     amp.set('mothers', mothers)
                     amplitude_number = amplitude_number + 1
                     amp.set('number', amplitude_number)
                     # Add the list with color indices to the amplitude
                     new_color_list = copy.copy(color_list)
                     if inter:
-                        new_color_list.append(coupl_key[0])
+                        new_color_list.append(color)
                         
                     amp.set('color_indices', new_color_list)
 
@@ -3291,7 +3336,7 @@
         """Return a list of (lorentz_name, conjugate, outgoing) with
         all lorentz structures used by this HelasMatrixElement."""
 
-        return [(wa.get('lorentz'), tuple(wa.get_conjugate_index()),
+        return [(tuple(wa.get('lorentz')), tuple(wa.get_conjugate_index()),
                  wa.find_outgoing_number()) for wa in \
                 self.get_all_wavefunctions() + self.get_all_amplitudes() \
                 if wa.get('interaction_id') != 0]
@@ -3967,7 +4012,7 @@
         coupling_list = []
 
         for me in self.get('matrix_elements'):
-            coupling_list.extend(me.get_used_couplings())
+            coupling_list.extend([c for l in me.get_used_couplings() for c in l])
 
         return list(set(coupling_list))
     

=== modified file 'madgraph/interface/cmd_interface.py'
--- madgraph/interface/cmd_interface.py	2011-04-11 19:58:39 +0000
+++ madgraph/interface/cmd_interface.py	2011-04-22 01:24:30 +0000
@@ -12,6 +12,7 @@
 # For more information, please visit: http://madgraph.phys.ucl.ac.be
 #
 ################################################################################
+from compiler.ast import Break
 """A user friendly command line interface to access MadGraph features.
    Uses the cmd package for command interpretation and tab completion.
 """
@@ -84,10 +85,10 @@
                                                     'display multiparticles'],
         'generate': ['add process PROCESS','output [OUTPUT_TYPE] [PATH]','draw .'],
         'add process':['output [OUTPUT_TYPE] [PATH]', 'display processes'],
-        'output':['launch','history PATH', 'exit'],
+        'output':['launch','open index.html','history PATH', 'exit'],
         'display': ['generate PROCESS', 'add process PROCESS', 'output [OUTPUT_TYPE] [PATH]'],
-        'draw': ['shell CMD'],
         'import proc_v4' : ['launch','exit'],
+        'launch': ['open index.html','exit'],
         'tutorial': ['generate PROCESS', 'import model MODEL', 'help TOPIC']
     }
     
@@ -280,9 +281,11 @@
         logger.info("-- display a the status of various internal state variables")
         logger.info("   for particles/interactions you can specify the name or id of the")
         logger.info("   particles/interactions to receive more details information.")
-        logger.info("   example display particles e+.")
+        logger.info("   Example: display particles e+.")
         logger.info("   For \"checks\", can specify only to see failed checks.")
-
+        logger.info("   For \"diagrams\", you can specify where the file will be written.")
+        logger.info("   Example: display diagrams ./")
+        
     def help_launch(self):
         """help for launch command"""
         _launch_parser.print_help()
@@ -291,8 +294,16 @@
         logger.info("syntax: tutorial [" + "|".join(self._tutorial_opts) + "]")
         logger.info("-- start/stop the tutorial mode")
 
+    def help_open(self):
+        logger.info("syntax: open FILE  ")
+        logger.info("-- open a file with the appropriate editor.")
+        logger.info('   If FILE belongs to index.html, param_card.dat, run_card.dat')
+        logger.info('   the path to the last created/used directory is used')
+        logger.info('   The program used to open those files can be chosen in the')
+        logger.info('   configuration file ./input/mg5_configuration.txt')
+        
     def help_output(self):
-        logger.info("syntax [" + "|".join(self._export_formats) + \
+        logger.info("syntax: output [" + "|".join(self._export_formats) + \
                     "] [path|.|auto] [options]")
         logger.info("-- Output any generated process(es) to file.")
         logger.info("   mode: Default mode is madevent. Default path is \'.\' or auto.")
@@ -380,13 +391,13 @@
         logger.info("   If FILEPATH is omitted, the history will be output to stdout.")
         logger.info("   \"clean\" will remove all entries from the history.")
 
-    def help_draw(self):
-        _draw_parser.print_help()
+#    def help_draw(self):
+#        _draw_parser.print_help()
 
     def help_set(self):
         logger.info("syntax: set %s argument" % "|".join(self._set_options))
         logger.info("-- set options for generation or output")
-        logger.info("   group_subprocesses_output True/False: ")
+        logger.info("   group_subprocesses True/False: ")
         logger.info("     Smart grouping of subprocesses into directories,")
         logger.info("     mirroring of initial states, and combination of")
         logger.info("     integration channels.")
@@ -397,7 +408,10 @@
         logger.info("     of the quarks given in multi_part_label.")
         logger.info("     These processes give negligible contribution to the")
         logger.info("     cross section but have subprocesses/channels.")
-
+        logger.info("   symmetry_max_time N")
+        logger.info("     (default 600) maximum time (in s) to find symmetric")
+        logger.info("     diagrams for each matrix element (0 means no timeout)")
+        
     def help_shell(self):
         logger.info("syntax: shell CMD (or ! CMD)")
         logger.info("-- run the shell command CMD and catch output")
@@ -491,8 +505,7 @@
         """
         
         if len(args) < 1:
-            self.help_draw()
-            raise self.InvalidCmd('\"draw\" command requires a directory path')
+            args.append('/tmp')
         
         if not self._curr_amps:
             raise self.InvalidCmd("No process generated, please generate a process!")
@@ -665,6 +678,8 @@
         mode = self.find_output_type(path)
         args[0] = mode
         args.append(path)
+        # inform where we are for future command
+        self._done_export = [path, mode]
         
     
     def find_output_type(self, path):
@@ -714,7 +729,7 @@
             raise self.InvalidCmd('Possible options for set are %s' % \
                                   self._set_options)
 
-        if args[0] in ['group_subprocesses_output']:
+        if args[0] in ['group_subprocesses']:
             if args[1] not in ['False', 'True']:
                 raise self.InvalidCmd('%s needs argument False or True' % \
                                       args[0])
@@ -727,7 +742,46 @@
             if args[1] not in ['DEBUG','INFO','WARNING','ERROR','CRITICAL']:
                 raise self.InvalidCmd('output_level needs ' + \
                                       'a valid level')       
+    
+    def check_open(self, args):
+        """ check the validity of the line """
         
+        if len(args) != 1:
+            self.help_open()
+            raise self.InvalidCmd('OPEN command requires exactly one argument')
+
+        if args[0].startswith('./'):
+            if not os.path.isfile(args[0]):
+                raise self.InvalidCmd('%s: not such file' % args[0])
+            return True
+
+        # if special : create the path.
+        if not self._done_export:
+            if not os.path.isfile(args[0]):
+                self.help_open()
+                raise self.InvalidCmd('No command \"output\" or \"launch\" used. Impossible to associate this name to a file')
+            else:
+                return True
+            
+        path = self._done_export[0]
+        if os.path.isfile(os.path.join(path,args[0])):
+            args[0] = os.path.join(path,args[0])
+        elif os.path.isfile(os.path.join(path,'Cards',args[0])):
+            args[0] = os.path.join(path,'Cards',args[0])
+        elif os.path.isfile(os.path.join(path,'HTML',args[0])):
+            args[0] = os.path.join(path,'HTML',args[0])
+        # special for card with _default define: copy the default and open it
+        elif '_card.dat' in args[0]:   
+            name = args[0].replace('_card.dat','_card_default.dat')
+            if os.path.isfile(os.path.join(path,'Cards', name)):
+                files.cp(path + '/Cards/' + name, path + '/Cards/'+ args[0])
+                args[0] = os.path.join(path,'Cards', args[0])
+            else:
+                raise self.InvalidCmd('No default path for this file')
+        elif not os.path.isfile(args[0]):
+            raise self.InvalidCmd('No default path for this file')
+                
+                
     def check_output(self, args):
         """ check the validity of the line"""
         
@@ -746,7 +800,7 @@
 
         if self._model_v4_path and \
                (self._export_format not in self._v4_export_formats or \
-                self._options['group_subprocesses_output']):
+                self._options['group_subprocesses']):
             text = " The Model imported (MG4 format) does not contain enough\n "
             text += " information for this type of output. In order to create\n"
             text += " output for " + args[0] + ", you have to use a UFO model.\n"
@@ -763,30 +817,37 @@
                 self.get_default_path()
             else:
                 self._export_dir = path
-
-        if not self._export_dir:
+        else:
             # No valid path
             self.get_default_path()
 
-        if self._done_export == (self._export_dir, self._export_format):
-            # We have already done export in this path
-            logger.info("Matrix elements already exported to directory %s" % \
-                                    self._export_dir)
-            return        
-
-        if self._export_format in ['madevent', 'standalone'] \
-               and not self._mgme_dir and \
-               os.path.realpath(self._export_dir) != os.path.realpath('.'):
-            raise MadGraph5Error, \
-                  "To generate a new MG4 directory, you need a valid MG_ME path"
-
         self._export_dir = os.path.realpath(self._export_dir)
 
             
     def get_default_path(self):
         """Set self._export_dir to the default (\'auto\') path"""
+
+        if self._export_format in ['madevent', 'standalone']:
+            # Detect if this script is launched from a valid copy of the Template,
+            # if so store this position as standard output directory
+            if 'TemplateVersion.txt' in os.listdir('.'):
+                #Check for ./
+                self._export_dir = os.path.realpath('.')
+                return
+            elif 'TemplateVersion.txt' in os.listdir('..'):
+                #Check for ../
+                self._export_dir = os.path.realpath('..')
+                return
+            elif self.stdin != sys.stdin:
+                #Check for position defined by the input files
+                input_path = os.path.realpath(self.stdin.name).split(os.path.sep)
+                print "Not standard stdin, use input path"
+                if input_path[-2] == 'Cards':
+                    self._export_dir = os.path.sep.join(input_path[:-2])
+                    if 'TemplateVersion.txt' in self._export_dir:
+                        return
         
-        if self._export_format.startswith('madevent'):
+        if self._export_format.startswith('madevent'):            
             name_dir = lambda i: 'PROC_%s_%s' % \
                                     (self._curr_model['name'], i)
             auto_path = lambda i: os.path.join(self.writing_dir,
@@ -885,7 +946,10 @@
         """ not authorize on web"""
         raise self.WebRestriction('\"save\" command not authorize online')
     
-
+    def check_open(self, args):
+        """ not authorize on web"""
+        raise self.WebRestriction('\"open\" command not authorize online')
+    
     def check_output(self, args):
         """ check the validity of the line"""
 
@@ -1141,7 +1205,37 @@
         # Filename if directory is not given
         if len(args) == 2:
             return self.path_completion(text)
-
+        
+    def complete_open(self, text, line, begidx, endidx): 
+        """ complete the open command """
+
+        args = split_arg(line[0:begidx])
+        
+        # Directory continuation
+        if os.path.sep in args[-1] + text:
+            return self.path_completion(text,
+                                    os.path.join('.',*[a for a in args if \
+                                                      a.endswith(os.path.sep)]))
+
+        possibility = []
+        if self._done_export:
+            path = self._done_export[0]
+            possibility = ['index.html']
+            if os.path.isfile(os.path.join(path,'README')):
+                possibility.append('README')
+            if os.path.isdir(os.path.join(path,'Cards')):
+                possibility += [f for f in os.listdir(os.path.join(path,'Cards')) 
+                                    if f.endswith('.dat')]
+            if os.path.isdir(os.path.join(path,'HTML')):
+                possibility += [f for f in os.listdir(os.path.join(path,'HTML')) 
+                                  if f.endswith('.html') and 'default' not in f]
+        else:
+            possibility.extend(['./','../'])
+        if os.path.exists('MG5_debug'):
+            possibility.append('MG5_debug')
+
+        return self.list_completion(text, possibility)
+       
     def complete_output(self, text, line, begidx, endidx,
                         possible_options = ['f', 'noclean', 'nojpeg'],
                         possible_options_full = ['-f', '-noclean', '-nojpeg']):
@@ -1191,7 +1285,7 @@
             return self.list_completion(text, self._set_options)
 
         if len(args) == 2:
-            if args[1] in ['group_subprocesses_output']:
+            if args[1] in ['group_subprocesses']:
                 return self.list_completion(text, ['False', 'True'])
             
             elif args[1] in ['ignore_six_quark_processes']:
@@ -1340,8 +1434,8 @@
 
     # Options and formats available
     _display_opts = ['particles', 'interactions', 'processes', 'diagrams', 
-                     'multiparticles', 'couplings', 'lorentz', 'checks',
-                     'parameters']
+                     'diagrams_text', 'multiparticles', 'couplings', 'lorentz', 
+                     'checks', 'parameters']
     _add_opts = ['process']
     _save_opts = ['model', 'processes']
     _tutorial_opts = ['start', 'stop']
@@ -1349,8 +1443,9 @@
     _import_formats = ['model_v4', 'model', 'proc_v4', 'command']
     _v4_export_formats = ['madevent', 'standalone', 'matrix'] 
     _export_formats = _v4_export_formats + ['standalone_cpp', 'pythia8']
-    _set_options = ['group_subprocesses_output',
+    _set_options = ['group_subprocesses',
                     'ignore_six_quark_processes',
+                    'symmetry_max_time',
                     'stdout_level']
     # Variables to store object information
     _curr_model = None  #base_objects.Model()
@@ -1391,21 +1486,6 @@
                              mgme_dir)
                 self._mgme_dir = MG4DIR
 
-        # Detect if this script is launched from a valid copy of the Template,
-        # if so store this position as standard output directory
-        if 'TemplateVersion.txt' in os.listdir('.'):
-            #Check for ./
-            self._export_dir = os.path.realpath('.')
-        elif 'TemplateVersion.txt' in os.listdir('..'):
-            #Check for ../
-            self._export_dir = os.path.realpath('..')
-        elif self.stdin != sys.stdin:
-            #Check for position defined by the input files
-            input_path = os.path.realpath(self.stdin.name).split(os.path.sep)
-            print "Not standard stdin, use input path"
-            if input_path[-2] == 'Cards':
-                self._export_dir = os.path.sep.join(input_path[:-2])
-                
         # Load the configuration file
         self.set_configuration()
         
@@ -1447,7 +1527,7 @@
 
             # Generate processes
             collect_mirror_procs = \
-                                 self._options['group_subprocesses_output']
+                                 self._options['group_subprocesses']
             ignore_six_quark_processes = \
                            self._options['ignore_six_quark_processes'] if \
                            "ignore_six_quark_processes" in self._options \
@@ -1513,6 +1593,9 @@
         #check the validity of the arguments
         self.check_display(args)
 
+        if args[0] == 'diagrams':
+            self.draw(' '.join(args[1:]))
+
         if args[0] == 'particles' and len(args) == 1:
             print "Current model contains %i particles:" % \
                     len(self._curr_model['particles'])
@@ -1589,7 +1672,7 @@
             for amp in self._curr_amps:
                 print amp.nice_string_processes()
 
-        elif args[0] == 'diagrams':
+        elif args[0] == 'diagrams_text':
             text = "\n".join([amp.nice_string() for amp in self._curr_amps])
             pydoc.pager(text)
 
@@ -1691,7 +1774,7 @@
                        "please run from a" + \
                        "\n\t         valid MG_ME directory.")
 
-    def do_draw(self, line):
+    def draw(self, line):
         """ draw the Feynman diagram for the given process """
 
         args = split_arg(line)
@@ -1728,6 +1811,7 @@
                          amp.get('process').nice_string())
             plot.draw(opt=options)
             logger.info("Wrote file " + filename)
+            self.exec_cmd('open %s' % filename)
 
         stop = time.time()
         logger.info('time to draw %s' % (stop - start)) 
@@ -1820,6 +1904,9 @@
         # Also reset _export_format and _export_dir
         self._export_format = None
 
+        # Remove previous generations from history
+        self.clean_history('generate', 'add process')
+
         # Call add process
         args = split_arg(line)
         args.insert(0, 'process')
@@ -1958,7 +2045,7 @@
             #                           leg.get('state') == False,
             #                           myleglist)) == 1
 
-            if overall_orders and self._options['group_subprocesses_output']:
+            if overall_orders and self._options['group_subprocesses']:
                 raise MadGraph5Error, \
                       "For grouped subprocess output, orders should be specified for each process (no overall orders after @N)"                
 
@@ -2134,7 +2221,18 @@
 
         if self.log:
             logger.info("History written to " + output_file.name)
-    
+
+    def clean_history(self,*arguments):
+        """Remove all commands in arguments from history"""
+
+        nline = 0
+        arguments = list(arguments) + ['display', 'open', 'launch', 'output']
+        while nline < len(self.history) - 1:
+            if any([self.history[nline].startswith(arg) for arg in arguments]):
+                self.history.pop(nline)
+            else:
+                nline += 1
+
     # Import files
     def do_import(self, line):
         """Import files with external formats"""
@@ -2146,7 +2244,9 @@
         if args[0].startswith('model'):
             self._model_v4_path = None
             # Clear history, amplitudes and matrix elements when a model is imported
-            self.history = self.history[-1:]
+            # Remove previous imports, generations and outputs from history
+            self.clean_history('import', 'generate', 'add process')
+            # Reset amplitudes and matrix elements
             self._curr_amps = diagram_generation.AmplitudeList()
             self._curr_matrix_elements = helas_objects.HelasMultiProcess()
             # Import model
@@ -2157,7 +2257,7 @@
                       helas_call_writers.FortranHelasCallWriter(\
                                                                self._curr_model)
                 # Automatically turn off subprocess grouping
-                self.do_set('group_subprocesses_output False')
+                self.do_set('group_subprocesses False')
             else:
                 self._curr_model = import_ufo.import_model(args[1])
                 self._curr_fortran_model = \
@@ -2167,7 +2267,7 @@
                       helas_call_writers.CPPUFOHelasCallWriter(\
                                                                self._curr_model)
                 # Automatically turn on subprocess grouping
-                self.do_set('group_subprocesses_output True')
+                self.do_set('group_subprocesses True')
 
             if '-modelname' not in args:
                 self._curr_model.pass_particles_name_in_mg_default()
@@ -2181,6 +2281,9 @@
             process_checks.store_aloha = []
             
         elif args[0] == 'command':
+            # Remove previous imports, generations and outputs from history
+            self.clean_history('import', 'generate', 'add process')
+
             if not os.path.isfile(args[1]):
                 raise MadGraph5Error("Path %s is not a valid pathname" % args[1])
             else:
@@ -2192,6 +2295,9 @@
         
         elif args[0] == 'proc_v4':
             
+            # Remove previous imports, generations and outputs from history
+            self.clean_history('import', 'generate', 'add process')
+
             if len(args) == 1 and self._export_dir:
                 proc_card = os.path.join(self._export_dir, 'Cards', \
                                                                 'proc_card.dat')
@@ -2327,7 +2433,10 @@
         """ assign all configuration variable from file 
             ./input/mg5_configuration.txt. assign to default if not define """
             
-        config = {'pythia8_path': './pythia8'}
+        self.configuration = {'pythia8_path': './pythia8',
+                              'web_browser':None,
+                              'eps_viewer':None,
+                              'text_editor':None}
         
         if not config_path:
             try:
@@ -2351,21 +2460,81 @@
             else:
                 name = name.strip()
                 value = value.strip()
-                config[name] = value
+                self.configuration[name] = value
+                if value.lower() == "none":
+                    self.configuration[name] = None
 
         # Treat each expected input
         # 1: Pythia8_path
         # try relative path
-        pythia8_dir = os.path.join(MG5DIR, config['pythia8_path'])
-        if not os.path.isfile(os.path.join(pythia8_dir, 'include', 'Pythia.h')):
-            if os.path.isfile(os.path.join(config['pythia8_path'], 'include', 'Pythia.h')):
-                pythia8_dir = config['pythia8_path']
+        for key in self.configuration:
+            if key == 'pythia8_path':
+                pythia8_dir = os.path.join(MG5DIR, self.configuration['pythia8_path'])
+                if not os.path.isfile(os.path.join(pythia8_dir, 'include', 'Pythia.h')):
+                    if os.path.isfile(os.path.join(self.configuration['pythia8_path'], 'include', 'Pythia.h')):
+                        pythia8_dir = self.configuration['pythia8_path']
+                    else:
+                        pythia8_dir = None
+                self.pythia8_path = pythia8_dir
+            elif key == 'text_editor':
+                # Treat text editor -> overwrites madgraph value
+                if self.configuration[key] and not misc.which(self.configuration[key]):
+                    logger.warning('Specified text editor %s not valid.' % \
+                                   self.configuration[key])
+                    self.configuration[key] = None
+                if self.configuration[key]:
+                    # All is good
+                    pass
+                elif os.environ.has_key('EDITOR'):
+                    self.configuration[key] = os.environ['EDITOR']
+                else:
+                    prog = ['vi', 'emacs', 'vim', 'gedit', 'nano']
+                    for p in prog:
+                        if misc.which(p):
+                            self.configuration[key] = p
+                            logger.warning(('Using default text editor \"%s\". ' % p) + \
+                                       'Set text_editor in ./input/mg5_configuration.txt')
+                            break
+                if not self.configuration[key]:
+                    logger.warning('No valid text editor found. ' + \
+                                   'Please set in ./input/mg5_configuration.txt') 
+            elif key == 'eps_viewer':
+                if self.configuration[key]:
+                    continue
+                prog = ['gv', 'ggv', 'evince']
+                for p in prog:
+                    if misc.which(p):
+                        self.configuration[key] = p
+                        logger.warning(('Using default eps viewer \"%s\". ' % p) + \
+                                    'Set eps_viewer in ./input/mg5_configuration.txt')
+                        break
+                if not self.configuration[key] and sys.platform != 'darwin':
+                    logger.warning('No valid eps viewer found. ' + \
+                                   'Please set in ./input/mg5_configuration.txt') 
+            elif key == 'web_browser':
+                if self.configuration[key]:
+                    continue
+                prog = ['firefox', 'chrome', 'safari','opera']
+                for p in prog:
+                    if misc.which(p):
+                        self.configuration[key] = p
+                        logger.warning(('Using default web browser \"%s\". ' % p) + \
+                                      'Set web_browser in ./input/mg5_configuration.txt')
+                        break
+                if not self.configuration[key] and sys.platform != 'darwin':
+                    logger.warning('No valid web browser found. ' + \
+                                   'Please set in ./input/mg5_configuration.txt') 
             else:
-                pythia8_dir = None
-        self.pythia8_path = pythia8_dir
-        
-        return config
-                
+                # Default: try to set parameter
+                try:
+                    self.do_set("%s %s" % (key, self.configuration[key]))
+                except MadGraph5Error:
+                    logger.warning("Option %s from config file not understood" \
+                                   % key)
+
+          
+        return self.configuration
+     
     def check_for_export_dir(self, filepath):
         """Check if the files is in a valid export directory and assign it to
         export path if if is"""
@@ -2396,16 +2565,28 @@
         # args is now MODE PATH
         
         if args[0].startswith('standalone'):
-            ext_program = launch_ext.SALauncher(args[1], self.timeout, **options)
+            ext_program = launch_ext.SALauncher(args[1], self.timeout,
+                                                configuration = self.configuration,
+                                                **options)
         elif args[0] == 'madevent':
-            ext_program = launch_ext.MELauncher(args[1], self.timeout, **options)
+            #check if this is a cross-section
+            if len(self._generate_info.split('>')[0].strip().split())>1:
+                ext_program = launch_ext.MELauncher(args[1], self.timeout,
+                                                    configuration = self.configuration,
+                                                    **options)
+            else:
+                # This is a width computation
+                ext_program = launch_ext.MELauncher(args[1], self.timeout, unit='GeV',
+                                                    configuration = self.configuration,
+                                                    **options)
         elif args[0] == 'pythia8':
-            ext_program = launch_ext.Pythia8Launcher(args[1], self.timeout, **options)
+            ext_program = launch_ext.Pythia8Launcher(args[1], self.timeout,
+                                                configuration = self.configuration,
+                                                **options)
         else:
             raise self.InvalidCmd , '%s cannot be run from MG5 interface' % args[0]
         
         
-        
         ext_program.run()
         
         
@@ -2478,10 +2659,10 @@
                 raise self.RWError('Could not load processes from file %s' % args[1])
         if self._model_v4_path:
             # Automatically turn off subprocess grouping
-            self.do_set('group_subprocesses_output False')
+            self.do_set('group_subprocesses False')
         else:
             # Automatically turn on subprocess grouping
-            self.do_set('group_subprocesses_output True')
+            self.do_set('group_subprocesses True')
         
     
     def do_save(self, line):
@@ -2527,22 +2708,40 @@
                             self._curr_model.get_particle(q).get('name') \
                             for q in self._options[args[0]]]))
             
-        elif args[0] == 'group_subprocesses_output':
+        elif args[0] == 'group_subprocesses':
             self._options[args[0]] = eval(args[1])
-            logger.info('Set group_subprocesses_output to %s' % \
+            logger.info('Set group_subprocesses to %s' % \
+                        str(self._options[args[0]]))
+            
+        elif args[0] == 'symmetry_max_time':
+            self._options[args[0]] = int(args[1])
+            logger.info('Set symmetry_max_time to %s' % \
                         str(self._options[args[0]]))
             
         elif args[0] == "stdout_level":
             logging.root.setLevel(eval('logging.' + args[1]))
             logging.getLogger('madgraph').setLevel(eval('logging.' + args[1]))
             logger.info('set output information to level: %s' % args[1])
-        
+    
+    def do_open(self, line):
+        """Open a text file/ eps file / html file"""
+        
+        args = split_arg(line)
+        # Check Argument validity and modify argument to be the real path
+        self.check_open(args)
+        file_path = args[0]
+        
+        launch_ext.open_file(file_path, self.configuration)
+                 
     def do_output(self, line):
         """Initialize a new Template or reinitialize one"""
 
         args = split_arg(line)
         # Check Argument validity
         self.check_output(args)
+
+        # Remove previous outputs from history
+        self.clean_history()
         
         noclean = '-noclean' in args
         force = '-f' in args 
@@ -2553,11 +2752,6 @@
         except:
             pass
             
-        if self._done_export == (self._export_dir, self._export_format):
-            logger.info('Matrix elements already exported to directory %s' % \
-                        self._export_dir)
-            return
-
         if not force and not noclean and os.path.isdir(self._export_dir)\
                and self._export_format in ['madevent', 'standalone']:
             # Don't ask if user already specified force or noclean
@@ -2568,13 +2762,14 @@
                 raise MadGraph5Error('Stopped by user request')
 
         group_subprocesses = self._export_format == 'madevent' and \
-                             self._options['group_subprocesses_output']
+                             self._options['group_subprocesses']
         # Make a Template Copy
         if self._export_format == 'madevent':
             if group_subprocesses:
                 self._curr_exporter = export_v4.ProcessExporterFortranMEGroup(\
                                       self._mgme_dir, self._export_dir,
-                                      not noclean)
+                                      not noclean,
+                                      self._options["symmetry_max_time"])
             else:
                 self._curr_exporter = export_v4.ProcessExporterFortranME(\
                                       self._mgme_dir, self._export_dir,
@@ -2596,6 +2791,12 @@
         # Perform export and finalize right away
         self.export(nojpeg, main_file_name)
 
+        # Automatically run finalize
+        self.finalize(nojpeg)
+            
+        # Remember that we have done export
+        self._done_export = (self._export_dir, self._export_format)
+
         # Reset _export_dir, so we don't overwrite by mistake later
         self._export_dir = None
 
@@ -2615,7 +2816,7 @@
             cpu_time1 = time.time()
             ndiags = 0
             if not self._curr_matrix_elements.get_matrix_elements():
-                if self._options['group_subprocesses_output']:
+                if self._options['group_subprocesses']:
                     cpu_time1 = time.time()
                     dc_amps = [amp for amp in self._curr_amps if isinstance(amp, \
                                         diagram_generation.DecayChainAmplitude)]
@@ -2780,12 +2981,12 @@
                 
         if self._export_format == 'pythia8':
             logger.info("- All necessary files for Pythia 8 generated.")
-            logger.info("  Please go to %s/examples and run" % path)
+            logger.info("- Run \"launch\" and select %s.cc," % filename)
+            logger.info("  or go to %s/examples and run" % path)
             logger.info("      make -f %s" % make_filename)
             logger.info("  (with process_name replaced by process name).")
             logger.info("  You can then run ./%s to produce events for the process" % \
                         filename)
-            logger.info("- Or run launch and select %s.cc." % filename)
 
         # Replace the amplitudes with the actual amplitudes from the
         # matrix elements, which allows proper diagram drawing also of
@@ -2794,12 +2995,6 @@
                [me.get('base_amplitude') for me in \
                 matrix_elements])
 
-        # Remember that we have done export
-        self._done_export = (self._export_dir, self._export_format)
-
-        # Automatically run finalize
-        self.finalize(nojpeg)
-            
     def finalize(self, nojpeg, online = False):
         """Make the html output, write proc_card_mg5.dat and create
         madevent.tar.gz for a MadEvent directory"""
@@ -2844,12 +3039,9 @@
         if self._export_format in ['madevent', 'standalone', 'standalone_cpp']:
             logger.info('Output to directory ' + self._export_dir + ' done.')
         if self._export_format == 'madevent':
-            logger.info('Please see ' + self._export_dir + '/README')
-            logger.info('for information about how to generate events from this process.')
-            logger.info('You can also use the launch command.')
-
-        #reinitialize to empty the default output dir
-        self._export_dir = None
+            logger.info('Type \"launch\" to generate events from this process, or see')
+            logger.info(self._export_dir + '/README')
+            logger.info('Run \"open index.html\" to see more information about this process.')
 
     def do_help(self, line):
         """ propose some usefull possible action """

=== modified file 'madgraph/interface/launch_ext_program.py'
--- madgraph/interface/launch_ext_program.py	2011-04-04 09:45:24 +0000
+++ madgraph/interface/launch_ext_program.py	2011-04-22 01:24:30 +0000
@@ -18,6 +18,7 @@
 import os
 import pydoc
 import re
+import sys
 import subprocess
 import time
 
@@ -38,40 +39,20 @@
     cards = [] # files can be modified (path from self.card_dir)
     force = False
     
-    def __init__(self, running_dir, card_dir='', timeout=None, **options):
+    def __init__(self, running_dir, card_dir='', timeout=None, configuration = {},
+                 **options):
         """ initialize an object """
         
         self.running_dir = running_dir
         self.card_dir = os.path.join(self.running_dir, card_dir)
         self.timeout = timeout
-        self.found_editor()
+        self.configuration = configuration
         
         #include/overwrite options
         for key,value in options.items():
             setattr(self, key, value)
             
             
-    def found_editor(self):
-        """ found a (shell) program for editing file """
-        
-        # let first try to use the prefer editor (if EDITOR is define)
-        # if not define use the first define in a pre-define list
-            
-        if os.environ.has_key('EDITOR') and misc.which(os.environ['EDITOR']):
-            self.editor =  os.environ['EDITOR']
-            return
-        
-        logger.info('INFO: You can choose your prefered editor by defining the shell variable EDITOR')
-        
-        possibility = ['vi', 'emacs', 'vim', 'gedit', 'nano']
-        for editor in possibility:
-            if misc.which(editor):
-                self.editor = editor
-                return
-        
-        logger.warning(
-         'No valid editor detected. Please configure the shell variable EDITOR')
-            
     def run(self):
         """ execute the main code """
 
@@ -94,8 +75,8 @@
         """edit a file"""
 
         path = os.path.realpath(path)
-        subprocess.call([self.editor, path], cwd=os.getcwd())
-        
+        open_file(path, self.configuration)
+            
     def ask(self, question, default, choices=[], path_info=[]):
         """ ask a question """
         
@@ -192,16 +173,16 @@
     
     cards = ['param_card.dat', 'run_card.dat']
 
-    def __init__(self, running_dir, timeout, **option):
+    def __init__(self, running_dir, timeout, unit='pb', **option):
         """ initialize the StandAlone Version"""
         
-        
         ExtLauncher.__init__(self, running_dir, './Cards', timeout, **option)
         self.executable = os.path.join('.', 'bin','generate_events')
 
         assert hasattr(self, 'cluster')
         assert hasattr(self, 'name')
-         
+        self.unit = unit
+        
         # Check for pythia-pgs directory
         if os.path.isdir(os.path.join(MG5DIR,'pythia-pgs')):
             self.pythia = os.path.join(MG5DIR,'pythia-pgs')
@@ -221,7 +202,7 @@
         
         # Assign a valid run name if not put in options
         if self.name == '':
-            for i in range(1000):
+            for i in range(1,1000):
                 path = os.path.join(self.running_dir, 'Events','run_%02i_banner.txt' % i)
                 if not os.path.exists(path):
                     self.name = 'run_%02i' % i
@@ -282,6 +263,13 @@
     def launch_program(self):
         """launch the main program"""
         
+        # Open the corresponding crossx.html page
+        os.system('touch %s' % os.path.join(self.running_dir,'RunWeb'))
+        subprocess.call([os.path.join('bin','gen_crossxhtml-pl')], 
+                         cwd=self.running_dir)
+        open_file(os.path.join(self.running_dir, 'HTML', 'crossx.html'),
+                  self.configuration)
+
         mode = str(self.cluster)
         if mode == "0":
             subprocess.call([self.executable, mode, self.name], 
@@ -308,7 +296,10 @@
         line = fsock.readline()
         cross, error = line.split()[0:2]
         
-        logger.info('The total cross-section is %s +- %s pb' % (cross, error))
+        if self.unit != 'GeV':
+            logger.info('The total cross-section is %s +- %s %s' % (cross, error, self.unit))
+        else:
+            logger.info('The width is %s +- %s GeV' % (cross, error))
         logger.info('more information in %s' 
                                  % os.path.join(self.running_dir, 'index.html'))
                 
@@ -429,9 +420,60 @@
 
         print "Output of the run is found at " + \
               os.path.realpath(os.path.join(self.running_dir, self.name))
-        
-        
-        
-        
-                    
-            
+
+#
+# Global function to open supported file types
+#       
+def open_file(file_path, configuration):
+    """Open file if supported type, using program given by the dict configuration"""
+
+    try:
+        extension = file_path.rsplit('.',1)[1]
+    except IndexError:
+        extension = ''
+
+    # For MAC:
+    if sys.platform == 'darwin':
+        if extension == 'html':
+            html_v = configuration['web_browser']
+            if misc.which(html_v):
+                os.system('%s %s &' % (html_v, file_path))
+            elif html_v:
+                os.system('open -a %s %s' % (html_v,file_path))
+            else:
+                os.system('open %s' % file_path)
+        elif extension in ['ps', 'eps']:
+            eps_v = configuration['eps_viewer']
+            if misc.which(eps_v):
+                os.system('%s %s &' % (eps_v, file_path))
+            elif eps_v:
+                os.system('open -a %s %s' % (eps_v,file_path))
+            else:
+                os.system('open %s' % file_path)
+        else:
+            if misc.which(configuration['text_editor']):
+                subprocess.call([configuration['text_editor'], file_path])
+            else:
+                logger.warning('Unable to find a valid editor for text file. Please specify one in ./input/mg5_configuration.txt')  
+    # For Linux
+    else:
+        if extension == 'html':
+            if configuration['web_browser']:
+                os.system('%s %s &' % 
+                            (configuration['web_browser'],file_path))
+            else:
+                logger.warning('No valid web browser found. ' + \
+                               'Please set in ./input/mg5_configuration.txt') 
+        elif extension in ['ps', 'eps']:
+            if configuration['eps_viewer']:
+                os.system('%s %s &' % 
+                               (configuration['eps_viewer'],file_path))
+            else:
+                logger.warning('No valid eps viewer found. ' + \
+                               'Please set in ./input/mg5_configuration.txt') 
+        else:
+            if misc.which(configuration['text_editor']):
+                subprocess.call([configuration['text_editor'], file_path])
+            else:
+                logger.warning('No valid text editor found. ' + \
+                               'Please set in ./input/mg5_configuration.txt')

=== modified file 'madgraph/interface/tutorial_text.py'
--- madgraph/interface/tutorial_text.py	2011-03-30 17:27:49 +0000
+++ madgraph/interface/tutorial_text.py	2011-04-22 01:24:30 +0000
@@ -152,33 +152,30 @@
 or from the shell:
 ./bin/mg5 my_mg5_cmd.dat
 
-It is also possible to display this file directly from MG5 by
-launching a shell command. For example:
-mg5>shell less my_mg5_cmd.dat
+It is also possible to display this file directly from MG5:
+mg5>open ./my_mg5_cmd.dat
 """
 
-shell = """
+
+
+open = """
+Note that in order to open some file, you might be need to use a shell command.
 Any shell command can be launched by MG5, by running \"shell\" or
 starting the line by an exclamation mark (!).
 
-The final command of the tutorial is draw. This allows you to draw and
+The final command of the tutorial is display diagrams. This allows you to draw and
 look at the diagrams for your processes (in eps format) before
 creating an output for a given format. This can be useful for a fast
-check of your process.  In order to draw diagrams, you need to specify
-a directory where the eps files will be written:
-mg5>draw .
+check of your process. For this last command, we will also show how combine
+different command in a single line: 
+mg5>generate p p > go go; display diagrams
 
 Note that when you run output [madevent_v4], the diagrams are
 automatically written to the matrix.ps files in subprocess
 directory, just like with MadGraph 4.
 """
 
-draw = """
-You can look at the diagrams for example by running
-mg5>!gv ./diagrams_0_gg_ttx.eps
-or on MacOS X
-mg5>!open ./diagrams_0_gg_ttx.eps
-
+display_diagrams = """
 This command was the last step of the tutorial. 
 Quit the tutorial by typing:
 mg5>tutorial stop

=== modified file 'madgraph/iolibs/export_cpp.py'
--- madgraph/iolibs/export_cpp.py	2011-04-12 16:37:14 +0000
+++ madgraph/iolibs/export_cpp.py	2011-04-22 01:24:30 +0000
@@ -20,6 +20,7 @@
 import glob
 import itertools
 import logging
+from math import fmod
 import os
 import re
 import shutil
@@ -231,7 +232,10 @@
         else:
             self.process_string = self.processes[0].base_string()
 
-        self.process_number = process_number
+        if process_number:
+            self.process_number = process_number
+        else:
+            self.process_number = self.processes[0].get('id')
 
         self.process_name = self.get_process_name()
         self.process_class = "CPPProcess"
@@ -240,7 +244,7 @@
         self.helas_call_writer = cpp_helas_call_writer
 
         if not isinstance(self.helas_call_writer, helas_call_writers.CPPUFOHelasCallWriter):
-            raise ProcessExporterCPPError, \
+            raise self.ProcessExporterCPPError, \
                 "helas_call_writer not CPPUFOHelasCallWriter"
 
         self.nexternal, self.ninitial = \
@@ -277,6 +281,23 @@
                         wf.set('number', wf_number)
                         self.wavefunctions.append(wf)
 
+            # Also combine amplitudes
+            self.amplitudes = helas_objects.HelasAmplitudeList()
+            amp_number = 0
+            for me in self.matrix_elements:
+                for iamp, amp in enumerate(me.get_all_amplitudes()):
+                    try:
+                        old_amp = \
+                               self.amplitudes[self.amplitudes.index(amp)]
+                        amp.set('number', old_amp.get('number'))
+                    except ValueError:
+                        amp_number += 1
+                        amp.set('number', amp_number)
+                        self.amplitudes.append(amp)
+            diagram = helas_objects.HelasDiagram({'amplitudes': self.amplitudes})
+            self.amplitudes = helas_objects.HelasMatrixElement({\
+                'diagrams': helas_objects.HelasDiagramList([diagram])})
+
     # Methods for generation of process files for C++
 
     def generate_process_files(self):
@@ -411,8 +432,11 @@
                           """// Calculate wavefunctions
                           void calculate_wavefunctions(const int perm[], const int hel[]);
                           static const int nwavefuncs = %d;
-                          std::complex<double> w[nwavefuncs][18];""" % \
-                                                    len(self.wavefunctions)
+                          std::complex<double> w[nwavefuncs][18];
+                          static const int namplitudes = %d;
+                          std::complex<double> amp[namplitudes];""" % \
+                          (len(self.wavefunctions),
+                           len(self.amplitudes.get_all_amplitudes()))
             replace_dict['all_matrix_definitions'] = \
                            "\n".join(["double matrix_%s();" % \
                                       me.get('processes')[0].shell_string().\
@@ -543,7 +567,7 @@
         return ret_lines
         
 
-    def get_calculate_wavefunctions(self, wavefunctions):
+    def get_calculate_wavefunctions(self, wavefunctions, amplitudes):
         """Return the lines for optimized calculation of the
         wavefunctions for all subprocesses"""
 
@@ -555,6 +579,9 @@
             self.helas_call_writer.get_wavefunction_calls(\
             helas_objects.HelasWavefunctionList(wavefunctions)))
 
+        replace_dict['amplitude_calls'] = "\n".join(\
+            self.helas_call_writer.get_amplitude_calls(amplitudes))
+
         file = read_template_file(self.process_wavefunction_template) % \
                 replace_dict
 
@@ -645,7 +672,7 @@
                 class_name)
             ret_lines.append("// Calculate wavefunctions for all processes")
             ret_lines.append(self.get_calculate_wavefunctions(\
-                self.wavefunctions))
+                self.wavefunctions, self.amplitudes))
             ret_lines.append("}")
         else:
             ret_lines.extend([self.get_sigmaKin_single_process(i, me) \
@@ -702,6 +729,7 @@
           matrix_element.get('processes')[0].shell_string().replace("0_", "")
         
 
+        # Wavefunction and amplitude calls
         if self.single_helicities:
             replace_dict['matrix_args'] = ""
             replace_dict['all_wavefunction_calls'] = "int i, j;"
@@ -710,8 +738,9 @@
             wavefunctions = matrix_element.get_all_wavefunctions()
             replace_dict['all_wavefunction_calls'] = \
                          """const int nwavefuncs = %d;
-                         std::complex<double> w[nwavefuncs][18];\n""" % len(wavefunctions)+ \
-                         self.get_calculate_wavefunctions(wavefunctions)
+                         std::complex<double> w[nwavefuncs][18];
+                         """ % len(wavefunctions)+ \
+                         self.get_calculate_wavefunctions(wavefunctions, [])
 
         # Process name
         replace_dict['process_class_name'] = class_name
@@ -728,9 +757,6 @@
         replace_dict['color_matrix_lines'] = \
                                      self.get_color_matrix_lines(matrix_element)
 
-        # The Helicity amplitude calls
-        replace_dict['amplitude_calls'] = "\n".join(\
-            self.helas_call_writer.get_amplitude_calls(matrix_element))
         replace_dict['jamp_lines'] = self.get_jamp_lines(color_amplitudes)
 
         file = read_template_file('cpp_process_matrix.inc') % \
@@ -987,8 +1013,11 @@
                           """// Calculate wavefunctions
                           void calculate_wavefunctions(const int perm[], const int hel[]);
                           static const int nwavefuncs = %d;
-                          std::complex<double> w[nwavefuncs][18];""" % \
-                                                    len(self.wavefunctions)
+                          std::complex<double> w[nwavefuncs][18];
+                          static const int namplitudes = %d;
+                          std::complex<double> amp[namplitudes];""" % \
+                          (len(self.wavefunctions),
+                           len(self.amplitudes.get_all_amplitudes()))
             replace_dict['all_matrix_definitions'] = \
                            "\n".join(["double matrix_%s();" % \
                                       me.get('processes')[0].shell_string().\
@@ -1269,7 +1298,7 @@
 
             if final_id_list and final_mirror_id_list or \
                not final_id_list and not final_mirror_id_list:
-                raise ProcessExporterCPPError,\
+                raise self.ProcessExporterCPPError,\
                       "Missing processes, or both process and mirror process"
 
 
@@ -1289,7 +1318,7 @@
                                            (i, l) in items if l > 0]).\
                                  replace('*1', ''))
                 if any([l>1 for (i, l) in items]):
-                    raise ProcessExporterCPPError,\
+                    raise self.ProcessExporterCPPError,\
                           "More than one process with identical " + \
                           "external particles is not supported"
 
@@ -1303,7 +1332,7 @@
                                            (i, l) in items if l > 0]).\
                                  replace('*1', ''))
                 if any([l>1 for (i, l) in items]):
-                    raise ProcessExporterCPPError,\
+                    raise self.ProcessExporterCPPError,\
                           "More than one process with identical " + \
                           "external particles is not supported"
 
@@ -1386,7 +1415,7 @@
 
                 color_flows = []
                 for color_flow_dict in color_flow_list:
-                    color_flows.append([color_flow_dict[l.get('number')][i] % 500 \
+                    color_flows.append([int(fmod(color_flow_dict[l.get('number')][i], 500)) \
                                         for (l,i) in itertools.product(legs, [0,1])])
 
                 # Write out colors for the selected color flow
@@ -1855,14 +1884,21 @@
             aloha_model.compute_subset(self.wanted_lorentz)
         else:
             aloha_model.compute_all(save=False)
+            
         for abstracthelas in dict(aloha_model).values():
-            aloha_writer = aloha_writers.ALOHAWriterForCPP(abstracthelas,
-                                                           self.dir_path)
-            header = aloha_writer.define_header()
-            template_h_files.append(self.write_function_declaration(\
-                                         aloha_writer, header))
-            template_cc_files.append(self.write_function_definition(\
-                                          aloha_writer, header))
+            h_rout, cc_rout = abstracthelas.write(output_dir=None, language='CPP', 
+                                                              compiler_cmd=False)
+
+            template_h_files.append(h_rout)
+            template_cc_files.append(cc_rout)
+            
+            #aloha_writer = aloha_writers.ALOHAWriterForCPP(abstracthelas,
+            #                                               self.dir_path)
+            #header = aloha_writer.define_header()
+            #template_h_files.append(self.write_function_declaration(\
+            #                             aloha_writer, header))
+            #template_cc_files.append(self.write_function_definition(\
+            #                              aloha_writer, header))
 
         replace_dict['function_declarations'] = '\n'.join(template_h_files)
         replace_dict['function_definitions'] = '\n'.join(template_cc_files)
@@ -1902,27 +1938,27 @@
 
         return template_files
 
-    def write_function_declaration(self, aloha_writer, header):
-        """Write the function declaration for the ALOHA routine"""
-
-        ret_lines = []
-        for line in aloha_writer.write_h(header).split('\n'):
-            if self.compiler_option_re.match(line) or self.namespace_re.match(line):
-                # Strip out compiler flags and namespaces
-                continue
-            ret_lines.append(line)
-        return "\n".join(ret_lines)
-
-    def write_function_definition(self, aloha_writer, header):
-        """Write the function definition for the ALOHA routine"""
-
-        ret_lines = []
-        for line in aloha_writer.write_cc(header).split('\n'):
-            if self.compiler_option_re.match(line) or self.namespace_re.match(line):
-                # Strip out compiler flags and namespaces
-                continue
-            ret_lines.append(line)
-        return "\n".join(ret_lines)
+#    def write_function_declaration(self, aloha_writer, header):
+#        """Write the function declaration for the ALOHA routine"""
+#
+#        ret_lines = []
+#        for line in aloha_writer.write_h(header).split('\n'):
+#            if self.compiler_option_re.match(line) or self.namespace_re.match(line):
+#                # Strip out compiler flags and namespaces
+#                continue
+#            ret_lines.append(line)
+#        return "\n".join(ret_lines)
+#
+#    def write_function_definition(self, aloha_writer, header):
+#        """Write the function definition for the ALOHA routine"""
+#
+#        ret_lines = []
+#        for line in aloha_writer.write_cc(header).split('\n'):
+#            if self.compiler_option_re.match(line) or self.namespace_re.match(line):
+#                # Strip out compiler flags and namespaces
+#                continue
+#            ret_lines.append(line)
+#        return "\n".join(ret_lines)
 
     def clean_line(self, line):
         """Strip a line of compiler options and namespace options."""

=== modified file 'madgraph/iolibs/export_python.py'
--- madgraph/iolibs/export_python.py	2010-12-01 01:33:28 +0000
+++ madgraph/iolibs/export_python.py	2011-04-22 01:24:30 +0000
@@ -354,10 +354,12 @@
             parameters.remove('ZERO')
 
         # Get all couplings used
-        couplings = list(set([func.get('coupling').replace('-', '') for func \
+
+        
+        couplings = list(set([c.replace('-', '') for func \
                               in matrix_element.get_all_wavefunctions() + \
-                              matrix_element.get_all_amplitudes()
-                              if func.get('mothers')]))
+                              matrix_element.get_all_amplitudes() for c in func.get('coupling')
+                              if func.get('mothers') ]))
         
         return "\n        ".join([\
                          "%(param)s = model.get(\'parameter_dict\')[\"%(param)s\"]"\

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2011-04-13 03:40:52 +0000
+++ madgraph/iolibs/export_v4.py	2011-04-22 01:24:30 +0000
@@ -35,6 +35,7 @@
 import madgraph.iolibs.template_files as template_files
 import madgraph.iolibs.ufo_expression_parsers as parsers
 import madgraph.various.diagram_symmetry as diagram_symmetry
+import madgraph.various.process_checks as process_checks
 
 
 import aloha.create_aloha as create_aloha
@@ -51,11 +52,14 @@
     """Class to take care of exporting a set of matrix elements to
     Fortran (v4) format."""
 
-    def __init__(self, mgme_dir = "", dir_path = "", clean = False):
+    def __init__(self, mgme_dir = "", dir_path = "", clean = False,
+                 symmetry_max_time = 600):
         """Initiate the ProcessExporterFortran with directory information"""
         self.mgme_dir = mgme_dir
         self.dir_path = dir_path
         self.clean = clean
+        self.model = None
+        self.symmetry_max_time = symmetry_max_time
 
     #===========================================================================
     # copy the Template in a new directory.
@@ -178,7 +182,7 @@
         """Make the copy/symbolic links"""
         model_path = self.dir_path + '/Source/MODEL/'
         if os.path.exists(os.path.join(model_path, 'ident_card.dat')):
-            ln(model_path + '/ident_card.dat', self.dir_path + '/Cards', log=False)
+            mv(model_path + '/ident_card.dat', self.dir_path + '/Cards')
         if os.path.exists(os.path.join(model_path, 'particles.dat')):
             ln(model_path + '/particles.dat', self.dir_path + '/SubProcesses')
             ln(model_path + '/interactions.dat', self.dir_path + '/SubProcesses')
@@ -187,7 +191,9 @@
         ln(model_path + '/coupl.inc', self.dir_path + '/Source')
         ln(model_path + '/coupl.inc', self.dir_path + '/SubProcesses')
         ln(self.dir_path + '/Source/run.inc', self.dir_path + '/SubProcesses', log=False)
+        ln(self.dir_path + '/Source/genps.inc', self.dir_path + '/SubProcesses', log=False)
         ln(self.dir_path + '/Source/maxconfigs.inc', self.dir_path + '/SubProcesses', log=False)
+        ln(self.dir_path + '/Source/maxparticles.inc', self.dir_path + '/SubProcesses', log=False)
 
     #===========================================================================
     # export the helas routine
@@ -258,7 +264,7 @@
         """Write the pmass.inc file for MG4"""
 
         model = matrix_element.get('processes')[0].get('model')
-
+        
         lines = []
         for wf in matrix_element.get_external_wavefunctions():
             mass = model.get('particle_dict')[wf.get('pdg_code')].get('mass')
@@ -728,7 +734,7 @@
     # generate_subprocess_directory_v4
     #===========================================================================
     def generate_subprocess_directory_v4(self, matrix_element,
-                                                    fortran_model):
+                                         fortran_model):
         """Generate the Pxxxxx directory for a subprocess in MG4 standalone,
         including the necessary matrix.f and nexternal.inc files"""
 
@@ -785,30 +791,6 @@
                      matrix_element.get('processes')[0].nice_string())
         plot.draw()
 
-        # Generate diagrams
-        filename = "matrix.ps"
-        plot = draw.MultiEpsDiagramDrawer(matrix_element.get('base_amplitude').\
-                                             get('diagrams'),
-                                          filename,
-                                          model=matrix_element.get('processes')[0].\
-                                             get('model'),
-                                          amplitude='')
-        logger.info("Generating Feynman diagrams for " + \
-                     matrix_element.get('processes')[0].nice_string())
-        plot.draw()
-
-        # Generate diagrams
-        filename = "matrix.ps"
-        plot = draw.MultiEpsDiagramDrawer(matrix_element.get('base_amplitude').\
-                                             get('diagrams'),
-                                          filename,
-                                          model=matrix_element.get('processes')[0].\
-                                             get('model'),
-                                          amplitude='')
-        logger.info("Generating Feynman diagrams for " + \
-                     matrix_element.get('processes')[0].nice_string())
-        plot.draw()
-
         linkfiles = ['check_sa.f', 'coupl.inc', 'makefile']
 
 
@@ -935,6 +917,9 @@
         cwd = os.getcwd()
         path = os.path.join(self.dir_path, 'SubProcesses')
 
+        if not self.model:
+            self.model = matrix_element.get('processes')[0].get('model')
+
         try:
              os.chdir(path)
         except OSError, error:
@@ -1064,6 +1049,7 @@
                      'reweight.f',
                      'run.inc',
                      'maxconfigs.inc',
+                     'maxparticles.inc',
                      'setcuts.f',
                      'setscales.f',
                      'sudakov.inc',
@@ -1108,6 +1094,11 @@
         self.write_maxconfigs_file(writers.FortranWriter(filename),
                                    matrix_elements)
         
+        # Write maxparticles.inc based on max of ME's/subprocess groups
+        filename = os.path.join(self.dir_path,'Source','maxparticles.inc')
+        self.write_maxparticles_file(writers.FortranWriter(filename),
+                                     matrix_elements)
+        
         # Touch "done" file
         os.system('touch %s/done' % os.path.join(self.dir_path,'SubProcesses'))
 
@@ -1388,6 +1379,27 @@
         return True
 
     #===========================================================================
+    # write_maxparticles_file
+    #===========================================================================
+    def write_maxparticles_file(self, writer, matrix_elements):
+        """Write the maxparticles.inc file for MadEvent"""
+
+        if isinstance(matrix_elements, helas_objects.HelasMultiProcess):
+            maxparticles = max([me.get_nexternal_ninitial()[0] for me in \
+                              matrix_elements.get('matrix_elements')])
+        else:
+            maxparticles = max([me.get_nexternal_ninitial()[0] \
+                              for me in matrix_elements])
+
+        lines = "integer max_particles\n"
+        lines += "parameter(max_particles=%d)" % maxparticles
+
+        # Write the file
+        writer.writelines(lines)
+
+        return True
+
+    #===========================================================================
     # write_configs_file
     #===========================================================================
     def write_configs_file(self, writer, matrix_element):
@@ -1843,16 +1855,8 @@
         # Update values in run_config.inc
         run_config = \
                 open(os.path.join(self.dir_path, 'Source', 'run_config.inc')).read()
-        #run_config = run_config.replace("min_events_channel = 1000",
-        #                                "min_events_channel = 4000")
-        #run_config = run_config.replace("min_events = 2000",
-        #                                "min_events = 4000")
-        run_config = run_config.replace("max_events = 2000",
-                                        "max_events = 8000")
         run_config = run_config.replace("ChanPerJob=5",
                                         "ChanPerJob=2")
-        run_config = run_config.replace("nhel_survey=1",
-                                        "nhel_survey=0")
         open(os.path.join(self.dir_path, 'Source', 'run_config.inc'), 'w').\
                                     write(run_config)
         # Update values in generate_events
@@ -1878,6 +1882,11 @@
             raise base_objects.PhysicsObject.PhysicsObjectError,\
                   "subproc_group object not SubProcessGroup"
 
+        if not self.model:
+            self.model = subproc_group.get('matrix_elements')[0].\
+                         get('processes')[0].get('model')
+            self.evaluator = process_checks.MatrixElementEvaluator(self.model)
+
         cwd = os.getcwd()
         path = os.path.join(self.dir_path, 'SubProcesses')
 
@@ -2033,7 +2042,9 @@
 
         # Find config symmetries and permutations
         symmetry, perms, ident_perms = \
-                  diagram_symmetry.find_symmetry(subproc_group)
+                  diagram_symmetry.find_symmetry(subproc_group,
+                                                 self.evaluator,
+                                                 self.symmetry_max_time)
 
         filename = 'symswap.inc'
         self.write_symswap_file(writers.FortranWriter(filename),
@@ -2066,6 +2077,7 @@
                      'reweight.f',
                      'run.inc',
                      'maxconfigs.inc',
+                     'maxparticles.inc',
                      'setcuts.f',
                      'setscales.f',
                      'sudakov.inc',
@@ -2441,11 +2453,6 @@
                 double complex gal(2)
                 common/weak/ gal
 
-                double precision DUM0
-                common/FRDUM0/ DUM0
-
-                double precision DUM1
-                common/FRDUM1/ DUM1
                 """        
         fsock.writelines(header)
         
@@ -2551,10 +2558,6 @@
                              gal(2) = 1d0
                          """)
 
-        fsock.write_comments("\nDefinition of DUM symbols\n")
-        fsock.writelines(""" DUM0 = 0
-                             DUM1 = 1
-                         """)
     
     def create_couplings(self):
         """ create couplings.f and all couplingsX.f """

=== modified file 'madgraph/iolibs/files.py'
--- madgraph/iolibs/files.py	2011-03-14 16:35:14 +0000
+++ madgraph/iolibs/files.py	2011-04-22 01:24:30 +0000
@@ -165,7 +165,11 @@
     starting_dir = format_path(starting_dir)
     if not name:
         name = os.path.split(file_pos)[1]    
-        
+
+    # Remove existing link if necessary
+    if os.path.islink(os.path.join(starting_dir, name)):
+        os.remove(os.path.join(starting_dir, name))
+
     try:
         os.symlink(os.path.relpath(file_pos, starting_dir), \
                         os.path.join(starting_dir, name))

=== modified file 'madgraph/iolibs/gen_infohtml.py'
--- madgraph/iolibs/gen_infohtml.py	2011-04-11 22:45:23 +0000
+++ madgraph/iolibs/gen_infohtml.py	2011-04-22 01:24:30 +0000
@@ -104,7 +104,7 @@
     </TR>
     <TR> 
         <TD ALIGN=CENTER> <A HREF="../Source/MODEL/interactions.dat">interactions</A></TD> 
-        <TD> Interactions file used for code generation.</TABLE><P><H3 align=center><A HREF="../index.html"> Main Page </A></TD>
+        <TD> Interactions file used for code generation.</TD>
     </TR>"""
         else:
             return ''

=== modified file 'madgraph/iolibs/group_subprocs.py'
--- madgraph/iolibs/group_subprocs.py	2011-04-11 03:10:18 +0000
+++ madgraph/iolibs/group_subprocs.py	2011-04-22 01:24:30 +0000
@@ -184,6 +184,12 @@
             if part.get('mass').lower() == 'zero' and part.get('color') != 1 \
                    and part.get('spin') == 2:
                 name += "q" # "j"
+            elif part.get('mass').lower() == 'zero' and part.get('color') == 1 \
+                   and part.get('spin') == 2:
+                if part.get('charge') == 0:
+                    name += "vl"
+                else:
+                    name += "l"
             else:
                 name += part.get_name().replace('~', 'x').\
                             replace('+', 'p').replace('-', 'm')

=== modified file 'madgraph/iolibs/helas_call_writers.py'
--- madgraph/iolibs/helas_call_writers.py	2011-03-31 18:03:05 +0000
+++ madgraph/iolibs/helas_call_writers.py	2011-04-22 01:24:30 +0000
@@ -16,6 +16,7 @@
 
 import madgraph.core.base_objects as base_objects
 import madgraph.core.helas_objects as helas_objects
+import aloha.aloha_writers as aloha_writers
 from madgraph import MadGraph5Error
 
 
@@ -224,210 +225,210 @@
 
         # Gluon 4-vertex division tensor calls ggT for the FR sm and mssm
 
-        key = ((3, 3, 5, 3), 'A')
+        key = ((3, 3, 5, 3), ('A',))
         call = lambda wf: \
                "CALL UVVAXX(W(1,%d),W(1,%d),%s,zero,zero,zero,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[0].get('number'),
                 FortranHelasCallWriter.sorted_mothers(wf)[1].get('number'),
-                wf.get('coupling'),
+                wf.get('coupling')[0],
                 wf.get('number'))
         self.add_wavefunction(key, call)
 
-        key = ((3, 5, 3, 1), 'A')
+        key = ((3, 5, 3, 1), ('A',))
         call = lambda wf: \
                "CALL JVTAXX(W(1,%d),W(1,%d),%s,zero,zero,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[0].get('number'),
                 FortranHelasCallWriter.sorted_mothers(wf)[1].get('number'),
-                wf.get('coupling'),
+                wf.get('coupling')[0],
                 wf.get('number'))
         self.add_wavefunction(key, call)
 
-        key = ((3, 3, 5), 'A')
+        key = ((3, 3, 5), ('A',))
         call = lambda amp: \
                "CALL VVTAXX(W(1,%d),W(1,%d),W(1,%d),%s,zero,AMP(%d))" % \
                (FortranHelasCallWriter.sorted_mothers(amp)[0].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[1].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[2].get('number'),
-                amp.get('coupling'),
+                amp.get('coupling')[0],
                 amp.get('number'))
         self.add_amplitude(key, call)
 
         # SM gluon 4-vertex components
 
-        key = ((3, 3, 3, 3, 1), 'gggg3')
+        key = ((3, 3, 3, 3, 1), ('gggg3',))
         call = lambda wf: \
                "CALL JGGGXX(W(1,%d),W(1,%d),W(1,%d),%s,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[1].get('number'),
                 FortranHelasCallWriter.sorted_mothers(wf)[0].get('number'),
                 FortranHelasCallWriter.sorted_mothers(wf)[2].get('number'),
-                wf.get('coupling'),
+                wf.get('coupling')[0],
                 wf.get('number'))
         self.add_wavefunction(key, call)
-        key = ((3, 3, 3, 3), 'gggg1')
+        key = ((3, 3, 3, 3), ('gggg1',))
         call = lambda amp: \
                "CALL GGGGXX(W(1,%d),W(1,%d),W(1,%d),W(1,%d),%s,AMP(%d))" % \
                (FortranHelasCallWriter.sorted_mothers(amp)[0].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[1].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[2].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[3].get('number'),
-                amp.get('coupling'),
+                amp.get('coupling')[0],
                 amp.get('number'))
         self.add_amplitude(key, call)
-        key = ((3, 3, 3, 3, 1), 'gggg2')
+        key = ((3, 3, 3, 3, 1), ('gggg2',))
         call = lambda wf: \
                "CALL JGGGXX(W(1,%d),W(1,%d),W(1,%d),%s,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[0].get('number'),
                 FortranHelasCallWriter.sorted_mothers(wf)[2].get('number'),
                 FortranHelasCallWriter.sorted_mothers(wf)[1].get('number'),
-                wf.get('coupling'),
+                wf.get('coupling')[0],
                 wf.get('number'))
         self.add_wavefunction(key, call)
-        key = ((3, 3, 3, 3), 'gggg2')
+        key = ((3, 3, 3, 3), ('gggg2',))
         call = lambda amp: \
                "CALL GGGGXX(W(1,%d),W(1,%d),W(1,%d),W(1,%d),%s,AMP(%d))" % \
                (FortranHelasCallWriter.sorted_mothers(amp)[2].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[0].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[1].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[3].get('number'),
-                amp.get('coupling'),
+                amp.get('coupling')[0],
                 amp.get('number'))
         self.add_amplitude(key, call)
-        key = ((3, 3, 3, 3, 1), 'gggg1')
+        key = ((3, 3, 3, 3, 1), ('gggg1',))
         call = lambda wf: \
                "CALL JGGGXX(W(1,%d),W(1,%d),W(1,%d),%s,W(1,%d))" % \
                (FortranHelasCallWriter.sorted_mothers(wf)[2].get('number'),
                 FortranHelasCallWriter.sorted_mothers(wf)[1].get('number'),
                 FortranHelasCallWriter.sorted_mothers(wf)[0].get('number'),
-                wf.get('coupling'),
+                wf.get('coupling')[0],
                 wf.get('number'))
         self.add_wavefunction(key, call)
-        key = ((3, 3, 3, 3), 'gggg3')
+        key = ((3, 3, 3, 3), ('gggg3',))
         call = lambda amp: \
                "CALL GGGGXX(W(1,%d),W(1,%d),W(1,%d),W(1,%d),%s,AMP(%d))" % \
                (FortranHelasCallWriter.sorted_mothers(amp)[1].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[2].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[0].get('number'),
                 FortranHelasCallWriter.sorted_mothers(amp)[3].get('number'),
-                amp.get('coupling'),
-                amp.get('number'))
-        self.add_amplitude(key, call)
-
-        # HEFT VVVS calls
-
-        key = ((1, 3, 3, 3, 3), '')
-        call = lambda wf: \
-               "CALL JVVSXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
-               (wf.get('mothers')[0].get('number'),
-                wf.get('mothers')[1].get('number'),
-                wf.get('mothers')[2].get('number'),
-                wf.get('coupling'),
-                wf.get('mass'),
-                wf.get('width'),
-                wf.get('number'))
-        self.add_wavefunction(key, call)
-
-        key = ((3, 3, 3, 1, 4), '')
-        call = lambda wf: \
-               "CALL HVVVXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
-               (wf.get('mothers')[0].get('number'),
-                wf.get('mothers')[1].get('number'),
-                wf.get('mothers')[2].get('number'),
-                wf.get('coupling'),
-                wf.get('mass'),
-                wf.get('width'),
-                wf.get('number'))
-        self.add_wavefunction(key, call)
-
-        key = ((1, 3, 3, 3), '')
-        call = lambda amp: \
-               "CALL VVVSXX(W(1,%d),W(1,%d),W(1,%d),W(1,%d),DUM1,%s,AMP(%d))" % \
-               (amp.get('mothers')[0].get('number'),
-                amp.get('mothers')[1].get('number'),
-                amp.get('mothers')[2].get('number'),
-                amp.get('mothers')[3].get('number'),
-                amp.get('coupling'),
-                amp.get('number'))
-        self.add_amplitude(key, call)
-
-        # HEFT VVVS calls
-
-        key = ((1, 3, 3, 3, 1), '')
-        call = lambda wf: \
-               "CALL JVVSXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
-               (wf.get('mothers')[0].get('number'),
-                wf.get('mothers')[1].get('number'),
-                wf.get('mothers')[2].get('number'),
-                wf.get('coupling'),
-                wf.get('mass'),
-                wf.get('width'),
-                wf.get('number'))
-        self.add_wavefunction(key, call)
-
-        key = ((3, 3, 3, 1, 4), '')
-        call = lambda wf: \
-               "CALL HVVVXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
-               (wf.get('mothers')[0].get('number'),
-                wf.get('mothers')[1].get('number'),
-                wf.get('mothers')[2].get('number'),
-                wf.get('coupling'),
-                wf.get('mass'),
-                wf.get('width'),
-                wf.get('number'))
-        self.add_wavefunction(key, call)
-
-        key = ((1, 3, 3, 3), '')
-        call = lambda amp: \
-               "CALL VVVSXX(W(1,%d),W(1,%d),W(1,%d),W(1,%d),DUM1,%s,AMP(%d))" % \
-               (amp.get('mothers')[0].get('number'),
-                amp.get('mothers')[1].get('number'),
-                amp.get('mothers')[2].get('number'),
-                amp.get('mothers')[3].get('number'),
-                amp.get('coupling'),
+                amp.get('coupling')[0],
+                amp.get('number'))
+        self.add_amplitude(key, call)
+
+        # HEFT VVVS calls
+
+        key = ((1, 3, 3, 3, 3), ('',))
+        call = lambda wf: \
+               "CALL JVVSXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
+               (wf.get('mothers')[0].get('number'),
+                wf.get('mothers')[1].get('number'),
+                wf.get('mothers')[2].get('number'),
+                wf.get('coupling')[0],
+                wf.get('mass'),
+                wf.get('width'),
+                wf.get('number'))
+        self.add_wavefunction(key, call)
+
+        key = ((3, 3, 3, 1, 4), ('',))
+        call = lambda wf: \
+               "CALL HVVVXX(W(1,%d),W(1,%d),W(1,%d),DUM1,%s,%s,%s,W(1,%d))" % \
+               (wf.get('mothers')[0].get('number'),
+                wf.get('mothers')[1].get('number'),
+                wf.get('mothers')[2].get('number'),
+                wf.get('coupling')[0],
+                wf.get('mass'),
+                wf.get('width'),
+                wf.get('number'))
+        self.add_wavefunction(key, call)
+
+        key = ((1, 3, 3, 3), ('',))
+        call = lambda amp: \
+               "CALL VVVSXX(W(1,%d),W(1,%d),W(1,%d),W(1,%d),DUM1,%s,AMP(%d))" % \
+               (amp.get('mothers')[0].get('number'),
+                amp.get('mothers')[1].get('number'),