Merge lp:~maddevelopers/mg5amcnlo/reduce_aloha_call into lp:~maddevelopers/mg5amcnlo/version_1_0_1

[Olivier Mattelaer]

- Merge multiple aloha call in a single one. (for same order/color)
- Add open command in cmd_interface
- Enlarge the symmetry possibility for aloha.
- small debug/improvment

Please check carefully the CPP output. (read at least one formated output)

[Johan Alwall (johan-alwall)]

Hello Olivier,

This is really cool! I do have a couple of questions about the code:
1) Did you add the check in the ALOHA routines if COUP == 0 (in that case exit with zero wf/amp)?
   This check should not be in the combined routines of course, only in the basic routines.

2) In aloha_writers, you do
        # Define which part of the routine should be called
        if 'C' in self.namestring:
            addon = 'C' +self.namestring.split('C',1)[1]
        else:
            addon = '_%s' % self.offshell
   What if a used calls the Lorentz something with a 'C' in the name?
   (this might apply in other places too)
Update: " alohafile_pattern = re.compile(r'''^[STFV\d_]*[_C\d]*_\d%s''' % file_ext)"
   So in fact you are not allowing the user freedom to select lorentz names outside of the
   FR defined ones? This of course solves the problem.

3) + # Two possible scheme FFV1C1_2_X or FFV1__FFV2C1_X
   Why use double "_" for the combined name instead of the slightly better FFV1_FFV2C1_X?

That's the first pass, I'll continue after lunch :-)

Johan

review: Needs Information

[Olivier Mattelaer (olivier-mattelaer)]

>Hello Olivier,
>
>This is really cool! I do have a couple of questions about the code:
>1) Did you add the check in the ALOHA routines if COUP == 0 (in that case exit with zero wf/amp)?
> This check should not be in the combined routines of course, only in the basic routines.

Indeed, this is not done yet.
In fact, I'm still quite puzzled about this point.
A much better way would be to have a python script which evaluates the couplings
and simplify the matrix.f when a COUP is equal to zero(by commenting line and defining amp to zero) --off course in a reversal way--

They are many advantages to do like this:
1) you can bypass the computation of the full amplitude. (and not only on the interactions with zero couplings
2) This is done once and for all. (compare to a point by point if statement)

The disadvantages are
1) request python on the node. (We are moving for more python so this should be fine)
2) request a new compilation for each param_card (since matrix.f is modified)
3) This is much more work...

I think the advantages (speed) are enough relevant in order to pass to this method.
I would just propose to have a tag in matrix.f in order to forbid such simplification.
This will help for some test.

Note that with such method, I would need a check for the coupling be equal to zero but only in the combined routine... so exactly the opposite

>2) In aloha_writers, you do
> # Define which part of the routine should be called
> if 'C' in self.namestring:
> addon = 'C' +self.namestring.split('C',1)[1]
> else:
> addon = '_%s' % self.offshell
> What if a used calls the Lorentz something with a 'C' in the name?
> (this might apply in other places too)
>Update: " alohafile_pattern = re.compile(r'''^[STFV\d_]*[_C\d]*_\d%s''' % file_ext)"
> So in fact you are not allowing the user freedom to select lorentz names outside of the
> FR defined ones? This of course solves the problem.

Indeed not FR pattern are allowed, and indeed, name with "C" might creates trouble.
I will think about that problem, but I don't think I can find a way working in 100% of the case.
(But I can reduce the problem to routines ending by C\d)

>3) + # Two possible scheme FFV1C1_2_X or FFV1__FFV2C1_X
> Why use double "_" for the combined name instead of the slightly better FFV1_FFV2C1_X?

Ok I can change it.

[Olivier Mattelaer (olivier-mattelaer)]

Hi valentin,all,

I propose that we discuss this not on launchpad and on the ticket but
rather directly by email.
I've already send my comment on this point.

Cheers,

Olivier

On 21-avr.-11, at 11:27, Valentin Hirschi wrote:

>
>
> On Thu, Apr 21, 2011 at 11:20 AM, Valentin Hirschi <<email address hidden>
> > wrote:
>
> On Thu, Apr 21, 2011 at 10:50 AM, Olivier Mattelaer <<email address hidden>
> > wrote:
> >Hello Olivier,
> >
> >This is really cool! I do have a couple of questions about the code:
> >1) Did you add the check in the ALOHA routines if COUP == 0 (in
> that case exit with zero wf/amp)?
> > This check should not be in the combined routines of course,
> only in the basic routines.
>
> Indeed, this is not done yet.
> In fact, I'm still quite puzzled about this point.
> A much better way would be to have a python script which evaluates
> the couplings
> and simplify the matrix.f when a COUP is equal to zero(by commenting
> line and defining amp to zero) --off course in a reversal way--
>
> They are many advantages to do like this:
> 1) you can bypass the computation of the full amplitude. (and not
> only on the interactions with zero couplings
> 2) This is done once and for all. (compare to a point by point if
> statement)
>
> Wait... I did not read this correctly the first time! Do you
> actually consider some if statement (with a condition being directly
> a true/false variable) an actual computing time cost to consider???
> (for the amplitude it's the same, you can have one general if
> statement englobing the whole thing).
> I'm no expert, but I believe that this cannot possibly slow the
> code! I hardly see this if statement taking more than 1/1000
> thousand of the time used to compute a HELAS amplitude call (which
> would be the only one here at the very least).
> For NLO this is for sure totally irrelevant, and I think it's still
> negligible for LO, but if this is indeed something to consider, then
> all what I said is void and this script is indeed the way to go.
>
> The disadvantages are
> 1) request python on the node. (We are moving for more python so
> this should be fine)
> 2) request a new compilation for each param_card (since matrix.f is
> modified)
> 3) This is much more work...
>
> I think the advantages (speed) are enough relevant in order to pass
> to this method.
> I would just propose to have a tag in matrix.f in order to forbid
> such simplification.
> This will help for some test.
>
> Note that with such method, I would need a check for the coupling be
> equal to zero but only in the combined routine... so exactly the
> opposite
>
> Hi all!
>
> I'm of course not very into what you guys are doing, so all what
> follows might be completely out of context and irrelevant, so sorry
> in advance for that.
>
> I just wanted to understand why you really insist on having matrix.f
> with hardcoded elimination of amplitudes being zero because of the
> coupling constants?
> Wouldn't be nice to have a subroutine in matrix.f like the one I
> have in MadLoop which is something like InitMatrixElement()?
>
> This subroutine reads the parameters intended ...

[Johan Alwall (johan-alwall)]

Hello Olivier,

Awesome job!

Some comments on the new version:
1) when I write open <TAB> (without any generation) it only gives me
MG5_debug, so gives no hint that ./ is possible. I would prefer
MG5_debug and ./ as completions. Fixed.
2) mg5>open ./PROC_sm_1/index.html
"Unable to find application named 'None'" - I fixed this too .

For the open and launch commands:
3) We should definitely not have the same code twice, so I put all the
open code in launch. Feel free to change this.
4) On Mac, you can't use the "open" command for text files, since
there's no way to make it wait for the return. So I only allow inline
editors (emacs is always defined). So I changed this.

3) In the code:
+ short_name, addon = name.split('C',1)
+ try:
+ addon = 'C' + str(int(addon))
+ except:
+ addon = ''
+ else:
+ name = short_name
I don't understand how this solves or improves the issue.
All we need is to ensure that the user doesn't use non-FR characters,
which is what they would assume anyway. Also not nice that this comes
multiple times - it's sufficiently complicated to warrant its own
routine (which also allows to modify it easily). But I would prefer to
go back to the previous case.

But I guess it doesn't matter at this point - just to keep in mind for
the future.

Cheers,
Johan

[Johan Alwall (johan-alwall)]

Hello Olivier,

Please don't do any more edits in this branch, but work directly in the version_1_0_1 branch from now on.

Cheers,
Johan

review: Approve

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2011-04-01 22:38:18 +0000
+++ aloha/aloha_writers.py	2011-04-20 16:57:29 +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,167 @@
     
     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
+        if 'C' in self.namestring:
+            addon = 'C' +self.namestring.split('C',1)[1]
+        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:
+        name, addon = name.split('C',1)
+        addon = 'C' + addon
+    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 +612,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 +641,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 +751,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 +764,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 +774,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 +784,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 +811,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 +864,149 @@
             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
+        if 'C' in self.namestring:
+            addon = 'C' +self.namestring.split('C',1)[1]
+        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 +1016,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 +1045,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 +1059,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 +1067,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 +1089,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 +1159,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 +1187,77 @@
             
         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
+        if 'C' in self.namestring:
+            addon = 'C' +self.namestring.split('C',1)[1]
+        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-20 16:57:29 +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-03-31 18:03:05 +0000
+++ aloha/create_aloha.py	2011-04-20 16:57:29 +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)
@@ -503,7 +534,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
@@ -520,7 +558,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
@@ -589,9 +627,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):
@@ -601,7 +639,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"""
         
@@ -690,7 +767,7 @@
     aloha_files = []
     
     # Identify the valid files
-    alohafile_pattern = re.compile(r'''^[STFV]*[_C\d]*_\d%s''' % file_ext)
+    alohafile_pattern = re.compile(r'''^[STFV\d_]*[_C\d]*_\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):

=== 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-20 16:57:29 +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-20 13:10:25 +0000
+++ madgraph/interface/cmd_interface.py	2011-04-20 16:57:29 +0000
@@ -2443,7 +2443,6 @@
             raise self.InvalidCmd , '%s cannot be run from MG5 interface' % args[0]
         
         
-        
         ext_program.run()
         
         

=== modified file 'madgraph/iolibs/export_cpp.py'
--- madgraph/iolibs/export_cpp.py	2011-04-15 19:21:44 +0000
+++ madgraph/iolibs/export_cpp.py	2011-04-20 16:57:29 +0000
@@ -232,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"
@@ -1881,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)
@@ -1928,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-20 16:57:29 +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-18 15:25:50 +0000
+++ madgraph/iolibs/export_v4.py	2011-04-20 16:57:29 +0000
@@ -2477,11 +2477,6 @@
                 double complex gal(2)
                 common/weak/ gal
 
-                double precision DUM0
-                common/FRDUM0/ DUM0
-
-                double precision DUM1
-                common/FRDUM1/ DUM1
                 """        
         fsock.writelines(header)
         
@@ -2587,10 +2582,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/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-20 16:57:29 +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'),
+                amp.get('mothers')[2].get('number'),
+                amp.get('mothers')[3].get('number'),
+                amp.get('coupling')[0],
+                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')[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'),
+                amp.get('mothers')[2].get('number'),
+                amp.get('mothers')[3].get('number'),
+                amp.get('coupling')[0],
                 amp.get('number'))
         self.add_amplitude(key, call)
 
         # Spin2 Helas Routine
-        key = ((-2, 2, 5), '')
+        key = ((-2, 2, 5), ('',))
         call = lambda amp: \
                "CALL IOTXXX(W(1,%d),W(1,%d),W(1,%d),%s,%s,AMP(%d))" % \
                (amp.get('mothers')[0].get('number'),
                 amp.get('mothers')[1].get('number'),
                 amp.get('mothers')[2].get('number'),
-                amp.get('coupling'),
+                amp.get('coupling')[0],
                 amp.get('mothers')[0].get('mass'),
                 amp.get('number'))
         self.add_amplitude(key, call)
         
-        key = ((-2, 2, 5, 3), '')
+        key = ((-2, 2, 5, 3), ('',))
         call = lambda wf: \
                "CALL UIOXXX(W(1,%d),W(1,%d),%s,%s,%s,%s,W(1,%d))" % \
                (wf.get('mothers')[0].get('number'),
                 wf.get('mothers')[1].get('number'),
-                wf.get('coupling'),
+                wf.get('coupling')[0],
                 wf.get('mothers')[0].get('mass'),
                 wf.get('mass'),
                 wf.get('width'),
                 wf.get('number'))
         self.add_wavefunction(key, call)
         
-        key = ((3,3,3,5),'')
+        key = ((3,3,3,5),('',))
         call = lambda amp: \
                "CALL VVVTXX(W(1,%d),W(1,%d),W(1,%d),W(1,%d),1d0,%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) 
  
-        key = ((3,3,5),'')
+        key = ((3,3,5),('',))
         call = lambda amp: \
                "CALL VVTXXX(W(1,%d),W(1,%d),W(1,%d),%s,%s,AMP(%d))" % \
                (amp.get('mothers')[0].get('number'),
                 amp.get('mothers')[1].get('number'),
                 amp.get('mothers')[2].get('number'),
-                amp.get('coupling'),
+                amp.get('coupling')[0],
                 amp.get('mothers')[0].get('mass'),
                 amp.get('number'))
         self.add_amplitude(key, call)  
@@ -569,19 +570,18 @@
 
             # If Lorentz structure is given, by default add this
             # to call name
-            addition = argument.get('lorentz')
+            lor_name = argument.get('lorentz')[0]
 
             # Take care of special case: WWWW or WWVV calls
-            if len(argument.get('lorentz')) > 3 and \
-                   argument.get('lorentz')[:2] == "WW":
-                if argument.get('lorentz')[:4] == "WWWW":
+            if len(lor_name) > 3 and lor_name[:2] == "WW":
+                if lor_name[:4] == "WWWW":
                     mother_letters = "WWWW"[:len(mother_letters)]
-                if argument.get('lorentz')[:4] == "WWVV":
+                if lor_name[:4] == "WWVV":
                     mother_letters = "W3W3"[:len(mother_letters)]
-                addition = argument.get('lorentz')[4:]
+                lor_name = lor_name[4:]
 
             call = call + mother_letters
-            call = call + addition
+            call = call + lor_name
 
             # Check if we need to append a charge conjugation flag
             if argument.needs_hermitian_conjugate():
@@ -641,7 +641,7 @@
                                      get('number'),
                                      FortranHelasCallWriter.sorted_mothers(wf)[1].\
                                      get('number'),
-                                     wf.get_with_flow('coupling'),
+                                     ','.join(wf.get_with_flow('coupling')),
                                      wf.get('mass'),
                                      wf.get('width'),
                                      wf.get('number'))
@@ -653,7 +653,7 @@
                                      get('number'),
                                      FortranHelasCallWriter.sorted_mothers(wf)[2].\
                                      get('number'),
-                                     wf.get_with_flow('coupling'),
+                                     ','.join(wf.get_with_flow('coupling')),
                                      wf.get('mass'),
                                      wf.get('width'),
                                      wf.get('number'))
@@ -668,7 +668,7 @@
                                      FortranHelasCallWriter.sorted_mothers(amp)[2].\
                                      get('number'),
 
-                                     amp.get('coupling'),
+                                     ','.join(amp.get('coupling')),
                                      amp.get('number'))
                 else:
                     call_function = lambda amp: call % \
@@ -680,7 +680,7 @@
                                      get('number'),
                                      FortranHelasCallWriter.sorted_mothers(amp)[3].\
                                      get('number'),
-                                     amp.get('coupling'),
+                                     ','.join(amp.get('coupling')),
                                      amp.get('number'))
 
         # Add the constructed function to wavefunction or amplitude dictionary
@@ -926,12 +926,14 @@
                 outgoing = 0
 
             # Check if we need to append a charge conjugation flag
+            l = [str(l) for l in argument.get('lorentz')]
             c_flag = '' 
             if argument.needs_hermitian_conjugate():
                 c_flag = "".join(['C%d' % i for i in \
                                   argument.get_conjugate_index()])
-
-            call = 'CALL %s%s_%s' % (argument.get('lorentz'), c_flag, outgoing) 
+            routine_name = aloha_writers.combine_name(
+                                        '%s%s' % (l[0], c_flag), l[1:], outgoing)
+            call = 'CALL %s' % (routine_name)
 
             # Add the wave function
             call = call + '('
@@ -946,7 +948,7 @@
                 #CALL L_4_011(W(1,%d),W(1,%d),%s,%s, %s, W(1,%d))
                 call_function = lambda wf: call % \
                     (tuple([mother.get('number') for mother in wf.get('mothers')]) + \
-                    (wf.get_with_flow('coupling'),
+                    (','.join(wf.get_with_flow('coupling')),
                                      wf.get('mass'),
                                      wf.get('width'),
                                      wf.get('number')))
@@ -956,7 +958,7 @@
                 call_function = lambda amp: call % \
                                 (tuple([mother.get('number') 
                                           for mother in amp.get('mothers')]) + \
-                                (amp.get('coupling'),
+                                (','.join(amp.get('coupling')),
                                 amp.get('number')))     
                      
         # Add the constructed function to wavefunction or amplitude dictionary
@@ -1061,12 +1063,14 @@
                 outgoing = 0
 
             # Check if we need to append a charge conjugation flag
+            l = [str(l) for l in argument.get('lorentz')]
             c_flag = '' 
             if argument.needs_hermitian_conjugate():
                 c_flag = "".join(['C%d' % i for i in \
                                   argument.get_conjugate_index()])
-
-            call = '%s%s_%s' % (argument.get('lorentz'), c_flag, outgoing)
+            routine_name = aloha_writers.combine_name(
+                                        '%s%s' % (l[0], c_flag), l[1:], outgoing)
+            call = '%s' % (routine_name)
 
             # Add the wave function
             call = call + '('
@@ -1081,8 +1085,8 @@
                 #CALL L_4_011(W(1,%d),W(1,%d),%s,%s, %s, W(1,%d))
                 call_function = lambda wf: call % \
                     (tuple([mother.get('number')-1 for mother in wf.get('mothers')]) + \
-                    (CPPUFOHelasCallWriter.format_coupling(\
-                                     wf.get_with_flow('coupling')),
+                    (','.join(CPPUFOHelasCallWriter.format_coupling(\
+                                     wf.get_with_flow('coupling'))),
                                      wf.get('mass'),
                                      wf.get('width'),
                                      wf.get('number')-1))
@@ -1092,8 +1096,8 @@
                 call_function = lambda amp: call % \
                                 (tuple([mother.get('number')-1
                                           for mother in amp.get('mothers')]) + \
-                                (CPPUFOHelasCallWriter.format_coupling(\
-                                 amp.get('coupling')),
+                                (','.join(CPPUFOHelasCallWriter.format_coupling(\
+                                 amp.get('coupling'))),
                                  amp.get('number')-1))
                 
         # Add the constructed function to wavefunction or amplitude dictionary
@@ -1103,13 +1107,18 @@
             self.add_amplitude(argument.get_call_key(), call_function)
 
     @staticmethod
-    def format_coupling(coupling):
+    def format_coupling(couplings):
         """Format the coupling so any minus signs are put in front"""
 
-        if coupling.startswith('-'):
-            return "-pars->" + coupling[1:]
-        else:
-            return "pars->" + coupling
+        output = []
+        for coupling in couplings:
+            if coupling.startswith('-'):
+                output.append("-pars->" + coupling[1:])
+            else:
+                output.append("pars->" + coupling)
+        
+        return output
+        
 
 #===============================================================================
 # PythonUFOHelasCallWriter
@@ -1229,16 +1238,20 @@
                 outgoing = 0
 
             # Check if we need to append a charge conjugation flag
+            l = [str(l) for l in argument.get('lorentz')]
             c_flag = '' 
             if argument.needs_hermitian_conjugate():
                 c_flag = "".join(['C%d' % i for i in \
                                   argument.get_conjugate_index()])
+            routine_name = aloha_writers.combine_name(
+                                        '%s%s' % (l[0], c_flag), l[1:], outgoing)
+
 
             if isinstance(argument, helas_objects.HelasWavefunction):
                 call = 'w[%d] = '
             else:
                 call = 'amp[%d] = '
-            call += '%s%s_%s' % (argument.get('lorentz'), c_flag, outgoing) 
+            call += '%s' % routine_name
 
             # Add the wave function
             call = call + '('
@@ -1255,7 +1268,7 @@
                                 ((wf.get('number')-1,) + \
                                  tuple([mother.get('number')-1 for mother in \
                                         wf.get('mothers')]) + \
-                                 (wf.get_with_flow('coupling'),
+                                 (','.join(wf.get_with_flow('coupling')),
                                   wf.get('mass'),
                                   wf.get('width')))
             else:
@@ -1265,7 +1278,7 @@
                                 ((amp.get('number')-1,) + \
                                  tuple([mother.get('number')-1 
                                         for mother in amp.get('mothers')]) + \
-                                 (amp.get('coupling'),))
+                                 (','.join(amp.get('coupling')),))
         
         # Add the constructed function to wavefunction or amplitude dictionary
         if isinstance(argument, helas_objects.HelasWavefunction):

=== modified file 'madgraph/iolibs/template_files/matrix_madevent_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_v4.inc	2011-04-18 02:53:53 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_v4.inc	2011-04-20 16:57:29 +0000
@@ -214,6 +214,8 @@
     REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
     COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR)
     COMPLEX*16 W(18,NWAVEFUNCS)
+	COMPLEX*16 DUM0,DUM1
+    DATA DUM0, DUM1/(0d0, 0d0), (1d0, 0d0)/
 C  
 C GLOBAL VARIABLES
 C  

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_v4.inc'
--- madgraph/iolibs/template_files/matrix_standalone_v4.inc	2011-02-08 16:27:15 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_v4.inc	2011-04-20 16:57:29 +0000
@@ -92,6 +92,8 @@
       REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR)
       COMPLEX*16 W(18,NWAVEFUNCS)
+      COMPLEX*16 DUM0,DUM1
+      DATA DUM0, DUM1/(0d0, 0d0), (1d0, 0d0)/
 C  
 C GLOBAL VARIABLES
 C  

=== modified file 'tests/acceptance_tests/test_cmd.py'
--- tests/acceptance_tests/test_cmd.py	2011-04-18 19:06:57 +0000
+++ tests/acceptance_tests/test_cmd.py	2011-04-20 16:57:29 +0000
@@ -913,7 +913,7 @@
         self.do('import model sm')
         self.do('define p g u d u~ d~')
         self.do('define j g u d u~ d~')
-        self.do('generate p p > w+ j')
+        self.do('generate p p > w+ j @2')
         self.do('output pythia8 %s' % self.out_dir)
         # Check that the needed files are generated
         files = ['Processes_sm/Sigma_sm_gq_wpq.h', 'Processes_sm/Sigma_sm_gq_wpq.cc',
@@ -933,7 +933,7 @@
             shutil.rmdir(self.out_dir)
 
         self.do('import model sm')
-        self.do('generate e+ e- > e+ e-')
+        self.do('generate e+ e- > e+ e- @2')
         self.do('output standalone_cpp %s' % self.out_dir)
 
         # Check that all needed src files are generated
@@ -955,21 +955,21 @@
         subprocess.call(['make', 'check'],
                         stdout=devnull, stderr=devnull, 
                         cwd=os.path.join(self.out_dir, 'SubProcesses',
-                                         'P0_Sigma_sm_epem_epem'))
+                                         'P2_Sigma_sm_epem_epem'))
 
         self.assertTrue(os.path.exists(os.path.join(self.out_dir,
                                                     'SubProcesses',
-                                                    'P0_Sigma_sm_epem_epem',
+                                                    'P2_Sigma_sm_epem_epem',
                                                     'check')))
 
         # Check that the output of check is correct 
         logfile = os.path.join(self.out_dir, 'SubProcesses',
-                               'P0_Sigma_sm_epem_epem', 'check.log')
+                               'P2_Sigma_sm_epem_epem', 'check.log')
 
         subprocess.call('./check', 
                         stdout=open(logfile, 'w'), stderr=devnull,
                         cwd=os.path.join(self.out_dir, 'SubProcesses',
-                                         'P0_Sigma_sm_epem_epem'), shell=True)
+                                         'P2_Sigma_sm_epem_epem'), shell=True)
 
         log_output = open(logfile, 'r').read()
         me_re = re.compile('Matrix element\s*=\s*(?P<value>[\d\.e\+-]+)\s*GeV',

=== modified file 'tests/parallel_tests/me_comparator.py'
--- tests/parallel_tests/me_comparator.py	2011-04-10 15:43:51 +0000
+++ tests/parallel_tests/me_comparator.py	2011-04-20 16:57:29 +0000
@@ -339,7 +339,7 @@
         the specified model, the specified maximal coupling orders and a certain
         energy for incoming particles (for decay, incoming particle is at rest).
         """
-
+        self.res_list = [] # ensure that to be void, and avoid pointer problem 
         self.proc_list = proc_list
         self.model = model
         self.orders = orders
@@ -373,7 +373,8 @@
 
         # Get the ME value
         for i, proc in enumerate(proc_list):
-            self.res_list.append(self.get_me_value(proc, i))
+            value = self.get_me_value(proc, i)
+            self.res_list.append(value)
 
         return self.res_list
 

=== modified file 'tests/parallel_tests/sample_script.py'
--- tests/parallel_tests/sample_script.py	2011-04-12 16:11:45 +0000
+++ tests/parallel_tests/sample_script.py	2011-04-20 16:57:29 +0000
@@ -46,6 +46,7 @@
     logging.getLogger('tutorial').setLevel(logging.ERROR)
         
     logging.basicConfig(level=logging.INFO)
+<<<<<<< TREE
     #my_proc_list = ['u u~ > g y $ g', 'u~ u > g y $ g ', 'y > u u~','Z >  u u~']
     #my_proc_list = ['t t > t t', 't t~ > t t~', 't t~ > z z', 't z > t z', 't~ t~ > t~ t~', 't~ z > t~ z', 'g g > y y', 'g y > g y', 'y y > g g', 'y y > z z', 'y y > a a', 'y z > t t~', 'y z > y z', 'a y > a y',' z z > t t~', 'z z > y y', 'a a > y y']
     #my_proc_list = ['t t~ > t t~','e+ e- > e+ e-','w+ w- > w+ w-', 't t~ > t t~ / a z h g', 'y > t t~', 'y > t t~ g', 't t~ > t t~ / z h g','t t~ > t t~ / a h g','t t~ > t t~ / z a g', 't t~ > t t~ / h ', 't t~ > t t~ / z','t t~ > z > t t~', 'y t > t z', 't t~ > z > t t~']
@@ -74,6 +75,17 @@
     #print len(my_proc_list)
     #print my_proc_list
     #
+=======
+
+    my_proc_list=['z z > z z','w+ w- > w+ w-','g g > g g', 'u u~ > g g g', 'u u~ > e+ e- z', 't t~ > g g'
+                  'e+ e- > e+ e- e+ e-', ' g g > t t~ g', 't t~ > h > t t~']
+    #my_proc_list1 = me_comparator.create_proc_list(['u', 'u~','t','t~','g','y','z','a'], initial=1,
+    #                                              final=2)
+    #my_proc_list = me_comparator.create_proc_list_enhanced(
+    #    fermion, fermion, boson,
+    #    initial=2, final_1=2, final_2 = 1)
+
+>>>>>>> MERGE-SOURCE
     #my_proc_list += me_comparator.create_proc_list(['w+','w-','z','a','x1+','x1-','n1'], initial=2,
     #                                              final=3)
     #my_proc_list += me_comparator.create_proc_list(['g','u','u~','go','ul','ul~','ur','ur~'], initial=2,
@@ -102,12 +114,17 @@
 
     # Create and setup a comparator
     my_comp = me_comparator.MEComparator()
-    my_comp.set_me_runners(my_mg5_ufo, my_mg4)
+    my_comp.set_me_runners(my_mg5, my_mg5_ufo, my_mg5_cpp, my_mg4)
 
     # Run the actual comparison
     my_comp.run_comparison(my_proc_list,
+<<<<<<< TREE
                            model=['RS','RS'],
                            orders={'QED':4, 'QCD':4, 'QTD':4}, energy=2000)
+=======
+                           model='sm',
+                           orders={'QED':99, 'QCD':99, 'QTD':4}, energy=2000)
+>>>>>>> MERGE-SOURCE
 
     # Do some cleanup
     #my_comp.cleanup()

=== modified file 'tests/unit_tests/core/test_helas_objects.py'
--- tests/unit_tests/core/test_helas_objects.py	2010-12-19 03:57:07 +0000
+++ tests/unit_tests/core/test_helas_objects.py	2011-04-20 16:57:29 +0000
@@ -408,8 +408,8 @@
                                              g,
                                              g]),
                       'color': [color.ColorString([color.f(0, 1, 2)]),
-                                color.ColorString([color.f(0, 1, 2)]),
-                                color.ColorString([color.f(0, 1, 2)])],
+                                color.ColorString([color.f(1, 2, 0)]),
+                                color.ColorString([color.f(2, 0, 1)])],
                       'lorentz':['gggg1', 'gggg2', 'gggg3'],
                       'couplings':{(0, 0):'GG',(1, 1):'GG',(2, 2):'GG'},
                       'orders':{'QCD':2}}))
@@ -1137,8 +1137,10 @@
         matrix_element = helas_objects.HelasMatrixElement(\
             myamplitude,
             0)
-
-        self.assertEqual(matrix_element.get('diagrams'), diagrams)
+        
+        self.assertEqual(matrix_element.get('diagrams')[0]['amplitudes'], amplitude1, 
+                         '%s != %s' %(matrix_element.get('diagrams')[0]['amplitudes'], amplitude1))
+        #self.assertEqual(matrix_element.get('diagrams'), diagrams)
 
     def test_generate_helas_diagrams_ae_ae(self):
         """Testing the helas diagram generation a e- > a e-
@@ -1888,9 +1890,9 @@
 
         helas_writer = helas_call_writers.FortranUFOHelasCallWriter(mymodel)
 
-        self.assertEqual(len(me.get_all_amplitudes()), 8)
+        self.assertEqual(len(me.get_all_amplitudes()), 2)
 
-        self.assertEqual(len(me.get_all_wavefunctions()), 8)
+        self.assertEqual(len(me.get_all_wavefunctions()), 6)
 
         for i, amp in enumerate(me.get_all_amplitudes()):
             self.assertEqual(amp.get('number'), i + 1)

=== modified file 'tests/unit_tests/iolibs/test_export_cpp.py'
--- tests/unit_tests/iolibs/test_export_cpp.py	2011-04-15 19:21:44 +0000
+++ tests/unit_tests/iolibs/test_export_cpp.py	2011-04-20 16:57:29 +0000
@@ -665,15 +665,30 @@
   oxxxxx(p[perm[2]], mME[2], hel[2], +1, w[2]); 
   ixxxxx(p[perm[3]], mME[3], hel[3], -1, w[3]); 
   FFV1_3(w[0], w[1], pars->GC_10, pars->ZERO, pars->ZERO, w[4]); 
-  FFV2_3(w[0], w[1], pars->GC_35, pars->MZ, pars->WZ, w[5]); 
-  FFV5_3(w[0], w[1], pars->GC_47, pars->MZ, pars->WZ, w[6]); 
-  FFV1_3(w[0], w[2], pars->GC_10, pars->ZERO, pars->ZERO, w[7]); 
-  FFV2_3(w[0], w[2], pars->GC_35, pars->MZ, pars->WZ, w[8]); 
-  FFV5_3(w[0], w[2], pars->GC_47, pars->MZ, pars->WZ, w[9]); 
-
+  FFV2_5_3(w[0], w[1], pars->GC_35, pars->GC_47, pars->MZ, pars->WZ, w[5]); 
+  FFV1_3(w[0], w[2], pars->GC_10, pars->ZERO, pars->ZERO, w[6]); 
+  FFV2_5_3(w[0], w[2], pars->GC_35, pars->GC_47, pars->MZ, pars->WZ, w[7]); 
+
+<<<<<<< TREE
+=======
+
+}
+double Sigma_sm_qqx_qqx::matrix_uux_uux() 
+{
+  int i, j; 
+  // Local variables
+  const int ngraphs = 4; 
+  const int ncolor = 2; 
+  std::complex<double> ztemp; 
+  std::complex<double> amp[ngraphs], jamp[ncolor]; 
+  // The color matrix;
+  static const double denom[ncolor] = {1, 1}; 
+  static const double cf[ncolor][ncolor] = {{9, 3}, {3, 9}}; 
+>>>>>>> MERGE-SOURCE
   // Calculate all amplitudes
   // Amplitude(s) for diagram number 0
   FFV1_0(w[3], w[2], w[4], pars->GC_10, amp[0]); 
+<<<<<<< TREE
   FFV2_0(w[3], w[2], w[5], pars->GC_35, amp[1]); 
   FFV5_0(w[3], w[2], w[5], pars->GC_47, amp[2]); 
   FFV2_0(w[3], w[2], w[6], pars->GC_35, amp[3]); 
@@ -683,6 +698,14 @@
   FFV5_0(w[3], w[1], w[8], pars->GC_47, amp[7]); 
   FFV2_0(w[3], w[1], w[9], pars->GC_35, amp[8]); 
   FFV5_0(w[3], w[1], w[9], pars->GC_47, amp[9]); 
+=======
+  // Amplitude(s) for diagram number 2
+  FFV2_5_0(w[3], w[2], w[5], pars->GC_35, pars->GC_47, amp[1]); 
+  // Amplitude(s) for diagram number 3
+  FFV1_0(w[3], w[1], w[6], pars->GC_10, amp[2]); 
+  // Amplitude(s) for diagram number 4
+  FFV2_5_0(w[3], w[1], w[7], pars->GC_35, pars->GC_47, amp[3]); 
+>>>>>>> MERGE-SOURCE
 
 
 }
@@ -699,10 +722,8 @@
   static const double cf[ncolor][ncolor] = {{9, 3}, {3, 9}}; 
 
   // Calculate color flows
-  jamp[0] = +1./6. * amp[0] - amp[1] - amp[2] - amp[3] - amp[4] + 1./2. *
-      amp[5];
-  jamp[1] = -1./2. * amp[0] - 1./6. * amp[5] + amp[6] + amp[7] + amp[8] +
-      amp[9];
+  jamp[0] = +1./6. * amp[0] - amp[1] + 1./2. * amp[2]; 
+  jamp[1] = -1./2. * amp[0] - 1./6. * amp[2] + amp[3]; 
 
   // Sum and square the color flows to get the matrix element
   double matrix = 0; 
@@ -1321,6 +1342,107 @@
         print "cd /tmp; g++ -c -I $PATH_TO_PYTHIA8/include Sigma_sm_qqx_qqx.cc.cc"
 
 
+<<<<<<< TREE
+=======
+    def test_write_process_h_file(self):
+        """Test writing the .h Pythia file for a matrix element"""
+
+        goal_string = \
+"""//==========================================================================
+// This file has been automatically generated for Pythia 8
+// MadGraph 5 v. %(version)s, %(date)s
+// By the MadGraph Development Team
+// Please visit us at https://launchpad.net/madgraph5
+//==========================================================================
+
+#ifndef Pythia8_Sigma_sm_qqx_qqx_H
+#define Pythia8_Sigma_sm_qqx_qqx_H
+
+#include <complex> 
+
+#include "SigmaProcess.h"
+#include "Parameters_sm.h"
+
+using namespace std; 
+
+namespace Pythia8 
+{
+//==========================================================================
+// A class for calculating the matrix elements for
+// Process: u u~ > u u~
+// Process: c c~ > c c~
+//--------------------------------------------------------------------------
+
+class Sigma_sm_qqx_qqx : public Sigma2Process 
+{
+  public:
+
+    // Constructor.
+    Sigma_sm_qqx_qqx() {}
+
+    // Initialize process.
+    virtual void initProc(); 
+
+    // Calculate flavour-independent parts of cross section.
+    virtual void sigmaKin(); 
+
+    // Evaluate sigmaHat(sHat).
+    virtual double sigmaHat(); 
+
+    // Select flavour, colour and anticolour.
+    virtual void setIdColAcol(); 
+
+    // Evaluate weight for decay angles.
+    virtual double weightDecay(Event& process, int iResBeg, int iResEnd); 
+
+    // Info on the subprocess.
+    virtual string name() const {return "q q~ > q q~ (sm)";}
+
+    virtual int code() const {return 10000;}
+
+    virtual string inFlux() const {return "qqbarSame";}
+
+    virtual int resonanceA() const {return 23;}
+    // Tell Pythia that sigmaHat returns the ME^2
+    virtual bool convertM2() const {return true;}
+
+  private:
+
+    // Private functions to calculate the matrix element for all subprocesses
+    // Calculate wavefunctions
+    void calculate_wavefunctions(const int perm[], const int hel[]); 
+    static const int nwavefuncs = 8; 
+    std::complex<double> w[nwavefuncs][18]; 
+    double matrix_uux_uux(); 
+
+    // Constants for array limits
+    static const int nexternal = 4; 
+    static const int nprocesses = 1; 
+
+    // Store the matrix element value from sigmaKin
+    double matrix_element[nprocesses]; 
+
+    // Color flows, used when selecting color
+    double * jamp2[nprocesses]; 
+
+    // Pointer to the model parameters
+    Parameters_sm * pars; 
+
+}; 
+
+}  // end namespace Pythia
+
+#endif  // Pythia8_Sigma_sm_qqx_qqx_H
+""" % misc.get_pkg_info()
+
+        self.pythia8_exporter.write_process_h_file(\
+            writers.CPPWriter(self.give_pos('test.h')))
+
+        #print open(self.give_pos('test.h')).read()
+
+        self.assertFileContains('test.h', goal_string)
+
+>>>>>>> MERGE-SOURCE
 #===============================================================================
 # ExportUFOModelPythia8Test
 #===============================================================================

=== modified file 'tests/unit_tests/iolibs/test_export_python.py'
--- tests/unit_tests/iolibs/test_export_python.py	2011-03-24 15:44:58 +0000
+++ tests/unit_tests/iolibs/test_export_python.py	2011-04-20 16:57:29 +0000
@@ -299,9 +299,9 @@
         #  
         # Process parameters
         #  
-        ngraphs = 10
+        ngraphs = 4
         nexternal = 4
-        nwavefuncs = 10
+        nwavefuncs = 8
         ncolor = 2
         ZERO = 0.
         #  
@@ -330,32 +330,24 @@
         w[4] = FFV1_3(w[0],w[1],GC_10,ZERO, ZERO)
         # Amplitude(s) for diagram number 1
         amp[0] = FFV1_0(w[3],w[2],w[4],GC_10)
-        w[5] = FFV2_3(w[0],w[1],GC_35,MZ, WZ)
-        w[6] = FFV5_3(w[0],w[1],GC_47,MZ, WZ)
+        w[5] = FFV2_5_3(w[0],w[1],GC_35,GC_47,MZ, WZ)
         # Amplitude(s) for diagram number 2
-        amp[1] = FFV2_0(w[3],w[2],w[5],GC_35)
-        amp[2] = FFV5_0(w[3],w[2],w[5],GC_47)
-        amp[3] = FFV2_0(w[3],w[2],w[6],GC_35)
-        amp[4] = FFV5_0(w[3],w[2],w[6],GC_47)
-        w[7] = FFV1_3(w[0],w[2],GC_10,ZERO, ZERO)
+        amp[1] = FFV2_5_0(w[3],w[2],w[5],GC_35,GC_47)
+        w[6] = FFV1_3(w[0],w[2],GC_10,ZERO, ZERO)
         # Amplitude(s) for diagram number 3
-        amp[5] = FFV1_0(w[3],w[1],w[7],GC_10)
-        w[8] = FFV2_3(w[0],w[2],GC_35,MZ, WZ)
-        w[9] = FFV5_3(w[0],w[2],GC_47,MZ, WZ)
+        amp[2] = FFV1_0(w[3],w[1],w[6],GC_10)
+        w[7] = FFV2_5_3(w[0],w[2],GC_35,GC_47,MZ, WZ)
         # Amplitude(s) for diagram number 4
-        amp[6] = FFV2_0(w[3],w[1],w[8],GC_35)
-        amp[7] = FFV5_0(w[3],w[1],w[8],GC_47)
-        amp[8] = FFV2_0(w[3],w[1],w[9],GC_35)
-        amp[9] = FFV5_0(w[3],w[1],w[9],GC_47)
+        amp[3] = FFV2_5_0(w[3],w[1],w[7],GC_35,GC_47)
 
         jamp = [None] * ncolor
-        jamp[0] = +1./6.*amp[0]-amp[1]-amp[2]-amp[3]-amp[4]+1./2.*amp[5]
-        jamp[1] = -1./2.*amp[0]-1./6.*amp[5]+amp[6]+amp[7]+amp[8]+amp[9]
+        jamp[0] = +1./6.*amp[0]-amp[1]+1./2.*amp[2]
+        jamp[1] = -1./2.*amp[0]-1./6.*amp[2]+amp[3]
 
         self.amp2[0]+=abs(amp[0]*amp[0].conjugate())
-        self.amp2[1]+=abs(amp[1]*amp[1].conjugate())+abs(amp[2]*amp[2].conjugate())+abs(amp[3]*amp[3].conjugate())+abs(amp[4]*amp[4].conjugate())
-        self.amp2[2]+=abs(amp[5]*amp[5].conjugate())
-        self.amp2[3]+=abs(amp[6]*amp[6].conjugate())+abs(amp[7]*amp[7].conjugate())+abs(amp[8]*amp[8].conjugate())+abs(amp[9]*amp[9].conjugate())
+        self.amp2[1]+=abs(amp[1]*amp[1].conjugate())
+        self.amp2[2]+=abs(amp[2]*amp[2].conjugate())
+        self.amp2[3]+=abs(amp[3]*amp[3].conjugate())
 
         matrix = 0.
         for i in range(ncolor):
@@ -373,9 +365,8 @@
         matrix_methods = exporter.get_python_matrix_methods()["0_uux_uux"].\
                           split('\n')
 
-        for iline in range(len(goal_method)):
-            self.assertEqual(matrix_methods[iline],
-                             goal_method[iline])
+
+        self.assertEqual(matrix_methods, goal_method)
         
 
     def test_run_python_matrix_element(self):

=== modified file 'tests/unit_tests/iolibs/test_export_v4.py'
--- tests/unit_tests/iolibs/test_export_v4.py	2011-04-18 02:53:53 +0000
+++ tests/unit_tests/iolibs/test_export_v4.py	2011-04-20 16:57:29 +0000
@@ -276,6 +276,8 @@
       REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR)
       COMPLEX*16 W(18,NWAVEFUNCS)
+      COMPLEX*16 DUM0,DUM1
+      DATA DUM0, DUM1/(0D0, 0D0), (1D0, 0D0)/
 C     
 C     GLOBAL VARIABLES
 C     
@@ -599,13 +601,21 @@
         amp2_lines = \
                  exporter.get_amp2_lines(matrix_elements[0],
                                           subprocess_group.get('diagram_maps')[0])
+                 
         self.assertEqual(amp2_lines,
+<<<<<<< TREE
                          ['AMP2(1)=AMP2(1)+AMP(1)*dconjg(AMP(1))',
                           'AMP2(2)=AMP2(2)+AMP(2)*dconjg(AMP(2))',
                           'AMP2(3)=AMP2(3)+AMP(3)*dconjg(AMP(3))+AMP(4)*dconjg(AMP(4))+AMP(5)*dconjg(AMP(5))+AMP(6)*dconjg(AMP(6))',
                           'AMP2(4)=AMP2(4)+AMP(7)*dconjg(AMP(7))',
                           'AMP2(5)=AMP2(5)+AMP(8)*dconjg(AMP(8))',
                           'AMP2(6)=AMP2(6)+AMP(9)*dconjg(AMP(9))+AMP(10)*dconjg(AMP(10))+AMP(11)*dconjg(AMP(11))+AMP(12)*dconjg(AMP(12))'])
+=======
+                         ['AMP2(1)=AMP2(1)+AMP(1)*dconjg(AMP(1))+AMP(2)*dconjg(AMP(2))', 
+                          'AMP2(3)=AMP2(3)+AMP(3)*dconjg(AMP(3))', 
+                          'AMP2(4)=AMP2(4)+AMP(4)*dconjg(AMP(4))+AMP(5)*dconjg(AMP(5))',
+                          'AMP2(6)=AMP2(6)+AMP(6)*dconjg(AMP(6))'])
+>>>>>>> MERGE-SOURCE
         
         # Test configs.inc
 
@@ -1758,14 +1768,14 @@
         #print '\n'.join(amp2_lines)
 
         self.assertEqual('\n'.join(amp2_lines),
-"""AMP2(5)=AMP2(5)+AMP(9)*dconjg(AMP(9))
-AMP2(6)=AMP2(6)+AMP(10)*dconjg(AMP(10))
-AMP2(7)=AMP2(7)+AMP(11)*dconjg(AMP(11))
-AMP2(8)=AMP2(8)+AMP(12)*dconjg(AMP(12))
-AMP2(9)=AMP2(9)+AMP(13)*dconjg(AMP(13))
-AMP2(10)=AMP2(10)+AMP(14)*dconjg(AMP(14))
-AMP2(11)=AMP2(11)+AMP(15)*dconjg(AMP(15))
-AMP2(12)=AMP2(12)+AMP(16)*dconjg(AMP(16))""")
+"""AMP2(5)=AMP2(5)+AMP(5)*dconjg(AMP(5))
+AMP2(6)=AMP2(6)+AMP(6)*dconjg(AMP(6))
+AMP2(7)=AMP2(7)+AMP(7)*dconjg(AMP(7))
+AMP2(8)=AMP2(8)+AMP(8)*dconjg(AMP(8))
+AMP2(9)=AMP2(9)+AMP(9)*dconjg(AMP(9))
+AMP2(10)=AMP2(10)+AMP(10)*dconjg(AMP(10))
+AMP2(11)=AMP2(11)+AMP(11)*dconjg(AMP(11))
+AMP2(12)=AMP2(12)+AMP(12)*dconjg(AMP(12))""")
         
         # Test configs.inc
 
@@ -3913,13 +3923,13 @@
 
         # Test get_used_lorentz
         # Wavefunctions
-        goal_lorentz_list = [('', (), 1), ('', (), 2), ('', (), 3),
-                             ('', (1,), 2),('', (), 3), ('', (1,), 1),
-                             ('', (), 2), ('', (), 3),('', (1,), 1),
-                             ('', (), 1), ('', (), 3),('', (1,), 2),
-                             ('', (), 3), ('', (), 3)]
+        goal_lorentz_list = [(('',), (), 1), (('',), (), 2), (('',), (), 3),
+                             (('',), (1,), 2),(('',), (), 3), (('',), (1,), 1),
+                             (('',), (), 2), (('',), (), 3),(('',), (1,), 1),
+                             (('',), (), 1), (('',), (), 3),(('',), (1,), 2),
+                             (('',), (), 3), (('',), (), 3)]
         # Amplitudes
-        goal_lorentz_list += [('', (), 0)] * 8
+        goal_lorentz_list += [(('',), (), 0)] * 8
         self.assertEqual(matrix_element.get_used_lorentz(),
                          goal_lorentz_list)
 
@@ -4741,34 +4751,21 @@
 
         matrix_element = helas_objects.HelasMatrixElement(myamplitude, gen_color=False)
 
-        myfortranmodel = helas_call_writers.FortranHelasCallWriter(mybasemodel)
-        self.assertEqual("\n".join(myfortranmodel.\
-                                   get_matrix_element_calls(matrix_element)),
+        myfortranmodel = helas_call_writers.FortranUFOHelasCallWriter(mybasemodel)
+        self.assertEqual(myfortranmodel.get_matrix_element_calls(matrix_element),
                          """CALL VXXXXX(P(0,1),zero,NHEL(1),-1*IC(1),W(1,1))
 CALL VXXXXX(P(0,2),zero,NHEL(2),-1*IC(2),W(1,2))
 CALL VXXXXX(P(0,3),zero,NHEL(3),+1*IC(3),W(1,3))
 CALL VXXXXX(P(0,4),zero,NHEL(4),+1*IC(4),W(1,4))
-CALL JVVL1X(W(1,1),W(1,2),G1,zero,zero,W(1,5))
-CALL JVVL2X(W(1,1),W(1,2),G2,zero,zero,W(1,6))
+CALL L1__L2_1(W(1,1),W(1,2),G1,G2,zero, zero, W(1,5))
 # Amplitude(s) for diagram number 1
-CALL VVVL1X(W(1,3),W(1,4),W(1,5),G1,AMP(1))
-CALL VVVL2X(W(1,3),W(1,4),W(1,5),G2,AMP(2))
-CALL VVVL1X(W(1,3),W(1,4),W(1,6),G1,AMP(3))
-CALL VVVL2X(W(1,3),W(1,4),W(1,6),G2,AMP(4))
-CALL JVVL1X(W(1,1),W(1,3),G1,zero,zero,W(1,7))
-CALL JVVL2X(W(1,1),W(1,3),G2,zero,zero,W(1,8))
+CALL L1__L2_0(W(1,3),W(1,4),W(1,5),G1,G2,AMP(1))
+CALL L1__L2_1(W(1,1),W(1,3),G1,G2,zero, zero, W(1,6))
 # Amplitude(s) for diagram number 2
-CALL VVVL1X(W(1,2),W(1,4),W(1,7),G1,AMP(5))
-CALL VVVL2X(W(1,2),W(1,4),W(1,7),G2,AMP(6))
-CALL VVVL1X(W(1,2),W(1,4),W(1,8),G1,AMP(7))
-CALL VVVL2X(W(1,2),W(1,4),W(1,8),G2,AMP(8))
-CALL JVVL1X(W(1,1),W(1,4),G1,zero,zero,W(1,9))
-CALL JVVL2X(W(1,1),W(1,4),G2,zero,zero,W(1,10))
+CALL L1__L2_0(W(1,2),W(1,4),W(1,6),G1,G2,AMP(2))
+CALL L1__L2_1(W(1,1),W(1,4),G1,G2,zero, zero, W(1,7))
 # Amplitude(s) for diagram number 3
-CALL VVVL1X(W(1,2),W(1,3),W(1,9),G1,AMP(9))
-CALL VVVL2X(W(1,2),W(1,3),W(1,9),G2,AMP(10))
-CALL VVVL1X(W(1,2),W(1,3),W(1,10),G1,AMP(11))
-CALL VVVL2X(W(1,2),W(1,3),W(1,10),G2,AMP(12))""")
+CALL L1__L2_0(W(1,2),W(1,3),W(1,7),G1,G2,AMP(3))""".split('\n'))
 
         exporter = export_v4.ProcessExporterFortranME()
 
@@ -4776,9 +4773,9 @@
         amp2_lines = \
                  exporter.get_amp2_lines(matrix_element)
         self.assertEqual(amp2_lines,
-                         ['AMP2(1)=AMP2(1)+AMP(1)*dconjg(AMP(1))+AMP(2)*dconjg(AMP(2))+AMP(3)*dconjg(AMP(3))+AMP(4)*dconjg(AMP(4))',
-                          'AMP2(2)=AMP2(2)+AMP(5)*dconjg(AMP(5))+AMP(6)*dconjg(AMP(6))+AMP(7)*dconjg(AMP(7))+AMP(8)*dconjg(AMP(8))',
-                          'AMP2(3)=AMP2(3)+AMP(9)*dconjg(AMP(9))+AMP(10)*dconjg(AMP(10))+AMP(11)*dconjg(AMP(11))+AMP(12)*dconjg(AMP(12))'])
+                         ['AMP2(1)=AMP2(1)+AMP(1)*dconjg(AMP(1))',
+                          'AMP2(2)=AMP2(2)+AMP(2)*dconjg(AMP(2))',
+                          'AMP2(3)=AMP2(3)+AMP(3)*dconjg(AMP(3))'])
 
         # Test configs file
         writer = writers.FortranWriter(self.give_pos('test'))
@@ -4963,45 +4960,23 @@
 CALL IXXXXX(P(0,4),Mx1p,NHEL(4),-1*IC(4),W(1,4))
 CALL OXXXXX(P(0,5),Mx1p,NHEL(5),+1*IC(5),W(1,5))
 CALL VVV1_2(W(1,1),W(1,3),GC_214,MW, WW, W(1,6))
-CALL FFV2C1_2(W(1,4),W(1,2),GC_422,Mneu1, Wneu1, W(1,7))
-CALL FFV3C1_2(W(1,4),W(1,2),GC_628,Mneu1, Wneu1, W(1,8))
+CALL FFV2C1_3_2(W(1,4),W(1,2),GC_422,GC_628,Mneu1, Wneu1, W(1,7))
 # Amplitude(s) for diagram number 1
-CALL FFV2_0(W(1,7),W(1,5),W(1,6),GC_422,AMP(1))
-CALL FFV3_0(W(1,7),W(1,5),W(1,6),GC_628,AMP(2))
-CALL FFV2_0(W(1,8),W(1,5),W(1,6),GC_422,AMP(3))
-CALL FFV3_0(W(1,8),W(1,5),W(1,6),GC_628,AMP(4))
-CALL FFV2_1(W(1,5),W(1,2),GC_422,Mneu1, Wneu1, W(1,9))
-CALL FFV3_1(W(1,5),W(1,2),GC_628,Mneu1, Wneu1, W(1,10))
+CALL FFV2_3_0(W(1,7),W(1,5),W(1,6),GC_422,GC_628,AMP(1))
+CALL FFV2_3_1(W(1,5),W(1,2),GC_422,GC_628,Mneu1, Wneu1, W(1,8))
 # Amplitude(s) for diagram number 2
-CALL FFV2C1_0(W(1,4),W(1,9),W(1,6),GC_422,AMP(5))
-CALL FFV3C1_0(W(1,4),W(1,9),W(1,6),GC_628,AMP(6))
-CALL FFV2C1_0(W(1,4),W(1,10),W(1,6),GC_422,AMP(7))
-CALL FFV3C1_0(W(1,4),W(1,10),W(1,6),GC_628,AMP(8))
-CALL FFV2C1_2(W(1,4),W(1,1),GC_422,Mneu1, Wneu1, W(1,11))
-CALL FFV3C1_2(W(1,4),W(1,1),GC_628,Mneu1, Wneu1, W(1,12))
-CALL VVV1_2(W(1,2),W(1,3),GC_214,MW, WW, W(1,13))
+CALL FFV2C1_3_0(W(1,4),W(1,8),W(1,6),GC_422,GC_628,AMP(2))
+CALL FFV2C1_3_2(W(1,4),W(1,1),GC_422,GC_628,Mneu1, Wneu1, W(1,9))
+CALL VVV1_2(W(1,2),W(1,3),GC_214,MW, WW, W(1,10))
 # Amplitude(s) for diagram number 3
-CALL FFV2_0(W(1,11),W(1,5),W(1,13),GC_422,AMP(9))
-CALL FFV3_0(W(1,11),W(1,5),W(1,13),GC_628,AMP(10))
-CALL FFV2_0(W(1,12),W(1,5),W(1,13),GC_422,AMP(11))
-CALL FFV3_0(W(1,12),W(1,5),W(1,13),GC_628,AMP(12))
+CALL FFV2_3_0(W(1,9),W(1,5),W(1,10),GC_422,GC_628,AMP(3))
 # Amplitude(s) for diagram number 4
-CALL FFV5_0(W(1,11),W(1,9),W(1,3),GC_418,AMP(13))
-CALL FFV5_0(W(1,11),W(1,10),W(1,3),GC_418,AMP(14))
-CALL FFV5_0(W(1,12),W(1,9),W(1,3),GC_418,AMP(15))
-CALL FFV5_0(W(1,12),W(1,10),W(1,3),GC_418,AMP(16))
-CALL FFV2_1(W(1,5),W(1,1),GC_422,Mneu1, Wneu1, W(1,14))
-CALL FFV3_1(W(1,5),W(1,1),GC_628,Mneu1, Wneu1, W(1,15))
+CALL FFV5_0(W(1,9),W(1,8),W(1,3),GC_418,AMP(4))
+CALL FFV2_3_1(W(1,5),W(1,1),GC_422,GC_628,Mneu1, Wneu1, W(1,11))
 # Amplitude(s) for diagram number 5
-CALL FFV2C1_0(W(1,4),W(1,14),W(1,13),GC_422,AMP(17))
-CALL FFV3C1_0(W(1,4),W(1,14),W(1,13),GC_628,AMP(18))
-CALL FFV2C1_0(W(1,4),W(1,15),W(1,13),GC_422,AMP(19))
-CALL FFV3C1_0(W(1,4),W(1,15),W(1,13),GC_628,AMP(20))
+CALL FFV2C1_3_0(W(1,4),W(1,11),W(1,10),GC_422,GC_628,AMP(5))
 # Amplitude(s) for diagram number 6
-CALL FFV5_0(W(1,7),W(1,14),W(1,3),GC_418,AMP(21))
-CALL FFV5_0(W(1,8),W(1,14),W(1,3),GC_418,AMP(22))
-CALL FFV5_0(W(1,7),W(1,15),W(1,3),GC_418,AMP(23))
-CALL FFV5_0(W(1,8),W(1,15),W(1,3),GC_418,AMP(24))""".split('\n')
+CALL FFV5_0(W(1,7),W(1,11),W(1,3),GC_418,AMP(6))""".split('\n')
 
         for i in range(len(goal)):
             self.assertEqual(result[i], goal[i])
@@ -5012,12 +4987,12 @@
         amp2_lines = \
                  exporter.get_amp2_lines(matrix_element)
         self.assertEqual(amp2_lines,
-                         ['AMP2(1)=AMP2(1)+AMP(1)*dconjg(AMP(1))+AMP(2)*dconjg(AMP(2))+AMP(3)*dconjg(AMP(3))+AMP(4)*dconjg(AMP(4))',
-                          'AMP2(2)=AMP2(2)+AMP(5)*dconjg(AMP(5))+AMP(6)*dconjg(AMP(6))+AMP(7)*dconjg(AMP(7))+AMP(8)*dconjg(AMP(8))',
-                          'AMP2(3)=AMP2(3)+AMP(9)*dconjg(AMP(9))+AMP(10)*dconjg(AMP(10))+AMP(11)*dconjg(AMP(11))+AMP(12)*dconjg(AMP(12))',
-                          'AMP2(4)=AMP2(4)+AMP(13)*dconjg(AMP(13))+AMP(14)*dconjg(AMP(14))+AMP(15)*dconjg(AMP(15))+AMP(16)*dconjg(AMP(16))',
-                          'AMP2(5)=AMP2(5)+AMP(17)*dconjg(AMP(17))+AMP(18)*dconjg(AMP(18))+AMP(19)*dconjg(AMP(19))+AMP(20)*dconjg(AMP(20))',
-                          'AMP2(6)=AMP2(6)+AMP(21)*dconjg(AMP(21))+AMP(22)*dconjg(AMP(22))+AMP(23)*dconjg(AMP(23))+AMP(24)*dconjg(AMP(24))'])
+                         ['AMP2(1)=AMP2(1)+AMP(1)*dconjg(AMP(1))',
+                          'AMP2(2)=AMP2(2)+AMP(2)*dconjg(AMP(2))',
+                          'AMP2(3)=AMP2(3)+AMP(3)*dconjg(AMP(3))',
+                          'AMP2(4)=AMP2(4)+AMP(4)*dconjg(AMP(4))',
+                          'AMP2(5)=AMP2(5)+AMP(5)*dconjg(AMP(5))',
+                          'AMP2(6)=AMP2(6)+AMP(6)*dconjg(AMP(6))'])
 
     def test_four_fermion_vertex_normal_fermion_flow(self):
         """Testing process u u > t t g with fermion flow (u~t)(u~t)
@@ -5550,24 +5525,21 @@
 
         result = helas_call_writers.FortranUFOHelasCallWriter(mybasemodel).\
                                    get_matrix_element_calls(matrix_element.get('core_processes')[0])
-        self.assertEqual("\n".join(result),
+        self.assertEqual(result,
                          """CALL OXXXXX(P(0,1),zero,NHEL(1),-1*IC(1),W(1,1))
 CALL IXXXXX(P(0,2),MT,NHEL(2),+1*IC(2),W(1,2))
 CALL VXXXXX(P(0,3),Mwp,NHEL(3),+1*IC(3),W(1,3))
 # Amplitude(s) for diagram number 1
-CALL FFS3_0(W(1,2),W(1,1),W(1,3),GC_108,AMP(1))
-CALL FFS4_0(W(1,2),W(1,1),W(1,3),GC_111,AMP(2))""")
+CALL FFS3_4_0(W(1,2),W(1,1),W(1,3),GC_108,GC_111,AMP(1))""".split('\n'))
         result = helas_call_writers.FortranUFOHelasCallWriter(mybasemodel).\
                                    get_matrix_element_calls(matrix_element.get('decay_chains')[0].get('core_processes')[0])
-        self.assertEqual("\n".join(result),
+        self.assertEqual(result,
                          """CALL VXXXXX(P(0,1),Mwp,NHEL(1),-1*IC(1),W(1,1))
 CALL IXXXXX(P(0,2),zero,NHEL(2),-1*IC(2),W(1,2))
 CALL OXXXXX(P(0,3),MT,NHEL(3),+1*IC(3),W(1,3))
-CALL FFS3_3(W(1,2),W(1,3),GC_108,Mwp, Wwp, W(1,4))
-CALL FFS4_3(W(1,2),W(1,3),GC_111,Mwp, Wwp, W(1,5))
+CALL FFS3_4_3(W(1,2),W(1,3),GC_108,GC_111,Mwp, Wwp, W(1,4))
 # Amplitude(s) for diagram number 1
-#
-#""")
+#""".split('\n'))
 
         matrix_elements = matrix_element.combine_decay_chain_processes()
 
@@ -5580,13 +5552,9 @@
 CALL IXXXXX(P(0,2),MT,NHEL(2),+1*IC(2),W(1,2))
 CALL IXXXXX(P(0,3),zero,NHEL(3),-1*IC(3),W(1,3))
 CALL OXXXXX(P(0,4),MT,NHEL(4),+1*IC(4),W(1,4))
-CALL FFS3_3(W(1,3),W(1,4),GC_108,Mwp, Wwp, W(1,5))
-CALL FFS4_3(W(1,3),W(1,4),GC_111,Mwp, Wwp, W(1,6))
+CALL FFS3_4_3(W(1,3),W(1,4),GC_108,GC_111,Mwp, Wwp, W(1,5))
 # Amplitude(s) for diagram number 1
-CALL FFS3_0(W(1,2),W(1,1),W(1,5),GC_108,AMP(1))
-CALL FFS3_0(W(1,2),W(1,1),W(1,6),GC_108,AMP(2))
-CALL FFS4_0(W(1,2),W(1,1),W(1,5),GC_111,AMP(3))
-CALL FFS4_0(W(1,2),W(1,1),W(1,6),GC_111,AMP(4))""".split('\n')
+CALL FFS3_4_0(W(1,2),W(1,1),W(1,5),GC_108,GC_111,AMP(1))""".split('\n')
 
         for i in range(len(goal)):
             self.assertEqual(result[i], goal[i])
@@ -5705,24 +5673,21 @@
 
         result = helas_call_writers.FortranUFOHelasCallWriter(mybasemodel).\
                                    get_matrix_element_calls(matrix_element.get('core_processes')[0])
-        self.assertEqual("\n".join(result),
+        self.assertEqual(result,
                          """CALL OXXXXX(P(0,1),zero,NHEL(1),-1*IC(1),W(1,1))
 CALL IXXXXX(P(0,2),MT,NHEL(2),+1*IC(2),W(1,2))
 CALL SXXXXX(P(0,3),+1*IC(3),W(1,3))
 # Amplitude(s) for diagram number 1
-CALL FFS3C1_0(W(1,2),W(1,1),W(1,3),GC_108,AMP(1))
-CALL FFS4C1_0(W(1,2),W(1,1),W(1,3),GC_111,AMP(2))""")
+CALL FFS3C1_4_0(W(1,2),W(1,1),W(1,3),GC_108,GC_111,AMP(1))""".split('\n'))
         result = helas_call_writers.FortranUFOHelasCallWriter(mybasemodel).\
                                    get_matrix_element_calls(matrix_element.get('decay_chains')[0].get('core_processes')[0])
-        self.assertEqual("\n".join(result),
+        self.assertEqual(result,
                          """CALL SXXXXX(P(0,1),-1*IC(1),W(1,1))
 CALL OXXXXX(P(0,2),zero,NHEL(2),+1*IC(2),W(1,2))
 CALL IXXXXX(P(0,3),MT,NHEL(3),-1*IC(3),W(1,3))
-CALL FFS3C1_3(W(1,3),W(1,2),GC_108,Msix1, Wsix1, W(1,4))
-CALL FFS4C1_3(W(1,3),W(1,2),GC_111,Msix1, Wsix1, W(1,5))
+CALL FFS3C1_4_3(W(1,3),W(1,2),GC_108,GC_111,Msix1, Wsix1, W(1,4))
 # Amplitude(s) for diagram number 1
-#
-#""")
+#""".split('\n'))
 
         matrix_elements = matrix_element.combine_decay_chain_processes()
 
@@ -5735,16 +5700,12 @@
 CALL IXXXXX(P(0,2),MT,NHEL(2),+1*IC(2),W(1,2))
 CALL OXXXXX(P(0,3),zero,NHEL(3),+1*IC(3),W(1,3))
 CALL IXXXXX(P(0,4),MT,NHEL(4),-1*IC(4),W(1,4))
-CALL FFS3C1_3(W(1,4),W(1,3),GC_108,Msix1, Wsix1, W(1,5))
-CALL FFS4C1_3(W(1,4),W(1,3),GC_111,Msix1, Wsix1, W(1,6))
+CALL FFS3C1_4_3(W(1,4),W(1,3),GC_108,GC_111,Msix1, Wsix1, W(1,5))
 # Amplitude(s) for diagram number 1
-CALL FFS3C1_0(W(1,2),W(1,1),W(1,5),GC_108,AMP(1))
-CALL FFS3C1_0(W(1,2),W(1,1),W(1,6),GC_108,AMP(2))
-CALL FFS4C1_0(W(1,2),W(1,1),W(1,5),GC_111,AMP(3))
-CALL FFS4C1_0(W(1,2),W(1,1),W(1,6),GC_111,AMP(4))""".split('\n')
+CALL FFS3C1_4_0(W(1,2),W(1,1),W(1,5),GC_108,GC_111,AMP(1))""".split('\n')
 
-        for i in range(len(goal)):
-            self.assertEqual(result[i], goal[i])
+        
+        self.assertEqual(result, goal)
 
     def test_matrix_multistage_decay_chain_process(self):
         """Test matrix.f for multistage decay chain
@@ -6033,8 +5994,8 @@
                                              g,
                                              g]),
                       'color': [color.ColorString([color.f(0, 1, 2)]),
-                                color.ColorString([color.f(0, 1, 2)]),
-                                color.ColorString([color.f(0, 1, 2)])],
+                                color.ColorString([color.f(2, 1, 0)]),
+                                color.ColorString([color.f(1, 0, 2)])],
                       'lorentz':['gggg1', 'gggg2', 'gggg3'],
                       'couplings':{(0, 0):'GG', (1, 1):'GG', (2, 2):'GG'},
                       'orders':{'QCD':2}}))
@@ -7760,7 +7721,7 @@
                                               gen_color=True)
 
         self.assertEqual(sum([len(diagram.get('amplitudes')) for diagram in \
-                          me.get('diagrams')]), 8)
+                          me.get('diagrams')]), 2)
 
         for i, amp in enumerate(me.get_all_amplitudes()):
             self.assertEqual(amp.get('number'), i + 1)
@@ -7770,8 +7731,8 @@
         exporter = export_v4.ProcessExporterFortranME()
 
         self.assertEqual(exporter.get_JAMP_lines(me),
-                         ["JAMP(1)=-AMP(1)-AMP(2)-AMP(3)-AMP(4)",
-                         "JAMP(2)=+AMP(5)+AMP(6)+AMP(7)+AMP(8)"])
+                         ["JAMP(1)=-AMP(1)",
+                         "JAMP(2)=+AMP(2)"])
 
     def test_generate_helas_diagrams_gg_gogo(self):
         """Testing the v4 helas diagram generation g g > go go,
@@ -8631,12 +8592,12 @@
 C     The process calculated in this file is: 
 C     Gamma(3,2,1)
 C     
-      SUBROUTINE FFV1_1(F2, V3, C, M1, W1, F1)
+      SUBROUTINE FFV1_1(F2, V3, COUP, M1, W1, F1)
       IMPLICIT NONE
       DOUBLE COMPLEX F1(6)
       DOUBLE COMPLEX F2(6)
       DOUBLE COMPLEX V3(6)
-      DOUBLE COMPLEX C
+      DOUBLE COMPLEX COUP
       DOUBLE COMPLEX DENOM
       DOUBLE PRECISION M1, W1
       DOUBLE PRECISION P1(0:3)
@@ -8653,10 +8614,9 @@
         abstract_M.write('/tmp','Fortran')
         
         self.assertTrue(os.path.exists('/tmp/FFV1_1.f'))
-        textfile = open('/tmp/FFV1_1.f','r')
+        textfile = open('/tmp/FFV1_1.f','r').read()
         split_sol = solution.split('\n')
-        for i in range(len(split_sol)):
-            self.assertEqual(split_sol[i]+'\n', textfile.readline())
+        self.assertEqual(split_sol, textfile.split('\n')[:len(split_sol)])
 
 
 class UFO_model_to_mg4_Test(unittest.TestCase):

=== modified file 'tests/unit_tests/iolibs/test_helas_call_writers.py'
--- tests/unit_tests/iolibs/test_helas_call_writers.py	2011-03-22 00:16:20 +0000
+++ tests/unit_tests/iolibs/test_helas_call_writers.py	2011-04-20 16:57:29 +0000
@@ -493,14 +493,14 @@
 
         wavefunctions = {}
         # IXXXXXX.Key: (spin, state)
-        key1 = ((-2, 0), '')
+        key1 = ((-2, 0), ('',))
         wavefunctions[key1] = \
                           lambda wf: 'CALL IXXXXX(P(0,%d),%s,NHEL(%d),%d*IC(%d),W(1,%d))' % \
                           (wf.get('number_external'), wf.get('mass'),
                            wf.get('number_external'), -(-1) ** wf.get_with_flow('is_part'),
                            wf.get('number_external'), wf.get('number'))
         # OXXXXXX.Key: (spin, state)
-        key2 = ((2, 0), '')
+        key2 = ((2, 0), ('',))
         wavefunctions[key2] = \
                           lambda wf: 'CALL OXXXXX(P(0,%d),%s,NHEL(%d),%d*IC(%d),W(1,%d))' % \
                           (wf.get('number_external'), wf.get('mass'),

=== modified file 'tests/unit_tests/various/test_4fermion_models.py'
--- tests/unit_tests/various/test_4fermion_models.py	2011-04-15 04:15:48 +0000
+++ tests/unit_tests/various/test_4fermion_models.py	2011-04-20 16:57:29 +0000
@@ -77,7 +77,6 @@
                                                               p)[0]
 
             
-
         self.assertAlmostEqual(values['scalar'], values['4ferm'], 3)
 
 #===============================================================================

=== modified file 'tests/unit_tests/various/test_aloha.py'
--- tests/unit_tests/various/test_aloha.py	2011-03-31 17:39:59 +0000
+++ tests/unit_tests/various/test_aloha.py	2011-04-20 16:57:29 +0000
@@ -16,6 +16,7 @@
 the output of the Feynman Rules."""
 from __future__ import division
 
+import os
 import time
 import aloha.aloha_object as aloha_obj
 import aloha.aloha_lib as aloha_lib
@@ -2798,7 +2799,7 @@
         # Check that full identification symmetry works
         helas_suite = create_aloha.AbstractALOHAModel('sm')
         helas_suite.look_for_symmetries()
-        solution = {'VVS1': {2: 1}, 'SSS1': {2: 1, 3: 2}, 'VVSS1': {2: 1, 4: 3}, 'VVVV2': {2: 1, 4: 3}, 'SSSS1': {2: 1, 3: 2, 4: 3}} 
+        solution = {'SSSS1': {2: 1, 3: 2, 4: 3}, 'VVVV4': {2: 1, 3: 2, 4: 3}, 'VVV1': {2: 1, 3: 2}, 'VVVV1': {2: 1, 3: 2, 4: 3}, 'VVVV3': {2: 1, 3: 2, 4: 3}, 'VVVV2': {2: 1, 4: 3}, 'SSS1': {2: 1, 3: 2}, 'VVSS1': {2: 1, 4: 3}, 'VVS1': {2: 1}}  
         self.assertEqual(solution, helas_suite.symmetries)
         
     def test_has_symmetries(self):
@@ -2824,6 +2825,246 @@
         
         base = helas_suite.has_symmetries('VVS1', 3, valid_output=(1, 2))
         self.assertEqual(base, None)   
+
+    def test_aloha_multiple_lorentz(self):
+        """ check if the detection of multiple lorentz work """
+        
+        helas_suite = create_aloha.AbstractALOHAModel('sm')
+        helas_suite.look_for_multiple_lorentz_interactions()
+        solution = {'FFV2': [('FFV3',), ('FFV4',), ('FFV5',)], 'FFS3': [('FFS4',)]}
+        self.assertEqual(solution, helas_suite.multiple_lor)
+        
+
+    def test_aloha_multiple_lorentz_and_symmetry(self):
+        """ check if the detection of multiple lorentz work """
+        
+        VVS1 = self.Lorentz(name = 'VVS1',
+                 spins = [ 3, 3, 1 ],
+                 structure = 'Metric(1,2)')
+
+        #VVS2 = self.Lorentz(name = 'VVS2',
+        #         spins = [ 3, 3, 1 ],
+        #         structure = 'Metric(2,1)')
+        
+        abstract = create_aloha.AbstractRoutineBuilder(VVS1).compute_routine(1)
+        abstract.add_symmetry(2)
+        abstract.add_combine(('VVS2',))
+        
+        text =  abstract.write(None, 'Fortran')
+
+        goal =""" subroutine VVS1_1(V2, S3, COUP, M1, W1, V1)
+implicit none 
+double complex V1(6)
+double complex V2(6)
+double complex S3(3)
+double complex COUP
+double complex denom
+double precision M1, W1
+double complex OM1
+double precision P1(0:3)
+
+V1(5)= V2(5)+S3(2)
+V1(6)= V2(6)+S3(3)
+P1(0) = - dble(V1(5))
+P1(1) = - dble(V1(6))
+P1(2) = - dimag(V1(6))
+P1(3) = - dimag(V1(5))
+OM1 = 0d0
+if (M1 .ne. 0d0) OM1=1d0/M1**2
+
+denom =1d0/(( (M1*( -M1+(0, 1)*W1))+( (P1(0)**2)-(P1(1)**2)-(P1(2)**2)-(P1(3)**2))))
+V1(1)= COUP*denom*(S3(1)*( (OM1*( (0, 1)*(V2(1)*P1(0))+(0, -1)*(V2(2)*P1(1))+(0, -1)*(V2(3)*P1(2))+(0, -1)*(V2(4)*P1(3)))*P1(0))+(0, -1)*V2(1)))
+V1(2)= COUP*denom*(S3(1)*( (OM1*( (0, 1)*(V2(1)*P1(0))+(0, -1)*(V2(2)*P1(1))+(0, -1)*(V2(3)*P1(2))+(0, -1)*(V2(4)*P1(3)))*P1(1))+(0, -1)*V2(2)))
+V1(3)= COUP*denom*(S3(1)*( (OM1*( (0, 1)*(V2(1)*P1(0))+(0, -1)*(V2(2)*P1(1))+(0, -1)*(V2(3)*P1(2))+(0, -1)*(V2(4)*P1(3)))*P1(2))+(0, -1)*V2(3)))
+V1(4)= COUP*denom*(S3(1)*( (OM1*( (0, 1)*(V2(1)*P1(0))+(0, -1)*(V2(2)*P1(1))+(0, -1)*(V2(3)*P1(2))+(0, -1)*(V2(4)*P1(3)))*P1(3))+(0, -1)*V2(4)))
+end
+
+
+
+ subroutine VVS1_2(V2, S3, COUP, M1, W1, V1)
+implicit none 
+double complex V1(6)
+double complex V2(6)
+double complex S3(3)
+double complex COUP
+double complex denom
+double precision M1, W1
+double complex OM1
+double precision P1(0:3)
+call VVS1_1(V2,S3,COUP,M1,W1,V1)
+end
+
+
+ subroutine VVS1_2_1(V2, S3, COUP1,COUP2, M1, W1, V1)
+implicit none 
+double complex V1(6)
+double complex V2(6)
+double complex S3(3)
+double complex COUP1,COUP2
+double complex denom
+double precision M1, W1
+double complex OM1
+double precision P1(0:3)
+ double complex TMP(6)
+ integer i
+
+ CALL VVS1_1(V2, S3, COUP1, M1, W1, V1)
+ CALL VVS2_1(V2, S3, COUP2, M1, W1, TMP)
+ do i=1,4
+                V1(i) = V1(i) + tmp(i)
+                enddo
+end
+
+ subroutine VVS1_2_2(V2, S3, COUP1,COUP2, M1, W1, V1)
+implicit none 
+double complex V1(6)
+double complex V2(6)
+double complex S3(3)
+double complex COUP1,COUP2
+double complex denom
+double precision M1, W1
+double complex OM1
+double precision P1(0:3)
+ double complex TMP(6)
+ integer i
+
+ CALL VVS1_1(V2, S3, COUP1, M1, W1, V1)
+ CALL VVS2_1(V2, S3, COUP2, M1, W1, TMP)
+ do i=1,4
+                V1(i) = V1(i) + tmp(i)
+                enddo
+end
+
+"""
+        self.assertEqual(text.split('\n'),goal.split('\n')) 
+   
+        text_h, text_cpp =  abstract.write(None, 'CPP')
+    
+        goal_h = """#ifndef VVS1_1_guard
+#define VVS1_1_guard
+#include <complex>
+using namespace std;
+
+void VVS1_1(complex<double> V2[],complex<double> S3[],complex<double> COUP, double M1, double W1, complex<double>V1[]);
+
+void VVS1_2(complex<double> V2[],complex<double> S3[],complex<double> COUP, double M1, double W1, complex<double>V1[]);
+
+#endif
+
+#ifndef VVS1_2_1_guard
+#define VVS1_2_1_guard
+#include <complex>
+using namespace std;
+
+void VVS1_2_1(complex<double> V2[],complex<double> S3[],complex<double> COUP1, complex <double>COUP2, double M1, double W1, complex<double>V1[]);
+
+void VVS1_2_2(complex<double> V2[],complex<double> S3[],complex<double> COUP1, complex <double>COUP2, double M1, double W1, complex<double>V1[]);
+
+#endif"""
+        goal_cpp = """#include "VVS1_1.h"
+
+void VVS1_1(complex<double> V2[],complex<double> S3[],complex<double> COUP, double M1, double W1, complex<double>V1[]){
+complex<double> denom;
+complex<double> OM1;
+double P1[4];
+V1[4]= V2[4]+S3[1];
+V1[5]= V2[5]+S3[2];
+P1[0] = -V1[4].real();
+P1[1] = -V1[5].real();
+P1[2] = -V1[5].imag();
+P1[3] = -V1[4].imag();
+OM1 = 0;
+if (M1 != 0) OM1= 1./pow(M1,2);
+denom =1./(( (M1*( -M1+complex<double>(0., 1.)*W1))+( (pow(P1[0],2))-(pow(P1[1],2))-(pow(P1[2],2))-(pow(P1[3],2)))));
+V1[0]= COUP*denom*(S3[0]*( (OM1*( complex<double>(0., 1.)*(V2[0]*P1[0])+complex<double>(0., -1.)*(V2[1]*P1[1])+complex<double>(0., -1.)*(V2[2]*P1[2])+complex<double>(0., -1.)*(V2[3]*P1[3]))*P1[0])+complex<double>(0., -1.)*V2[0]));
+V1[1]= COUP*denom*(S3[0]*( (OM1*( complex<double>(0., 1.)*(V2[0]*P1[0])+complex<double>(0., -1.)*(V2[1]*P1[1])+complex<double>(0., -1.)*(V2[2]*P1[2])+complex<double>(0., -1.)*(V2[3]*P1[3]))*P1[1])+complex<double>(0., -1.)*V2[1]));
+V1[2]= COUP*denom*(S3[0]*( (OM1*( complex<double>(0., 1.)*(V2[0]*P1[0])+complex<double>(0., -1.)*(V2[1]*P1[1])+complex<double>(0., -1.)*(V2[2]*P1[2])+complex<double>(0., -1.)*(V2[3]*P1[3]))*P1[2])+complex<double>(0., -1.)*V2[2]));
+V1[3]= COUP*denom*(S3[0]*( (OM1*( complex<double>(0., 1.)*(V2[0]*P1[0])+complex<double>(0., -1.)*(V2[1]*P1[1])+complex<double>(0., -1.)*(V2[2]*P1[2])+complex<double>(0., -1.)*(V2[3]*P1[3]))*P1[3])+complex<double>(0., -1.)*V2[3]));
+}
+
+void VVS1_2(complex<double> V2[],complex<double> S3[],complex<double> COUP, double M1, double W1, complex<double>V1[]){
+VVS1_1(V2,S3,COUP,M1,W1,V1);
+}
+
+#include "VVS1_2_1.h"
+
+void VVS1_2_1(complex<double> V2[],complex<double> S3[],complex<double> COUP1, complex<double>COUP2, double M1, double W1, complex<double>V1[]){
+complex<double> tmp[6];
+ int i = 0;
+
+VVS1_1(V2, S3, COUP1, M1, W1, V1);
+VVS2_1(V2, S3, COUP2, M1, W1, tmp);
+ while (i < 4)
+                {
+                V1[i] = V1[i] + tmp[i];
+                i++;
+                }
+}
+
+#include "VVS1_2_2.h"
+
+void VVS1_2_2(complex<double> V2[],complex<double> S3[],complex<double> COUP1, complex<double>COUP2, double M1, double W1, complex<double>V1[]){
+complex<double> tmp[6];
+ int i = 0;
+
+VVS1_1(V2, S3, COUP1, M1, W1, V1);
+VVS2_1(V2, S3, COUP2, M1, W1, tmp);
+ while (i < 4)
+                {
+                V1[i] = V1[i] + tmp[i];
+                i++;
+                }
+}
+
+"""
+        
+        self.assertEqual(text_h.split('\n'),goal_h.split('\n'))
+        self.assertEqual(text_cpp.split('\n'),goal_cpp.split('\n'))
+        
+        
+        text =  abstract.write(None, 'Python')
+        goal = """import wavefunctions
+def VVS1_1(V2, S3, COUP, M1, W1):
+    V1 = wavefunctions.WaveFunction(size=6)
+    V1[4] = V2[4]+S3[1]
+    V1[5] = V2[5]+S3[2]
+    P1 = [-complex(V1[4]).real, \\
+            - complex(V1[5]).real, \\
+            - complex(V1[5]).imag, \\
+            - complex(V1[4]).imag]
+    OM1 = 0.0
+    if (M1): OM1=1.0/M1**2
+    denom =1.0/(( (M1*( -M1+1j*W1))+( (P1[0]**2)-(P1[1]**2)-(P1[2]**2)-(P1[3]**2))))
+    V1[0]= COUP*denom*(S3[0]*( (OM1*( 1j*(V2[0]*P1[0])-1j*(V2[1]*P1[1])-1j*(V2[2]*P1[2])-1j*(V2[3]*P1[3]))*P1[0])-1j*V2[0]))
+    V1[1]= COUP*denom*(S3[0]*( (OM1*( 1j*(V2[0]*P1[0])-1j*(V2[1]*P1[1])-1j*(V2[2]*P1[2])-1j*(V2[3]*P1[3]))*P1[1])-1j*V2[1]))
+    V1[2]= COUP*denom*(S3[0]*( (OM1*( 1j*(V2[0]*P1[0])-1j*(V2[1]*P1[1])-1j*(V2[2]*P1[2])-1j*(V2[3]*P1[3]))*P1[2])-1j*V2[2]))
+    V1[3]= COUP*denom*(S3[0]*( (OM1*( 1j*(V2[0]*P1[0])-1j*(V2[1]*P1[1])-1j*(V2[2]*P1[2])-1j*(V2[3]*P1[3]))*P1[3])-1j*V2[3]))
+    return V1
+    
+    
+def VVS1_2(V2, S3, COUP, M1, W1):
+    return VVS1_1(V2,S3,COUP,M1,W1)
+
+def VVS1_2_1(V2, S3, COUP1,COUP2, M1, W1):
+
+
+    V1 = VVS1_1(V2, S3, COUP1, M1, W1)
+    tmp = VVS2_1(V2, S3, COUP2, M1, W1)
+    for i in range(4):
+        V1[i] += tmp[i]
+    return V1
+
+def VVS1_2_2(V2, S3, COUP1,COUP2, M1, W1):
+
+
+    V1 = VVS1_1(V2, S3, COUP1, M1, W1)
+    tmp = VVS2_1(V2, S3, COUP2, M1, W1)
+    for i in range(4):
+        V1[i] += tmp[i]
+    return V1
+
+"""
+        self.assertEqual(text.split('\n'),goal.split('\n'))
         
         
     def test_full_sm_aloha(self):
@@ -2852,17 +3093,77 @@
             spin_solution = spin_index[helas.spins[output_part -1]]
             self.assertEqual(abstract.expr.numerator.nb_spin, spin_solution, \
                              error % name)
+            
+    def test_multiple_lorentz_subset(self):
+        """test if we create the correct set of routine/files for multiple lorentz"""
+        
+        helas_suite = create_aloha.AbstractALOHAModel('sm')
+        requested_routines=[(('FFS1','FFS2') , (), 0), 
+                            (('FFV1',) , (), 0), 
+                            (('FFV1','FFV2') , (1,), 0)]
+        
+        helas_suite.compute_subset(requested_routines)
+
+        # Check that the 5 base routines are created
+        # FFS1, FFS2, FFV1, FFV1C1, FFV2C1
+        self.assertEqual(len(helas_suite), 5)
+        
+        # Check that FFS1 and FFV1C1 are correctly connected to the associate
+        #lorentz
+        linked = helas_suite[('FFS1',0)].combined
+        self.assertEqual(linked, [('FFS2',)])
+        linked = helas_suite[('FFV1C1',0)].combined
+        self.assertEqual(linked, [('FFV2',)])        
+        linked = helas_suite[('FFV1',0)].combined
+        self.assertEqual(linked, [])
+        
+        # Check that the file are correctly written
+        os.system('rm -r /tmp/mg5 &> /dev/null; mkdir /tmp/mg5 &> /dev/null')
+        helas_suite.write('/tmp/mg5', 'Fortran')
+        
+        content = set(os.listdir('/tmp/mg5'))
+        self.assertEqual(content, set(['FFS1_0.f', 'FFS2_0.f', 'FFV1_0.f',
+                                       'FFV1C1_0.f','FFV2C1_0.f', 'FFS1_2_0.f',
+                                       'FFV1C1_2_0.f']))
+        
+        # Check the content of FFV1__FFV2C1_0.f
+        fsock = open('/tmp/mg5/FFV1C1_2_0.f')
+        goal = """C     This File is Automatically generated by ALOHA 
+C     
+      SUBROUTINE FFV1C1_2_0(F1,F2,V3,COUP1,COUP2,VERTEX)
+      IMPLICIT NONE
+      DOUBLE COMPLEX F1(6)
+      DOUBLE COMPLEX F2(6)
+      DOUBLE COMPLEX V3(6)
+      DOUBLE COMPLEX COUP1,COUP2
+      DOUBLE COMPLEX VERTEX
+      DOUBLE COMPLEX TMP
+
+
+      CALL FFV1C1_0(F1, F2, V3, COUP1, VERTEX)
+      CALL FFV2C1_0(F1, F2, V3, COUP2, TMP)
+      VERTEX = VERTEX + TMP
+      END
+
+
+"""
+        self.assertEqual(fsock.read().split('\n'), goal.split('\n'))
+        
+        
+        
+        
+        
     
     def test_mssm_subset_creation(self):
         """ test the creation of subpart of ALOHA routines 
         including clash routines """
         helas_suite = create_aloha.AbstractALOHAModel('mssm')
         
-        requested_routines=[('FFV1' , (), 0), 
-                            ('FFV1', (), 2),
-                            ('FFV1', (1,), 0),
-                            ('FFV2', (1,), 3),
-                            ('VVV1', (), 3)]
+        requested_routines=[(('FFV1',) , (), 0), 
+                            (('FFV1',), (), 2),
+                            (('FFV1',), (1,), 0),
+                            (('FFV2',), (1,), 3),
+                            (('VVV1',), (), 3)]
         
         helas_suite.compute_subset(requested_routines)        
         self.assertEqual(len(helas_suite), 5)
@@ -3008,7 +3309,7 @@
     and test_export_pythia"""
     
     
-    def test_reorder_call_listFFVV(self):
+    def old_test_reorder_call_listFFVV(self):
         
         FFVV = UFOLorentz(name = 'FFVV',
                spins = [ 2, 2, 3, 3])
@@ -3021,7 +3322,7 @@
         new_call = writer.reorder_call_list(call_list, 1, 2)
         self.assertEqual(['F2', 'V3', 'V4'], new_call)
 
-    def test_reorder_call_listFVVV(self):
+    def old_test_reorder_call_listFVVV(self):
         FVVV = UFOLorentz(name = 'FVVV',
                spins = [ 2, 3, 3, 3])
         
@@ -3041,7 +3342,7 @@
         new_call = writer.reorder_call_list(call_list, 2, 4)
         self.assertEqual(['F1', 'V4', 'V3'], new_call)                  
     
-    def test_reorder_call_listVVVV(self):
+    def old_test_reorder_call_listVVVV(self):
         VVVV = UFOLorentz(name = 'VVVV',
                spins = [ 3, 3, 3, 3])
     
@@ -3092,7 +3393,7 @@
         new_call = writer.reorder_call_list(call_list, 3, 4)
         self.assertEqual(['V1', 'V2', 'V4'], new_call)
 
-    def test_reorder_call_listUVVS(self):
+    def old_test_reorder_call_listUVVS(self):
         UVVS = UFOLorentz(name = 'UVVS',
                spins = [ 2, 3, 3, 1])
     
@@ -3144,7 +3445,7 @@
         """ test that python writer works """
         
         solution ="""import wavefunctions
-def SSS1_1(S2, S3, C, M1, W1):
+def SSS1_1(S2, S3, COUP, M1, W1):
     S1 = wavefunctions.WaveFunction(size=3)
     S1[1] = S2[1]+S3[1]
     S1[2] = S2[2]+S3[2]
@@ -3153,12 +3454,17 @@
             - complex(S1[2]).imag, \\
             - complex(S1[1]).imag]
     denom =1.0/(( (M1*( -M1+1j*W1))+( (P1[0]**2)-(P1[1]**2)-(P1[2]**2)-(P1[3]**2))))
-    S1[0]= C*denom*1j*(S3[0]*S2[0])
+    S1[0]= COUP*denom*1j*(S3[0]*S2[0])
     return S1
-def SSS1_2(S2, S3, C, M1, W1):
-    return SSS1_1(S2,S3,C,M1,W1)
-def SSS1_3(S2, S3, C, M1, W1):
-    return SSS1_1(S3,S2,C,M1,W1)"""
+    
+    
+def SSS1_2(S2, S3, COUP, M1, W1):
+    return SSS1_1(S2,S3,COUP,M1,W1)
+
+def SSS1_3(S2, S3, COUP, M1, W1):
+    return SSS1_1(S2,S3,COUP,M1,W1)
+
+"""
         
         SSS = UFOLorentz(name = 'SSS1',
                  spins = [ 1, 1, 1 ],
@@ -3172,8 +3478,7 @@
         
         split_solution = solution.split('\n')
         split_routine = routine.split('\n')
-        for i,line in enumerate(split_routine):
-            self.assertEqual(split_solution[i],line)
+        self.assertEqual(split_solution, split_routine)
         self.assertEqual(len(split_routine), len(split_solution))
         
         
@@ -3191,14 +3496,16 @@
 
         
         solution = """import wavefunctions
-def FFV2C1_0(F1,F2,V3,C):
-    vertex = C*( (F1[3]*( (F2[0]*( -1j*V3[1]-V3[2]))+(F2[1]*( 1j*V3[0]+1j*V3[3]))))+(F1[2]*( (F2[0]*( 1j*V3[0]-1j*V3[3]))+(F2[1]*( -1j*V3[1]+V3[2])))))
-    return vertex""" 
+def FFV2C1_0(F1,F2,V3,COUP):
+    vertex = COUP*( (F1[3]*( (F2[0]*( -1j*V3[1]-V3[2]))+(F2[1]*( 1j*V3[0]+1j*V3[3]))))+(F1[2]*( (F2[0]*( 1j*V3[0]-1j*V3[3]))+(F2[1]*( -1j*V3[1]+V3[2])))))
+    return vertex
+    
+    
+""" 
 
         split_solution = solution.split('\n')
         split_routine = routine.split('\n')
-        for i,line in enumerate(split_routine):
-            self.assertEqual(split_solution[i],line)
+        self.assertEqual(split_solution,split_routine)
         self.assertEqual(len(split_routine), len(split_solution))