Merge lp:~maddevelopers/mg5amcnlo/multi_fermion_vertices into lp:~maddevelopers/mg5amcnlo/combine_v0_6_0

[Johan Alwall]

- Added functionality for multi-fermion vertices in MG5
- Modified order in ALOHA calls to use the order in the UFO interaction (change both in MG5 and in ALOHA)
- Include Modification for Spin2 in ALOHA

[lp:~maddevelopers/mg5amcnlo/multi_fermion_vertices]

71. By Johan Alwall on 2010-12-19

Fixed documentation texts for the t-channel tests in test_4fermion_models.py

72. By Johan Alwall on 2010-12-19

Merged with the latest version of branch combine_v0_6_0

[Olivier Mattelaer (olivier-mattelaer)]

thanks for the merging, this will help a lot (since the diff is more clear now :-) )

[Johan Alwall (johan-alwall)]

??? Can you get on Skype?

On Tue, Dec 21, 2010 at 11:24 AM, Olivier Mattelaer
<email address hidden> wrote:
> thanks for the merging, this will help a lot (since the diff is more clear now :-) )
> --
> https://code.launchpad.net/~maddevelopers/madgraph5/multi_fermion_vertices/+merge/44150
> You proposed lp:~maddevelopers/madgraph5/multi_fermion_vertices for merging.
>

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

Don't use this to contact me. I'm never receiving the email....

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/multi_fermion_vertices]

73. By Johan Alwall on 2010-12-30

Merged with latest trunk 0.6.1

74. By Johan Alwall on 2011-01-13

Fixed configs for multivertices with additional radiation

75. By Johan Alwall on 2011-01-19

Fixed small bug in configs.inc

76. By Johan Alwall on 2011-01-20

Merged with trunk v. 0.6.1.2

77. By Johan Alwall on 2011-01-21

Fixed AMP2 writing for pure multiparticle vertex diagrams, fixed format of couplings

78. By Johan Alwall on 2011-01-21

Updated VERSION

79. By Johan Alwall on 2011-01-21

Added tests for AMP2 expressions also for 4-fermion vertex processes

[Olivier Mattelaer (olivier-mattelaer)]

So Johan,

I think that we can make the merge. What do think?

review: Approve

[lp:~maddevelopers/mg5amcnlo/multi_fermion_vertices]

80. By Olivier Mattelaer on 2011-01-21

update value for the SM coupling estimation

81. By Olivier Mattelaer on 2011-01-26

merge ALOHA_SPIN2 and resolve conflicts (in tests only)

[Johan Alwall (johan-alwall)]

Yeah, I think so too. There is one more thing I would like to fix,
related to the order of decaying particles in decay chains (I have got
one question about this, so others are bound to have the same
problem). I'll fix that in the multi_fermion_vertices brach, and let
you know so you can take a look.

[lp:~maddevelopers/mg5amcnlo/multi_fermion_vertices]

82. By Johan Alwall on 2011-01-27

Allow decay chain specification with order of decays exchanged, e.g. p p > t t~, t~ > w- b~, t > w+ b

=== removed file '.project'
--- .project	2010-08-16 15:41:39 +0000
+++ .project	1970-01-01 00:00:00 +0000
@@ -1,17 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<projectDescription>
-	<name>MG-new_color_stuff</name>
-	<comment></comment>
-	<projects>
-	</projects>
-	<buildSpec>
-		<buildCommand>
-			<name>org.python.pydev.PyDevBuilder</name>
-			<arguments>
-			</arguments>
-		</buildCommand>
-	</buildSpec>
-	<natures>
-		<nature>org.python.pydev.pythonNature</nature>
-	</natures>
-</projectDescription>

=== modified file 'aloha/aloha_lib.py'
--- aloha/aloha_lib.py	2010-11-18 14:21:11 +0000
+++ aloha/aloha_lib.py	2011-01-27 02:53:44 +0000
@@ -591,6 +591,7 @@
         self.prefactor *= obj.prefactor
         if obj in self:
             index = self.index(obj)
+            self[index] = self[index].copy()
             self[index].power += 1
         else:
             obj.prefactor = 1
@@ -1292,6 +1293,7 @@
     def __mul__(self, obj):
         """multiplication performing directly the einstein/spin sommation.
         """
+
         if not hasattr(obj, 'vartype') or not self.vartype:
             out = LorentzObjectRepresentation({}, self.lorentz_ind, self.spin_ind)
             for ind in out.listindices():
@@ -1388,10 +1390,11 @@
                 factor = self.get_rep(self_ind) 
                 factor *= obj.get_rep(obj_ind)
                 
+                
                 if factor:
                     #compute the prefactor due to the lorentz contraction
                     factor *= (-1) ** (len(l_value) - l_value.count(0))
-                    out += factor                         
+                    out += factor                        
         return out
 
     def combine_indices(self, l_dict, s_dict):

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2010-09-12 01:46:50 +0000
+++ aloha/aloha_writers.py	2011-01-27 02:53:44 +0000
@@ -34,8 +34,23 @@
         #prepare the necessary object
         self.collect_variables() # Look for the different variables
         self.make_all_lists()    # Compute the expression for the call ordering
-                                 #the definition of objects,...
-
+                                 #the definition of objects,..
+                                 
+    def pass_to_HELAS(self, indices, start=0):
+        """find the Fortran HELAS position for the list of index""" 
+        
+        
+        if len(indices) == 1:
+            return indices[0] + start
+        
+        ind_name = self.obj.numerator.lorentz_ind 
+        if ind_name == ['I3', 'I2']:
+            return  4 * indices[1] + indices[0] + start 
+        elif len(indices) == 2: 
+            return  4 * indices[0] + indices[1] + start 
+        else:
+            raise Exception                                 
+                                 
     def collect_variables(self):
         """Collects Momenta,Mass,Width into lists"""
          
@@ -146,90 +161,28 @@
         if outgoing is None:
             outgoing = self.offshell
 
-        # particle type counter
-        nb_type = {'S':0, 'F':0, 'V':0, 'T':0}        
         call_arg = [] #incoming argument of the routine
         
-        # update the type counter + make call_arg for amplitude
-        for index,spin in enumerate(self.particles):
-            nb_type[spin] += 1
-            call_arg.append('%s%d' % (spin, index +1))
-            
-        # reorder call_arg if not amplitude
-        if outgoing:
-            part_pos = outgoing -1 
-            out_type = self.particles[part_pos]
-            
-            #order is FVST #look at the border of the cycling move
-            # start/stop are the index of the group of spin where to perform
-            #cycling ordering.
-            if out_type == 'F':
-                start = 0
-                stop = nb_type['F']
-            elif out_type == 'V':
-                start = nb_type['F']
-                stop = start + nb_type['V']
-            elif out_type == 'S':
-                start = nb_type['F'] + nb_type['V']
-                stop = start + nb_type['S']
-            elif out_type == 'T':
-                start = nb_type['F'] + nb_type['V']+ nb_type['S']
-                stop = start + nb_type['T']
-            else:
-                raise NotImplemented, 'Only type FVST are supported' 
-            
-            #reorganize the order and suppress the output from this part
-            call_arg = self.new_order(call_arg, part_pos, start, stop)
-        
+
+        call_arg = ['%s%d' % (spin, index +1) 
+                                 for index,spin in enumerate(self.particles)
+                                 if outgoing != index +1]
+                
         return call_arg
 
     def reorder_call_list(self, call_list, old, new):
-        """ restore the cycling ordering """
+        """ 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 = []
-        #before first same spin
-        for i in range(len(spins)):
-            if spins[i] == spins[old]:
-                start_spin = i
-                break
-            new_call.append(call_list[i])
-        #find end of same spin area
-        stop_spin = len(spins)
-        for i in range(new+1, len(spins)):
-            if spins[i] != spins[old]:
-                stop_spin = i
-                break
-        l_old = range(old + 1, stop_spin) + range(start_spin, old)
-        l_new = range(new + 1, stop_spin) + range(start_spin, new)
-        old_index = l_new.index(old)
-        l_new[old_index] = new
-        for i in range(len(l_old)):
-            part_nb = l_old[l_new.index(l_old[i])] +1
-            new_call.append('%s%d' % (spin, part_nb))
-
-        for i in range(stop_spin, len(spins)):
-            new_call.append(call_list[i-1])           
+        new_call = call_list[:]
+        val = new_call.pop(old)
+        new_call.insert(new - 1, val)
         return new_call
     
-    @ staticmethod
-    def new_order(call_list, remove, start, stop):
-        """ create the new order for the calling using cycling order"""
-        
-        assert(start <= remove <= stop <= len(call_list))
-        
-        new_list= call_list[:start]
-        for i in range(remove+1, stop):
-            new_list.append(call_list[i])
-        for i in range(start, remove):
-            new_list.append(call_list[i])
-        new_list += call_list[stop:]
-        
-        return new_list
         
     def make_momentum_conservation(self):
         """ compute the sign for the momentum conservation """
@@ -444,6 +397,7 @@
             out = '%.9f' % number
             out = self.zero_pattern.sub('', out)
         return out
+        
     
     def define_expression(self):
         OutString = ''
@@ -466,15 +420,13 @@
             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)
             OutString = OutString + string + '\n'
-            counter = 1
             for ind in numerator.listindices():
-                string = '%s(%d)= C*denom*' % (OffShellParticle, counter)
+                string = '%s(%d)= C*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)
                 string = re.sub('(?P<num>[0-9])[Jj]\.', '\g<num>*(0d0,1d0)', string)
                 OutString = OutString + string + '\n' 
-                counter += 1
         return OutString 
     
     def define_symmetry(self, new_nb):
@@ -674,13 +626,11 @@
             string = 'denom =' + '1./(' + denom + ')'
             string = string.replace('+-', '-')
             OutString = OutString + string + ';\n'
-            counter = 0
             for ind in numerator.listindices():
-                string = '%s[%d]= C*denom*' % (OffShellParticle, counter)
+                string = '%s[%d]= C*denom*' % (OffShellParticle, self.pass_to_HELAS(ind))
                 string += self.write_obj(numerator.get_rep(ind))
                 string = string.replace('+-', '-')
                 OutString = OutString + string + ';\n' 
-                counter += 1
         OutString = re.sub('(?P<variable>[A-Za-z]+[0-9]\[*[0-9]*\]*)\*\*(?P<num>[0-9])','pow(\g<variable>,\g<num>)',OutString)
         return OutString 
 
@@ -824,13 +774,11 @@
             string = 'denom =' + '1.0/(' + denom + ')'
             string = string.replace('+-', '-')
             OutString += string + '\n'
-            counter = 0
             for ind in numerator.listindices():
-                string = '%s[%d]= C*denom*' % (self.outname, counter)
+                string = '%s[%d]= C*denom*' % (self.outname, self.pass_to_HELAS(ind))
                 string += self.write_obj(numerator.get_rep(ind))
                 string = string.replace('+-', '-')
                 OutString += string + '\n' 
-                counter += 1
         return OutString 
     
     def define_foot(self):

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2010-10-11 20:32:01 +0000
+++ aloha/create_aloha.py	2011-01-27 02:53:44 +0000
@@ -27,6 +27,7 @@
 sys.path.append(root_path)
 from aloha.aloha_object import *
 import aloha.aloha_writers as aloha_writers
+import aloha.aloha_lib as aloha_lib
 try:
     import madgraph.iolibs.files as files
 except:
@@ -36,6 +37,7 @@
 logger = logging.getLogger('ALOHA')
 
 _conjugate_gap = 50
+_spin2_mult = 1000
 
 class ALOHAERROR(Exception): pass
 
@@ -91,10 +93,10 @@
         self.routine_kernel = None
         
     
-    def compute_routine(self, mode):
+    def compute_routine(self, mode, factorize=True):
         """compute the expression and return it"""
         self.outgoing = mode
-        self.expr = self.compute_aloha_high_kernel(mode)
+        self.expr = self.compute_aloha_high_kernel(mode, factorize)
         return self.define_simple_output()
     
     def define_all_conjugate_builder(self, pair_list):
@@ -115,9 +117,10 @@
         one is conjugated"""
         
         new_builder = copy.copy(self)
+        new_builder.conjg = self.conjg[:]
         try:
             for index in pairs:
-               new_builder.apply_conjugation(index) 
+                new_builder.apply_conjugation(index) 
         except TypeError:
             new_builder.apply_conjugation(pairs) 
         return new_builder
@@ -125,7 +128,8 @@
     def apply_conjugation(self, pair=1):
         """ apply conjugation on self object"""
         
-        old_id = pair
+        
+        old_id = 2 * pair - 1
         new_id = _conjugate_gap + old_id
         
         if self.routine_kernel is None:
@@ -138,7 +142,8 @@
 
         if pair:
             self.name += str(pair)
-        self.conjg.append(pair) 
+        self.conjg.append(pair)
+
     
     
     def define_simple_output(self):
@@ -148,11 +153,10 @@
         return AbstractRoutine(self.expr, self.outgoing, self.spins, self.name, \
                                                                         infostr)
         
-    def compute_aloha_high_kernel(self, mode):
+    def compute_aloha_high_kernel(self, mode, factorize=True):
         """compute the abstract routine associate to this mode """
         # reset tag for particles
         aloha_lib.USE_TAG=set()
-
         #multiply by the wave functions
         nb_spinor = 0
         if not self.routine_kernel:
@@ -171,7 +175,6 @@
             aloha_lib.USE_TAG = set(self.kernel_tag) 
         for (i, spin ) in enumerate(self.spins):
             id = i + 1
-            
             #Check if this is the outgoing particle
             if id == self.outgoing:
                 if spin == 1: 
@@ -182,9 +185,9 @@
                         id += _conjugate_gap
                     nb_spinor += 1
                     if nb_spinor %2:
-                        lorentz *= SpinorPropagator(id, 'I2', i + 1)
+                        lorentz *= SpinorPropagator(id, 'I2', self.outgoing)
                     else:
-                        lorentz *= SpinorPropagator('I2', id, i + 1) 
+                        lorentz *= SpinorPropagator('I2', id, self.outgoing) 
                 elif spin == 3 :
                     lorentz *= VectorPropagator(id, 'I2', id)
                 elif spin == 5 :
@@ -199,13 +202,15 @@
                 elif spin == 2:
                     # shift the tag if we multiply by C matrices
                     if (id+1) // 2 in self.conjg:
-                        id += _conjugate_gap
+                        spin_id = id + _conjugate_gap
+                    else:
+                        spin_id = id
                     nb_spinor += 1
-                    lorentz *= Spinor(id, i + 1)
+                    lorentz *= Spinor(spin_id, id)
                 elif spin == 3:        
                     lorentz *= Vector(id, id)
                 elif spin == 5:
-                    lorentz *= Spin2(10*id+1, 10*id+2, 'I2', 'I3', id)
+                    lorentz *= Spin2(1 * _spin2_mult + id, 2 * _spin2_mult + id, id)
                 else:
                     raise self.AbstractALOHAError(
                                 'The spin value %s is not supported yet' % spin)                    
@@ -221,15 +226,19 @@
           
         #lorentz = lorentz.simplify()
         lorentz = lorentz.expand()
-
-        if self.spins[self.outgoing-1] == 5:
+        if self.outgoing and self.spins[self.outgoing-1] == 5:
             if not self.aloha_lib:
                 AbstractRoutineBuilder.load_library()
             lorentz *= self.aloha_lib[('Spin2Prop', id)]
+            aloha_lib.USE_TAG.add('OM%d' % id)
+            aloha_lib.USE_TAG.add('P%d' % id)
+            
+            
 
         
         lorentz = lorentz.simplify()
-        lorentz = lorentz.factorize()
+        if factorize:
+            lorentz = lorentz.factorize()
         
         lorentz.tag = set(aloha_lib.USE_TAG)
         #raise
@@ -269,8 +278,12 @@
     @classmethod
     def load_library(cls):
     # load the library
-        fsock = open(os.path.join(aloha_path, 'ALOHALib.pkl'), 'r')
-        cls.aloha_lib = cPickle.load(fsock)
+        try:
+            fsock = open(os.path.join(aloha_path, 'ALOHALib.pkl'), 'r')
+        except IOError:
+            cls.aloha_lib = create_library()
+        else:
+            cls.aloha_lib = cPickle.load(fsock)
         
 
 class AbstractALOHAModel(dict):
@@ -638,22 +651,22 @@
     lib = {} # key: (name, part_nb, special) -> object
     for i in range(1, 10):
         logger.info('step %s/9' % i)
-        lib[('Scalar', i)] = create( Scalar(i) )
-        lib[('ScalarProp', i)] = complex(0,1)
-        lib[('Denom', i )] = create( DenominatorPropagator(i) )
-        lib[('Spinor', i )] = create( Spinor(i, i) )
-        lib[('SpinorProp', i, 0)] = create( SpinorPropagator(i, 'I2', i) )
-        lib[('SpinorProp', i, 1)] = create( SpinorPropagator('I2', i, i) )
-        lib[('Vector', i)] = create( Vector(i+1, i+1) )
-        lib[('VectorProp', i)] = create( VectorPropagator(i,'I2', i) )
-        lib[('Spin2', i )] = create( Spin2(10*i+1, 10*i+2, i) )
-        lib[('Spin2Prop',i)] = create( Spin2Propagator(10*i+1, \
-                                            10*i+2,'I2','I3', i) )
-    logger.info('writing')         
-    fsock = open('./ALOHALib.pkl','wb')
+        #lib[('Scalar', i)] = create( Scalar(i) )
+        #lib[('ScalarProp', i)] = complex(0,1)
+        #lib[('Denom', i )] = create( DenominatorPropagator(i) )
+        #lib[('Spinor', i )] = create( Spinor(i, i) )
+        #lib[('SpinorProp', i, 0)] = create( SpinorPropagator(i, 'I2', i) )
+        #lib[('SpinorProp', i, 1)] = create( SpinorPropagator('I2', i, i) )
+        #lib[('Vector', i)] = create( Vector(i+1, i+1) )
+        #lib[('VectorProp', i)] = create( VectorPropagator(i,'I2', i) )
+        #lib[('Spin2', i )] = create( Spin2(10*i+1, 10*i+2, i) )
+        lib[('Spin2Prop',i)] = create( Spin2Propagator(_spin2_mult + i, \
+                                            2 * _spin2_mult + i,'I2','I3', i) )
+    logger.info('writing Spin2 lib')         
+    fsock = open(os.path.join(aloha_path, 'ALOHALib.pkl'),'wb')
     cPickle.dump(lib, fsock, -1)
-    logger.info('done')
-    
+    return lib
+
 if '__main__' == __name__:       
     logging.basicConfig(level=0)
     #create_library()

=== modified file 'aloha/template_files/wavefunctions.py'
--- aloha/template_files/wavefunctions.py	2010-09-05 16:14:23 +0000
+++ aloha/template_files/wavefunctions.py	2011-01-27 02:53:44 +0000
@@ -1,224 +1,355 @@
-import cmath 
+from __future__ import division 
 import math
+from math import sqrt, pow
 
 class WaveFunction(list):
-	"""a objet for a WaveFunction"""
-	
-	spin_to_size={0:1,
-				  1:3,
-				  2:6,
-				  3:6}
-	
-	def __init__(self, spin= None, size=None):
-		"""Init the list with zero value"""
-		
-		if spin:
-			size = self.spin_to_size[spin]
-		list.__init__(self, [0]*size)
-		
+    """a objet for a WaveFunction"""
+    
+    spin_to_size={0:1,
+                  1:3,
+                  2:6,
+                  3:6,
+                  5:18}
+    
+    def __init__(self, spin= None, size=None):
+        """Init the list with zero value"""
+        
+        if spin:
+            size = self.spin_to_size[spin]
+        list.__init__(self, [0]*size)
+        
 
 def ixxxxx(p,fmass,nhel,nsf):
-	"""Defines an inflow fermion."""
-	
-	fi = WaveFunction(2)
-	
-	fi[4] = complex(p[0]*nsf,p[3]*nsf)
-	fi[5] = complex(p[1]*nsf,p[2]*nsf) 
-	
-	nh = nhel*nsf 
-
-	if (fmass != 0.):
-		pp = min(p[0],math.sqrt(math.pow(p[1],2)+math.pow(p[2],2)+math.pow(p[3],2)))
-		if (pp == 0.): 
-			sqm = math.sqrt(abs(fmass)) 
-			ip = (1+nh)/2 
-			im = (1-nh)/2 
-
-			fi[0] = ip*sqm
-			fi[1] = im*nsf*sqm
-			fi[2] = ip*nsf*sqm
-			fi[3] = im*sqm
-
-		else:
-			sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
-			omega = [math.sqrt(p[0]+pp),fmass/(math.sqrt(p[0]+pp))]
-			ip = (1+nh)/2
-			im = (1-nh)/2
-			sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
-			pp3 = max(pp+p[3],0.)
-			if (pp3 == 0.):
-				chi1 = complex(-nh,0.) 
-			else:
-				chi1 = complex(nh*p[1]/math.sqrt(2.*pp*pp3),\
-				p[2]/math.sqrt(2.*pp*pp3))
-			chi = [complex(math.sqrt(pp3*0.5/pp)),chi1]
-
-			fi[0] = sfomeg[0]*chi[im]
-			fi[1] = sfomeg[0]*chi[ip]
-			fi[2] = sfomeg[1]*chi[im]
-			fi[3] = sfomeg[1]*chi[ip] 
-	
-	else: 
-		sqp0p3 = math.sqrt(max(p[0]+p[3],0.))*nsf
-		if (sqp0p3 == 0.):
-			chi1 = complex(-nhel*math.sqrt(2.*p[0]),0.)
-		else:
-			chi1 = complex(nh*p[1]/sqp0p3,p[2]/sqp0p3)
-		chi = [complex(sqp0p3,0.),chi1]
-		if (nh == 1):
-			fi[0] = complex(0.,0.)
-			fi[1] = complex(0.,0.)
-			fi[2] = chi[0]
-			fi[3] = chi[1] 
-		else:
-			fi[0] = chi[1]
-			fi[1] = chi[0]
-			fi[2] = complex(0.,0.)
-			fi[3] = complex(0.,0.) 
-	
-	return fi 
+    """Defines an inflow fermion."""
+    
+    fi = WaveFunction(2)
+    
+    fi[4] = complex(p[0]*nsf,p[3]*nsf)
+    fi[5] = complex(p[1]*nsf,p[2]*nsf) 
+    
+    nh = nhel*nsf 
+
+    if (fmass != 0.):
+        pp = min(p[0],sqrt(p[1]**2 + p[2]**2 + p[3]**2 ))
+        if (pp == 0.): 
+            sqm = sqrt(abs(fmass)) 
+            ip = (1+nh)/2 
+            im = (1-nh)/2 
+
+            fi[0] = ip*sqm
+            fi[1] = im*nsf*sqm
+            fi[2] = ip*nsf*sqm
+            fi[3] = im*sqm
+
+        else:
+            sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
+            omega = [sqrt(p[0]+pp),fmass/(sqrt(p[0]+pp))]
+            ip = (1+nh)//2
+            im = (1-nh)//2
+            sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
+            pp3 = max(pp+p[3],0.)
+            if (pp3 == 0.):
+                chi1 = complex(-nh,0.) 
+            else:
+                chi1 = complex(nh*p[1]/sqrt(2.*pp*pp3),\
+                p[2]/sqrt(2.*pp*pp3))
+            chi = [complex(sqrt(pp3*0.5/pp)),chi1]
+
+            fi[0] = sfomeg[0]*chi[im]
+            fi[1] = sfomeg[0]*chi[ip]
+            fi[2] = sfomeg[1]*chi[im]
+            fi[3] = sfomeg[1]*chi[ip] 
+    
+    else: 
+        sqp0p3 = sqrt(max(p[0]+p[3],0.))*nsf
+        if (sqp0p3 == 0.):
+            chi1 = complex(-nhel*sqrt(2.*p[0]),0.)
+        else:
+            chi1 = complex(nh*p[1]/sqp0p3,p[2]/sqp0p3)
+        chi = [complex(sqp0p3,0.),chi1]
+        if (nh == 1):
+            fi[0] = complex(0.,0.)
+            fi[1] = complex(0.,0.)
+            fi[2] = chi[0]
+            fi[3] = chi[1] 
+        else:
+            fi[0] = chi[1]
+            fi[1] = chi[0]
+            fi[2] = complex(0.,0.)
+            fi[3] = complex(0.,0.) 
+    
+    return fi 
 
 def oxxxxx(p,fmass,nhel,nsf):
-	""" initialize an outgoing fermion"""
-	
-	fo = WaveFunction(2)
-	 	
-	fo[4] = complex(p[0]*nsf,p[3]*nsf)
-	fo[5] = complex(p[1]*nsf,p[2]*nsf)
-
-	nh = nhel*nsf
-
-	if (fmass != 0.):
-		pp = min(p[0],math.sqrt(math.pow(p[1],2)+math.pow(p[2],2)+math.pow(p[3],2)))
-		if (pp == 0.): 
-			sqm = math.sqrt(abs(fmass)) 
-			ip = -((1+nh)/2) 
-			im = (1-nh)/2 
-
-			fo[0] = im*sqm
-			fo[1] = ip*nsf*sqm
-			fo[2] = im*nsf*sqm
-			fo[3] = ip*sqm
-
-		else:
-			sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
-			omega = [math.sqrt(p[0]+pp),fmass/(math.sqrt(p[0]+pp))]
-			ip = (1+nh)/2
-			im = (1-nh)/2
-			sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
-			pp3 = max(pp+p[3],0.)
-			if (pp3 == 0.):
-				chi1 = complex(-nh,0.) 
-			else:
-				chi1 = complex(nh*p[1]/math.sqrt(2.*pp*pp3),\
-				-p[2]/math.sqrt(2.*pp*pp3))
-			chi = [complex(math.sqrt(pp3*0.5/pp)),chi1]
-
-			fo[0] = sfomeg[1]*chi[im]
-			fo[1] = sfomeg[1]*chi[ip]
-			fo[2] = sfomeg[0]*chi[im]
-			fo[3] = sfomeg[0]*chi[ip] 
-			
-	else: 
-		sqp0p3 = math.sqrt(max(p[0]+p[3],0.))*nsf
-		if (sqp0p3 == 0.):
-			chi1 = complex(-nhel*math.sqrt(2.*p[0]),0.)
-		else:
-			chi1 = complex(nh*p[1]/sqp0p3,-p[2]/sqp0p3)
-		chi = [complex(sqp0p3,0.),chi1]
-		if (nh == 1):
-			fo[0] = chi[0]
-			fo[1] = chi[1]
-			fo[2] = complex(0.,0.)
-			fo[3] = complex(0.,0.) 
-		else:
-			fo[0] = complex(0.,0.)
-			fo[1] = complex(0.,0.)
-			fo[2] = chi[1]
-			fo[3] = chi[0] 
-	
-	return fo
+    """ initialize an outgoing fermion"""
+    
+    fo = WaveFunction(2)
+         
+    fo[4] = complex(p[0]*nsf,p[3]*nsf)
+    fo[5] = complex(p[1]*nsf,p[2]*nsf)
+
+    nh = nhel*nsf
+
+    if (fmass != 0.):
+        pp = min(p[0],sqrt(p[1]**2 + p[2]**2 + p[3]**2 ))
+        if (pp == 0.): 
+            sqm = sqrt(abs(fmass)) 
+            ip = -((1+nh)/2) 
+            im = (1-nh)/2 
+
+            fo[0] = im*sqm
+            fo[1] = ip*nsf*sqm
+            fo[2] = im*nsf*sqm
+            fo[3] = ip*sqm
+
+        else:
+            sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
+            omega = [sqrt(p[0]+pp),fmass/(sqrt(p[0]+pp))]
+            ip = (1+nh)//2
+            im = (1-nh)//2
+            sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
+            pp3 = max(pp+p[3],0.)
+            if (pp3 == 0.):
+                chi1 = complex(-nh,0.) 
+            else:
+                chi1 = complex(nh*p[1]/sqrt(2.*pp*pp3),\
+                -p[2]/sqrt(2.*pp*pp3))
+            chi = [complex(sqrt(pp3*0.5/pp)),chi1]
+
+            fo[0] = sfomeg[1]*chi[im]
+            fo[1] = sfomeg[1]*chi[ip]
+            fo[2] = sfomeg[0]*chi[im]
+            fo[3] = sfomeg[0]*chi[ip] 
+            
+    else: 
+        sqp0p3 = sqrt(max(p[0]+p[3],0.))*nsf
+        if (sqp0p3 == 0.):
+            chi1 = complex(-nhel*sqrt(2.*p[0]),0.)
+        else:
+            chi1 = complex(nh*p[1]/sqp0p3,-p[2]/sqp0p3)
+        chi = [complex(sqp0p3,0.),chi1]
+        if (nh == 1):
+            fo[0] = chi[0]
+            fo[1] = chi[1]
+            fo[2] = complex(0.,0.)
+            fo[3] = complex(0.,0.) 
+        else:
+            fo[0] = complex(0.,0.)
+            fo[1] = complex(0.,0.)
+            fo[2] = chi[1]
+            fo[3] = chi[0] 
+    
+    return fo
 
 def vxxxxx(p,vmass,nhel,nsv):
-	""" initialize a vector wavefunction. nhel=4 is for checking BRST"""
-	
-	vc = WaveFunction(3)
-	
-	sqh = math.sqrt(0.5)
-	nsvahl = nsv*abs(nhel)
-	pt2 = math.pow(p[1],2)+math.pow(p[2],2)
-	pp = min(p[0],math.sqrt(pt2+pow(p[3],2)))
-	pt = min(pp,math.sqrt(pt2))
-
-	vc[4] = complex(p[0]*nsv,p[3]*nsv)
-	vc[5] = complex(p[1]*nsv,p[2]*nsv)
-
-	if (nhel == 4):
-		if (vmass == 0.):
-			vc[0] = 1.
-			vc[1]=p[1]/p[0]
-			vc[2]=p[2]/p[0]
-			vc[3]=p[3]/p[0]
-		else:
-			vc[0] = p[0]/vmass
-			vc[1] = p[1]/vmass
-			vc[2] = p[2]/vmass
-			vc[3] = p[3]/vmass
-		
-		return vc 
-
-	if (vmass != 0.):
-		hel0 = 1.-abs(nhel) 
-
-		if (pp == 0.):
-			vc[0] = complex(0.,0.)
-			vc[1] = complex(-nhel*sqh,0.)
-			vc[2] = complex(0.,nsvahl*sqh) 
-			vc[3] = complex(hel0,0.)
-
-		else:
-			emp = p[0]/(vmass*pp)
-			vc[0] = complex(hel0*pp/vmass,0.)
-			vc[3] = complex(hel0*p[3]*emp+nhel*pt/pp*sqh)
-			if (pt != 0.):
-				pzpt = p[3]/(pp*pt)*sqh*nhel
-				vc[1] = complex(hel0*p[1]*emp-p[1]*pzpt, \
-					-nsvahl*p[2]/pt*sqh)
-				vc[2] = complex(hel0*p[2]*emp-p[2]*pzpt, \
-					nsvahl*p[1]/pt*sqh) 
-			else:
-				vc[1] = complex(-nhel*sqh,0.)
-				vc[2] = complex(0.,nsvahl*sign(sqh,p[3]))
-	else: 
-		pp = p[0]
-		pt = math.sqrt(math.pow(p[1],2)+math.pow(p[2],2))
-		vc[0] = complex(0.,0.)
-		vc[3] = complex(nhel*pt/pp*sqh)
-		if (pt != 0.):
-			pzpt = p[3]/(pp*pt)*sqh*nhel
-			vc[1] = complex(-p[1]*pzpt,-nsv*p[2]/pt*sqh)
-			vc[2] = complex(-p[2]*pzpt,nsv*p[1]/pt*sqh)
-		else:
-			vc[1] = complex(-nhel*sqh,0.)
-			vc[2] = complex(0.,nsv*sign(sqh,p[3]))
-	
-	return vc
+    """ initialize a vector wavefunction. nhel=4 is for checking BRST"""
+    
+    vc = WaveFunction(3)
+    
+    sqh = sqrt(0.5)
+    nsvahl = nsv*abs(nhel)
+    pt2 = p[1]**2 + p[2]**2 
+    pp = min(p[0],sqrt(pt2 + p[3]**2))
+    pt = min(pp,sqrt(pt2))
+
+    vc[4] = complex(p[0]*nsv,p[3]*nsv)
+    vc[5] = complex(p[1]*nsv,p[2]*nsv)
+
+    if (nhel == 4):
+        if (vmass == 0.):
+            vc[0] = 1.
+            vc[1]=p[1]/p[0]
+            vc[2]=p[2]/p[0]
+            vc[3]=p[3]/p[0]
+        else:
+            vc[0] = p[0]/vmass
+            vc[1] = p[1]/vmass
+            vc[2] = p[2]/vmass
+            vc[3] = p[3]/vmass
+        
+        return vc 
+
+    if (vmass != 0.):
+        hel0 = 1.-abs(nhel) 
+
+        if (pp == 0.):
+            vc[0] = complex(0.,0.)
+            vc[1] = complex(-nhel*sqh,0.)
+            vc[2] = complex(0.,nsvahl*sqh) 
+            vc[3] = complex(hel0,0.)
+
+        else:
+            emp = p[0]/(vmass*pp)
+            vc[0] = complex(hel0*pp/vmass,0.)
+            vc[3] = complex(hel0*p[3]*emp+nhel*pt/pp*sqh)
+            if (pt != 0.):
+                pzpt = p[3]/(pp*pt)*sqh*nhel
+                vc[1] = complex(hel0*p[1]*emp-p[1]*pzpt, \
+                    -nsvahl*p[2]/pt*sqh)
+                vc[2] = complex(hel0*p[2]*emp-p[2]*pzpt, \
+                    nsvahl*p[1]/pt*sqh) 
+            else:
+                vc[1] = complex(-nhel*sqh,0.)
+                vc[2] = complex(0.,nsvahl*sign(sqh,p[3]))
+    else: 
+        pp = p[0]
+        pt = sqrt(p[1]**2 + p[2]**2)
+        vc[0] = complex(0.,0.)
+        vc[3] = complex(nhel*pt/pp*sqh)
+        if (pt != 0.):
+            pzpt = p[3]/(pp*pt)*sqh*nhel
+            vc[1] = complex(-p[1]*pzpt,-nsv*p[2]/pt*sqh)
+            vc[2] = complex(-p[2]*pzpt,nsv*p[1]/pt*sqh)
+        else:
+            vc[1] = complex(-nhel*sqh,0.)
+            vc[2] = complex(0.,nsv*sign(sqh,p[3]))
+    
+    return vc
 
 def sign(x,y):
-	"""Fortran's sign transfer function"""
-	if (y < 0.):
-		return -abs(x) 
-	else:
-		return abs(x) 
-	
+    """Fortran's sign transfer function"""
+    if (y < 0.):
+        return -abs(x) 
+    else:
+        return abs(x) 
+    
 
 def sxxxxx(p,nss):
-	"""initialize a scalar wavefunction"""
-	
-	sc = WaveFunction(1)
-	
-	sc[0] = complex(1.,0.)
-	sc[1] = complex(p[0]*nss,p[3]*nss)
-	sc[2] = complex(p[1]*nss,p[2]*nss)
-	return sc
+    """initialize a scalar wavefunction"""
+    
+    sc = WaveFunction(1)
+    
+    sc[0] = complex(1.,0.)
+    sc[1] = complex(p[0]*nss,p[3]*nss)
+    sc[2] = complex(p[1]*nss,p[2]*nss)
+    return sc
+
+
+def txxxxx(p, tmass, nhel, nst):
+    """ initialize a tensor wavefunction"""
+    
+    tc = WaveFunction(5)
+    
+    sqh = sqrt(0.5)
+    sqs = sqrt(1/6)
+
+    pt2 = p[1]**2 + p[2]**2
+    pp = min(p[0],sqrt(pt2+p[3]**2))
+    pt = min(pp,sqrt(pt2))
+
+    ft = {}
+    ft[(4,0)] = complex(p[0], p[3]) * nst
+    ft[(5,0)] = complex(p[1], p[2]) * nst
+
+    if ( nhel >= 0 ): 
+        #construct eps+
+        ep = [0] * 4
+        
+        if ( pp == 0 ):
+            #ep[0] = 0
+            ep[1] = -sqh
+            ep[2] = complex(0, nst*sqh)
+            #ep[3] = 0
+        else:            
+            #ep[0] = 0
+            ep[3] = pt/pp*sqh
+            if (pt != 0):
+               pzpt = p[3]/(pp*pt)*sqh
+               ep[1] = complex( -p[1]*pzpt , -nst*p[2]/pt*sqh )
+               ep[2] = complex( -p[2]*pzpt ,  nst*p[1]/pt*sqh )
+            else:
+               ep[1] = -sqh 
+               ep[2] = complex( 0 , nst*sign(sqh,p[3]) )
+            
+         
+     
+    if ( nhel <= 0 ): 
+        #construct eps-
+        em = [0] * 4
+        if ( pp == 0 ):
+            #em[0] = 0
+            em[1] = sqh 
+            em[2] = cpmplex( 0 , nst*sqh )
+            #em[3] = 0
+        else:
+            #em[0] = 0
+            em[3] = -pt/pp*sqh
+            if pt:
+               pzpt = -p[3]/(pp*pt)*sqh
+               em[1] = complex( -p[1]*pzpt , -nst*p[2]/pt*sqh )
+               em[2] = complex( -p[2]*pzpt ,  nst*p[1]/pt*sqh )
+            else:
+               em[1] = sqh
+               em[2] = complex( 0 , nst*sign(sqh,p[3]) )
+            
+    
+    if ( abs(nhel) <= 1 ):  
+        #construct eps0
+        e0 = [0] * 4
+        if ( pp == 0 ):
+            #e0[0] = dcmplx( rZero )
+            #e0[1] = dcmplx( rZero )
+            #e0[2] = dcmplx( rZero )
+            e0[3] = 1
+        else:
+            emp = p[0]/(tmass*pp)
+            e0[0] = pp/tmass 
+            e0[3] = p[3]*emp
+            if pt:
+               e0[1] = p[1]*emp 
+               e0[2] = p[2]*emp 
+            #else:
+            #   e0[1] = dcmplx( rZero )
+            #   e0[2] = dcmplx( rZero )
+
+    if nhel == 2:
+        for j in range(4):
+            for i in range(4):         
+                ft[(i,j)] = ep[i]*ep[j]
+    elif nhel == -2:
+        for j in range(4):
+            for i in range(4):         
+                ft[(i,j)] = em[i]*em[j]
+    elif tmass == 0:
+        for j in range(4):
+            for i in range(4):         
+                ft[(i,j)] = 0
+    elif nhel == 1:
+        for j in range(4):
+            for i in range(4): 
+                 ft[(i,j)] = sqh*( ep[i]*e0[j] + e0[i]*ep[j] )
+    elif nhel == 0:
+        for j in range(4):
+            for i in range(4):       
+                ft[(i,j)] = sqs*( ep[i]*em[j] + em[i]*ep[j] + 2 *e0[i]*e0[j] )
+    elif nhel == -1:
+        for j in range(4):
+            for i in range(4): 
+            	ft[(i,j)] = sqh*( em[i]*e0[j] + e0[i]*em[j] )
+
+    else:
+    	raise Exception, 'invalid helicity TXXXXXX' 
+
+    tc[0] = ft[(0,0)]
+    tc[1] = ft[(0,1)]
+    tc[2] = ft[(0,2)]
+    tc[3] = ft[(0,3)]
+    tc[4] = ft[(1,0)]
+    tc[5] = ft[(1,1)]
+    tc[6] = ft[(1,2)]
+    tc[7] = ft[(1,3)]
+    tc[8] = ft[(2,0)]
+    tc[9] = ft[(2,1)]
+    tc[10] = ft[(2,2)]
+    tc[11] = ft[(2,3)]
+    tc[12] = ft[(3,0)]
+    tc[13] = ft[(3,1)]
+    tc[14] = ft[(3,2)]
+    tc[15] = ft[(3,3)]
+    tc[16] = ft[(4,0)]
+    tc[17] = ft[(5,0)]
+  
+    return tc
+
+    
 

=== modified file 'bin/create_release.py'
--- bin/create_release.py	2010-12-21 14:49:59 +0000
+++ bin/create_release.py	2011-01-27 02:53:44 +0000
@@ -51,12 +51,15 @@
 import optparse
 import logging
 import logging.config
+import time
 
 import os
 import os.path as path
 import shutil
 import subprocess
 
+from datetime import date
+
 # Get the parent directory (mg root) of the script real path (bin)
 # and add it to the current PYTHONPATH
 
@@ -80,14 +83,25 @@
                     format="%(message)s")
 
 # 0. check that all modification are commited in this directory
-#
+#    and that the date is up-to-date
 diff_result = subprocess.Popen(["bzr", "diff"], stdout=subprocess.PIPE).communicate()[0] 
 
 if diff_result:
-    logging.error("Directory is not up-to-date. The release follow the last commited version.")
+    logging.warning("Directory is not up-to-date. The release follow the last commited version.")
     answer = raw_input('Do you want to continue anyway? (y/n)')
     if answer != 'y':
         exit()
+release_date = date.fromtimestamp(time.time())
+for line in file('VERSION'):
+    if 'version' in line:
+        logging.info(line)
+    if 'date' in line:
+        if not str(release_date.year) in line or not str(release_date.month) in line or \
+                                                           not str(release_date.day) in line:
+            logging.warning("WARNING: The release time information is : %s" % line)
+            answer = raw_input('Do you want to continue anyway? (y/n)')
+            if answer != 'y':
+                exit()
 
 # 1. bzr branch the present directory to a new directory
 #    MadGraph5_vVERSION

=== modified file 'bin/mg5'
--- bin/mg5	2010-11-19 21:40:50 +0000
+++ bin/mg5	2011-01-27 02:53:44 +0000
@@ -19,21 +19,21 @@
 and call immediately the command line interface scripts"""
 
 import sys
-
 if not sys.version_info[0] == 2 or sys.version_info[1] < 6:
     sys.exit('MadGraph 5 works only with python 2.6 or later (but not python 3.X).\n\
                Please upgrate your version of python.')
 
+import os
+
 # tab completion and history
 try:
+    import rlcompleter
     import readline
 except ImportError:
     print "For tab completion and history, install module readline."
 else:
-    import rlcompleter
-    readline.parse_and_bind("tab: complete")
+    readline.parse_and_bind("bind ^I rl_complete")
     try:
-        import os
         history_file = os.path.join(os.environ['HOME'], '.pythonhistory')
         readline.read_history_file(history_file)
     except:
@@ -44,13 +44,12 @@
 import logging.config
 
 import subprocess
-import os.path as path
 
 # Get the parent directory (mg root) of the script real path (bin)
 # and add it to the current PYTHONPATH
 
-root_path = path.split(path.dirname(path.realpath( __file__ )))[0]
-sys.path.append(root_path)
+root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0]
+sys.path.insert(0, root_path)
 
 
 # Write out nice usage message if called with -h or --help
@@ -71,17 +70,19 @@
     args = ''
 
 # Check if optimize mode is (and should be) activated
-if __debug__ and not options.debug:
-    if not os.path.exists(os.path.join(root_path, 'bin','create_release.py')):
+if __debug__ and not options.debug and \
+    (not os.path.exists(os.path.join(root_path, 'bin','create_release.py')) or options.web):
         subprocess.call([sys.executable] + ['-O'] + sys.argv)
         sys.exit()
-    else:
+
+if __debug__:
         print 'Running MG5 in debug mode'
 
+
 # Set logging level according to the logging level given by options
 #logging.basicConfig(level=vars(logging)[options.logging])
 try:
-    logging.config.fileConfig(path.join(root_path, 'madgraph', 'interface', '.mg5_logging.conf'))
+    logging.config.fileConfig(os.path.join(root_path, 'madgraph', 'interface', '.mg5_logging.conf'))
     logging.root.setLevel(eval('logging.' + options.logging))
     logging.getLogger('madgraph').setLevel(eval('logging.' + options.logging))
 except:

=== modified file 'madgraph/VERSION'
--- madgraph/VERSION	2010-12-17 01:47:59 +0000
+++ madgraph/VERSION	2011-01-27 02:53:44 +0000
@@ -1,2 +1,2 @@
-version = 0.6.0
-date = 2010-12-17
+version = 0.6.2d
+date = 2011-01-21

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2010-12-15 21:03:07 +0000
+++ madgraph/core/base_objects.py	2011-01-27 02:53:44 +0000
@@ -341,16 +341,22 @@
         if spin == 1:
             # Scalar
             return [ 0 ]
-        if spin == 2:
+        elif spin == 2:
             # Spinor
             return [ -1, 1 ]
-        if spin == 3 and self.get('mass').lower() == 'zero':
+        elif spin == 3 and self.get('mass').lower() == 'zero':
             # Massless vector
             return [ -1, 1 ]
-        if spin == 3:
+        elif spin == 3:
             # Massive vector
             return [ -1, 0, 1 ]
-
+        elif spin == 5 and self.get('mass').lower() == 'zero':
+            # Massless tensor
+            return [-2, -1, 1, 2]
+        elif spin == 5:
+            # Massive tensor
+            return [-2, -1, 0, 1, 2]
+        
         raise self.PhysicsObjectError, \
               "No helicity state assignment for spin %d particles" % spin
 

=== modified file 'madgraph/core/color_amp.py'
--- madgraph/core/color_amp.py	2010-12-19 23:58:53 +0000
+++ madgraph/core/color_amp.py	2011-01-27 02:53:44 +0000
@@ -378,7 +378,7 @@
                 my_cf[0].append(color_algebra.K6Bar(indices[1],
                                                     indices[2],
                                                     indices[3]))
-            if indices[0] == 8:
+            if abs(indices[0]) == 8:
                 # Add a T which contracts the octet to a
                 # triplet-antitriplet pair
                 my_cf[0].append(color_algebra.T(indices[1],
@@ -450,12 +450,12 @@
                 res_dict[leg_num] = [0, 0]
 
                 # Raise an error if external legs contain non supported repr
-                if leg_repr not in [1, 3, -3, 6, -6, 8]:
+                if abs(leg_repr) not in [1, 3, 6, 8]:
                     raise ColorBasis.ColorBasisError, \
         "Particle ID=%i has an unsupported color representation" % leg_repr
 
                 # Build the fake indices replacements for octets
-                if leg_repr == 8:
+                if abs(leg_repr) == 8:
                     fake_repl.append((leg_repr, leg_num,
                                       offset1 + leg_num,
                                       offset2 + leg_num))

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2010-12-17 01:47:59 +0000
+++ madgraph/core/helas_objects.py	2011-01-27 02:53:44 +0000
@@ -443,16 +443,10 @@
             # For boson, set state to intermediate
             self.set('state', 'intermediate')
         else:
-            # For fermion, set state to same as other fermion (in the right way)
-            mother_fermions = filter(lambda wf: wf.is_fermion(),
-                                     self.get('mothers'))
-            if len(mother_fermions) != 1:
-                raise self.PhysicsObjectError, \
-                      """Multifermion vertices not implemented.
-                      Please decompose your vertex into 2-fermion
-                      vertices to get fermion flow correct."""
+            # For fermion, set state to same as other fermion (in the
+            # right way)
+            mother = self.find_mother_fermion()
 
-            mother = mother_fermions[0]
             if self.get('self_antipart'):
                 self.set('state', mother.get_with_flow('state'))
                 self.set('is_part', mother.get_with_flow('is_part'))
@@ -499,12 +493,15 @@
         # Use the HelasWavefunctionList helper function
         # Have to keep track of wavefunction number, since we might
         # need to add new wavefunctions.
-        wf_number = self.get('mothers').check_and_fix_fermion_flow(\
-                                   wavefunctions,
-                                   diagram_wavefunctions,
-                                   external_wavefunctions,
-                                   self.get_with_flow('state'),
-                                   wf_number)
+        self.set('mothers', self.get('mothers').sort_by_pdg_codes(\
+            self.get('pdg_codes'), self.get_anti_pdg_code())[0])
+
+        wf_number = self.get('mothers').\
+                         check_and_fix_fermion_flow(wavefunctions,
+                                                    diagram_wavefunctions,
+                                                    external_wavefunctions,
+                                                    self,
+                                                    wf_number)
 
         return self, wf_number
 
@@ -550,22 +547,14 @@
                 flip_sign = found_majorana
         else:
             # Follow fermion flow up through tree
-            fermion_mother = filter(lambda wf: wf.is_fermion() and
-                                     wf.get_with_flow('state') == \
-                                     self.get_with_flow('state'),
-                                     mothers)
+            fermion_mother = self.find_mother_fermion()
 
-            if len(fermion_mother) > 1:
-                raise self.PhysicsObjectError, \
-                      "6-fermion vertices not yet implemented"
-            if len(fermion_mother) == 0:
-                # Fermion flow already resolved for mothers
-                fermion_mother = filter(lambda wf: wf.is_fermion(),
-                                    mothers)
-                new_mother = fermion_mother[0]
+            if fermion_mother.get_with_flow('state') != \
+                                           self.get_with_flow('state'):
+                new_mother = fermion_mother
             else:
                 # Perform recursion by calling on mother
-                new_mother, wf_number = fermion_mother[0].\
+                new_mother, wf_number = fermion_mother.\
                                         check_majorana_and_flip_flow(\
                                            found_majorana,
                                            wavefunctions,
@@ -586,7 +575,7 @@
                         not self.get('self_antipart')
 
             # Replace old mother with new mother
-            mothers[mothers.index(fermion_mother[0])] = new_mother
+            mothers[mothers.index(fermion_mother)] = new_mother
 
         # Flip sign if needed
         if flip_flow or flip_sign:
@@ -683,49 +672,40 @@
             else:
                 return []
 
-        # Pick out fermion mothers
-        out_fermions = filter(lambda wf: wf.get_with_flow('state') == \
-                              'outgoing', self.get('mothers'))
-        in_fermions = filter(lambda wf: wf.get_with_flow('state') == \
-                             'incoming', self.get('mothers'))
+        # Pick out fermion mother
+        fermion_mother = None
+        if self.is_fermion():
+            fermion_mother = self.find_mother_fermion()
+
+        other_fermions = [wf for wf in self.get('mothers') if \
+                          wf.is_fermion() and wf != fermion_mother]
+
         # Pick out bosons
         bosons = filter(lambda wf: wf.is_boson(), self.get('mothers'))
 
-        if self.is_boson() and len(in_fermions) + len(out_fermions) > 2\
-               or self.is_fermion() and \
-               len(in_fermions) + len(out_fermions) > 1:
-            raise self.PhysicsObjectError, \
-                  "Multifermion vertices not implemented"
-
         fermion_number_list = []
 
-        for boson in bosons:
-            # Bosons return a list [n1,n2,...]
-            fermion_number_list.extend(boson.get_fermion_order())
-
         if self.is_fermion():
             # Fermions return the result N from their mother
             # and the list from bosons, so [N,[n1,n2,...]]
-            fermion_mother = filter(lambda wf: wf.is_fermion(),
-                                    self.get('mothers'))[0]
             mother_list = fermion_mother.get_fermion_order()
             fermion_number_list.extend(mother_list[1])
+
+        # If there are fermion line pairs, append them as
+        # [NI,NO,n1,n2,...]
+        fermion_numbers = [f.get_fermion_order() for f in other_fermions]
+        for iferm in range(0, len(fermion_numbers), 2):
+            fermion_number_list.append(fermion_numbers[iferm][0])
+            fermion_number_list.append(fermion_numbers[iferm+1][0])
+            fermion_number_list.extend(fermion_numbers[iferm][1])
+            fermion_number_list.extend(fermion_numbers[iferm+1][1])
+
+        for boson in bosons:
+            # Bosons return a list [n1,n2,...]
+            fermion_number_list.extend(boson.get_fermion_order())
+
+        if self.is_fermion():
             return [mother_list[0], fermion_number_list]
-        elif in_fermions and out_fermions:
-            # Combine the incoming and outgoing fermion numbers
-            # and add the bosonic numbers: [NI,NO,n1,n2,...]
-            in_list = in_fermions[0].get_fermion_order()
-            out_list = out_fermions[0].get_fermion_order()
-            # Combine to get [N1,N2,n1,n2,...]
-            fermion_number_list.append(in_list[0])
-            fermion_number_list.append(out_list[0])
-            fermion_number_list.extend(in_list[1])
-            fermion_number_list.extend(out_list[1])
-        elif len(in_fermions) != len(out_fermions):
-            raise self.PhysicsObjectError, \
-                  "Error: %d incoming fermions != %d outgoing fermions" % \
-                  (len(in_fermions), len(out_fermions))
-
 
         return fermion_number_list
 
@@ -763,6 +743,18 @@
                state_number[self.get('state')] * \
                self.get('spin')
 
+    def find_mother_fermion(self):
+        """Return the fermion mother which is fermion flow connected to
+        this fermion"""
+
+        if not self.is_fermion():
+            return None
+
+        part_number = self.find_outgoing_number()
+        mother_number = (part_number-1)//2*2
+
+        return HelasMatrixElement.sorted_mothers(self)[mother_number]
+
     def find_outgoing_number(self):
         "Return the position of the resulting particles in the interactions"
         # First shot: just the index in the interaction
@@ -800,7 +792,7 @@
 
         # Check if we need to append a charge conjugation flag
         if self.needs_hermitian_conjugate():
-            res.append('C')
+            res.append(self.get_conjugate_index())
 
         return (tuple(res), self.get('lorentz'))
 
@@ -815,6 +807,8 @@
         if not mothers:
             return vertices
 
+        mothers = self.get('mothers')
+
         # Add vertices for all mothers
         for mother in mothers:
             # This is where recursion happens
@@ -890,7 +884,8 @@
         """Returns two lists of vertices corresponding to the s- and
         t-channels that can be traced from this wavefunction, ordered
         from the outermost s-channel and in/down towards the highest
-        number initial state leg."""
+        number initial state leg. mother_leg corresponds to self but with
+        correct leg number = min(final state mothers)."""
 
         schannels = base_objects.VertexList()
         tchannels = base_objects.VertexList()
@@ -916,15 +911,7 @@
             # leg of a decay process. Add vertex and continue stepping
             # down towards external initial state
             legs = base_objects.LegList()
-
-            if ninitial == 1 or init_mothers[0].get('number_external') == 2 \
-                   or init_mothers[0].get('leg_state') == True:
-                # This is an s-channel or a leg on its way towards the
-                # final vertex
-                mothers = final_mothers + init_mothers
-            else:
-                # This is a t-channel leg going up towards leg number 1
-                mothers = init_mothers + final_mothers
+            mothers = final_mothers + init_mothers
 
             for mother in mothers:
                 legs.append(base_objects.Leg({
@@ -934,19 +921,20 @@
                     'from_group': False
                     }))
 
-            if ninitial == 1 or init_mothers[0].get('number_external') == 2 \
-                   or init_mothers[0].get('leg_state') == True:
-                # For decay processes or if this is an s-channel leg
-                # or we are going towards external leg 2, mother leg
-                # is one of the mothers
+            if init_mothers[0].get('number_external') == 1 and \
+                   not init_mothers[0].get('leg_state') and \
+                   ninitial > 1:
+                # If this is t-channel going towards external leg 1,
+                # mother_leg is resulting wf
+                legs.append(mother_leg)
+            else:
+                # For decay processes or if init_mother is an s-channel leg
+                # or we are going towards external leg 2, mother_leg
+                # is one of the mothers (placed next-to-last)
                 legs.insert(-1, mother_leg)
-                # Also need to switch direction of the resulting s-channel
+                # Need to switch direction of the resulting s-channel
                 legs[-1].set('id', init_mothers[0].get_anti_pdg_code())
-            else:
-                # If the mother is going
-                # towards external leg 1, mother leg is resulting wf
-                legs.append(mother_leg)
-
+                
             # Renumber resulting leg according to minimum leg number
             legs[-1].set('number', min([l.get('number') for l in legs[:-1]]))
 
@@ -954,30 +942,32 @@
                 'id': self.get('interaction_id'),
                 'legs': legs})
 
-            # Add s- and t-channels from further down
+            # Add s- and t-channels from init_mother
             new_mother_leg = legs[-1]
-            if init_mothers[0].get('number_external') == 1:
-                # If we are going towards external leg 1, mother of
-                # next vertex is legs[0]
-                new_mother_leg = legs[0]
+            if init_mothers[0].get('number_external') == 1 and \
+                   not init_mothers[0].get('leg_state') and \
+                   ninitial > 1:
+                # Mother of next vertex is init_mothers[0]
+                # (next-to-last in legs)
+                new_mother_leg = legs[-2]
 
             mother_s, tchannels = \
                       init_mothers[0].get_s_and_t_channels(ninitial,
                                                            new_mother_leg)
-            schannels.extend(mother_s)
-
             if ninitial == 1 or init_mothers[0].get('leg_state') == True:
-                # This leg is s-channel
+                # This vertex is s-channel
                 schannels.append(vertex)
             elif init_mothers[0].get('number_external') == 1:
-                # If we are going towards external leg 1, add to
-                # t-channels, at end
+                # If init_mothers is going towards external leg 1, add
+                # to t-channels, at end
                 tchannels.append(vertex)
             else:
-                # If we are going towards external leg 2, add to
+                # If init_mothers is going towards external leg 2, add to
                 # t-channels, at start
                 tchannels.insert(0, vertex)
 
+            schannels.extend(mother_s)
+
         elif len(init_mothers) == 2:
             # This is a t-channel junction. Start with the leg going
             # towards external particle 1, and then do external
@@ -996,7 +986,7 @@
                     'state': mother.get('leg_state'),
                     'from_group': False
                     }))
-            legs.insert(-1, mother_leg)
+            legs.insert(0, mother_leg)
 
             # Renumber resulting leg according to minimum leg number
             legs[-1].set('number', min([l.get('number') for l in legs[:-1]]))
@@ -1007,7 +997,7 @@
 
             # Add s- and t-channels going down towards leg 1
             mother_s, tchannels = \
-                      init_mothers1.get_s_and_t_channels(ninitial, legs[0])
+                      init_mothers1.get_s_and_t_channels(ninitial, legs[1])
             schannels.extend(mother_s)
 
             # Add vertex
@@ -1019,6 +1009,7 @@
             schannels.extend(mother_s)
             tchannels.extend(mother_t)
 
+
         return schannels, tchannels
 
     def get_conjugate_index(self):
@@ -1039,6 +1030,18 @@
         else:
             return ()
 
+    def get_vertex_leg_numbers(self):
+        """Get a list of the number of legs in vertices in this diagram"""
+
+        if not self.get('mothers'):
+            return []
+
+        vertex_leg_numbers = [len(self.get('mothers')) + 1]
+        for mother in self.get('mothers'):
+            vertex_leg_numbers.extend(mother.get_vertex_leg_numbers())
+
+        return vertex_leg_numbers
+
     # Overloaded operators
 
     def __eq__(self, other):
@@ -1100,7 +1103,7 @@
                                    wavefunctions,
                                    diagram_wavefunctions,
                                    external_wavefunctions,
-                                   my_state,
+                                   my_wf,
                                    wf_number,
                                    force_flip_flow=False,
                                    number_to_wavefunctions=[]):
@@ -1117,90 +1120,85 @@
         wavefunction number.
         """
 
-        # Clash is defined by whether the mothers have N(incoming) !=
-        # N(outgoing) after this state has been subtracted
-        mother_states = [ wf.get_with_flow('state') for wf in \
-                          self ]
-        try:
-            mother_states.remove(my_state)
-        except ValueError:
-            pass
-
-        Nincoming = len(filter(lambda state: state == 'incoming',
-                               mother_states))
-        Noutgoing = len(filter(lambda state: state == 'outgoing',
-                               mother_states))
-
-        if Nincoming == Noutgoing:
+        # Clash is defined by whether any of the fermion lines are clashing
+        fermion_mother = None
+
+        # Keep track of clashing fermion wavefunctions
+        clashes = []
+        
+        # First check the fermion mother on the same fermion line
+        if my_wf and my_wf.is_fermion():
+            fermion_mother = my_wf.find_mother_fermion()
+            if my_wf.get_with_flow('state') != \
+                   fermion_mother.get_with_flow('state'):
+                clashes.append([fermion_mother])
+
+        # Now check all other fermions
+        other_fermions = [w for w in self if \
+                          w.is_fermion() and w != fermion_mother]
+
+        for iferm in range(0, len(other_fermions), 2):
+            if other_fermions[iferm].get_with_flow('state') == \
+               other_fermions[iferm+1].get_with_flow('state'):
+                clashes.append([other_fermions[iferm],
+                                other_fermions[iferm+1]])
+
+        if not clashes:
             return wf_number
 
-        fermion_mothers = filter(lambda wf: wf.is_fermion(),
-                                 self)
-        if my_state in ['incoming', 'outgoing'] and len(fermion_mothers) > 1\
-           or len(fermion_mothers) > 2:
-            raise self.PhysicsObjectListError, \
-                      """Multifermion vertices not implemented.
-                      Please decompose your vertex into 2-fermion
-                      vertices to get fermion flow correct."""
-
         # If any of the mothers have negative fermionflow, we need to
         # take this mother first.
-        neg_fermionflow_mothers = [mother for mother in fermion_mothers if \
-                                    mother.get('fermionflow') < 0]
-
-        if not neg_fermionflow_mothers:
-            neg_fermionflow_mothers = fermion_mothers
-
-        for mother in neg_fermionflow_mothers:
-            if Nincoming > Noutgoing and \
-               mother.get_with_flow('state') == 'outgoing' or \
-               Nincoming < Noutgoing and \
-               mother.get_with_flow('state') == 'incoming' or \
-               Nincoming == Noutgoing:
-                # This is not a problematic leg
-                continue
-
-            # Call recursive function to check for Majorana fermions
-            # and flip fermionflow if found
-
-            found_majorana = False
-            new_mother, wf_number = mother.check_majorana_and_flip_flow(\
-                                                found_majorana,
-                                                wavefunctions,
-                                                diagram_wavefunctions,
-                                                external_wavefunctions,
-                                                wf_number,
-                                                force_flip_flow,
-                                                number_to_wavefunctions)
-
-            # Replace old mother with new mother
-            self[self.index(mother)] = new_mother
-
-            # Update counters
-            mother_states = [ wf.get_with_flow('state') for wf in \
-                             self ]
-            try:
-                mother_states.remove(my_state)
-            except ValueError:
-                pass
-
-            Nincoming = len(filter(lambda state: state == 'incoming',
-                                       mother_states))
-            Noutgoing = len(filter(lambda state: state == 'outgoing',
-                                       mother_states))
-
-        if Nincoming != Noutgoing:
-            # No Majorana fermion in any relevant legs - try again,
-            # but simply use the first relevant leg
-            force_flip_flow = True
-            wf_number = self.check_and_fix_fermion_flow(\
-                                   wavefunctions,
-                                   diagram_wavefunctions,
-                                   external_wavefunctions,
-                                   my_state,
-                                   wf_number,
-                                   force_flip_flow,
-                                   number_to_wavefunctions)
+        for clash in clashes:
+            neg_fermionflow_mothers = [m for m in clash if \
+                                       m.get('fermionflow') < 0]
+
+            if not neg_fermionflow_mothers:
+                neg_fermionflow_mothers = clash
+
+            for mother in neg_fermionflow_mothers:
+
+                # Call recursive function to check for Majorana fermions
+                # and flip fermionflow if found
+
+                found_majorana = False
+                state_before = mother.get_with_flow('state')
+                new_mother, wf_number = mother.check_majorana_and_flip_flow(\
+                                                    found_majorana,
+                                                    wavefunctions,
+                                                    diagram_wavefunctions,
+                                                    external_wavefunctions,
+                                                    wf_number,
+                                                    force_flip_flow,
+                                                    number_to_wavefunctions)
+
+                if new_mother.get_with_flow('state') == state_before:
+                    # Fermion flow was not flipped, try next mother
+                    continue
+
+                # Replace old mother with new mother
+                mother_index = self.index(mother)
+                self[self.index(mother)] = new_mother
+                clash_index = clash.index(mother)
+                clash[clash.index(mother)] = new_mother
+
+                # Fermion flow was flipped, abort loop
+                break
+            
+            if len(clash) == 1 and clash[0].get_with_flow('state') != \
+                   my_wf.get_with_flow('state') or \
+                   len(clash) == 2 and clash[0].get_with_flow('state') == \
+                   clash[1].get_with_flow('state'):
+                # No Majorana fermion in any relevant legs - try again,
+                # but simply use the first relevant leg
+                force_flip_flow = True
+                wf_number = self.check_and_fix_fermion_flow(\
+                                       wavefunctions,
+                                       diagram_wavefunctions,
+                                       external_wavefunctions,
+                                       my_wf,
+                                       wf_number,
+                                       force_flip_flow,
+                                       number_to_wavefunctions)
 
         return wf_number
 
@@ -1231,13 +1229,16 @@
         order of the pdg codes given. my_pdg_code is the pdg code of
         the daughter wavefunction (or 0 if daughter is amplitude)."""
 
+        if not pdg_codes:
+            return self, 0
+
         pdg_codes = copy.copy(pdg_codes)
 
         # Remove the argument wavefunction code from pdg_codes
 
         my_index = -1
         if my_pdg_code:
-            # Find the last index instead of the first, to work with UFO models
+            # Remember index of my code
             my_index = pdg_codes.index(my_pdg_code)
             pdg_codes.pop(my_index)
         
@@ -1262,7 +1263,7 @@
         if mothers:
             raise base_objects.PhysicsObject.PhysicsObjectError
 
-        return sorted_mothers, my_index
+        return HelasWavefunctionList(sorted_mothers), my_index
 
     @staticmethod
     def extract_wavefunctions(mothers):
@@ -1482,11 +1483,14 @@
         HelasWavefunction.check_and_fix_fermion_flow.
         """
 
+        self.set('mothers', self.get('mothers').sort_by_pdg_codes(\
+            self.get('pdg_codes'), 0)[0])
+                 
         return self.get('mothers').check_and_fix_fermion_flow(\
                                    wavefunctions,
                                    diagram_wavefunctions,
                                    external_wavefunctions,
-                                   "Nostate",
+                                   None,
                                    wf_number)
 
 
@@ -1510,7 +1514,7 @@
 
         # Check if we need to append a charge conjugation flag
         if self.needs_hermitian_conjugate():
-            res.append('C')
+            res.append(self.get_conjugate_index())
 
         return (tuple(res), self.get('lorentz'))
 
@@ -1519,38 +1523,26 @@
         to this amplitude"""
 
         # Pick out fermion mothers
-        out_fermions = filter(lambda wf: wf.get_with_flow('state') == \
-                              'outgoing', self.get('mothers'))
-        in_fermions = filter(lambda wf: wf.get_with_flow('state') == \
-                             'incoming', self.get('mothers'))
+        fermions = [wf for wf in self.get('mothers') if wf.is_fermion()]
+
         # Pick out bosons
         bosons = filter(lambda wf: wf.is_boson(), self.get('mothers'))
 
-        if len(in_fermions) + len(out_fermions) > 2:
-            raise self.PhysicsObjectError, \
-                  "Multifermion vertices not implemented"
-
         fermion_number_list = []
 
+        # If there are fermion line pairs, append them as
+        # [NI,NO,n1,n2,...]
+        fermion_numbers = [f.get_fermion_order() for f in fermions]
+        for iferm in range(0, len(fermion_numbers), 2):
+            fermion_number_list.append(fermion_numbers[iferm][0])
+            fermion_number_list.append(fermion_numbers[iferm+1][0])
+            fermion_number_list.extend(fermion_numbers[iferm][1])
+            fermion_number_list.extend(fermion_numbers[iferm+1][1])
+
         for boson in bosons:
             # Bosons return a list [n1,n2,...]
             fermion_number_list.extend(boson.get_fermion_order())
 
-        if in_fermions and out_fermions:
-            # Fermions return the result N from their mother
-            # and the list from bosons, so [N,[n1,n2,...]]
-            in_list = in_fermions[0].get_fermion_order()
-            out_list = out_fermions[0].get_fermion_order()
-            # Combine to get [N1,N2,n1,n2,...]
-            fermion_number_list.append(in_list[0])
-            fermion_number_list.append(out_list[0])
-            fermion_number_list.extend(in_list[1])
-            fermion_number_list.extend(out_list[1])
-        elif len(in_fermions) != len(out_fermions):
-            raise self.PhysicsObjectError, \
-                  "Error: %d incoming fermions != %d outgoing fermions" % \
-                  (len(in_fermions), len(out_fermions))
-
         self['fermionfactor'] = self.sign_flips_to_order(fermion_number_list)
 
     def sign_flips_to_order(self, fermions):
@@ -1671,17 +1663,13 @@
 
             # Create vertex
             legs = base_objects.LegList()
-            for mother in [init_mothers1] + final_mothers + [init_mothers2]:
+            for mother in final_mothers + [init_mothers1] + [init_mothers2]:
                 legs.append(base_objects.Leg({
                     'id': mother.get_pdg_code(),
                     'number': mother.get('number_external'),
                     'state': mother.get('leg_state'),
                     'from_group': False
                     }))
-            # Sort the legs
-            legs = base_objects.LegList(\
-                   sorted(legs[:-1], lambda l1, l2: l1.get('number') - \
-                          l2.get('number')) + legs[-1:])
             # Renumber resulting leg according to minimum leg number
             legs[-1].set('number', min([l.get('number') for l in legs[:-1]]))
 
@@ -1691,7 +1679,7 @@
 
             # Add s- and t-channels going down towards leg 1
             mother_s, tchannels = \
-                      init_mothers1.get_s_and_t_channels(ninitial, legs[0])
+                      init_mothers1.get_s_and_t_channels(ninitial, legs[-2])
 
             schannels.extend(mother_s)
 
@@ -1704,6 +1692,41 @@
             schannels.extend(mother_s)
             tchannels.extend(mother_t)
 
+        # Split up multiparticle vertices using fake s-channel propagators
+        multischannels = [(i, v) for (i, v) in enumerate(schannels) \
+                          if len(v.get('legs')) > 3]
+        multitchannels = [(i, v) for (i, v) in enumerate(tchannels) \
+                          if len(v.get('legs')) > 3]
+
+        for channel in multischannels + multitchannels:
+            newschannels = []
+            vertex = channel[1]
+            while len(vertex.get('legs')) > 3:
+                # Pop the first two legs and create a new
+                # s-channel from them
+                popped_legs = \
+                           base_objects.LegList([vertex.get('legs').pop(0) \
+                                                    for i in [0,1]])
+                popped_legs.append(base_objects.Leg({'id': 21,
+                    'number': min([l.get('number') for l in popped_legs]),
+                    'state': True,
+                    'from_group': False}))
+
+                new_vertex = base_objects.Vertex({
+                    'id': vertex.get('id'),
+                    'legs': popped_legs})
+
+                # Insert the new s-channel before this vertex
+                if channel in multischannels:
+                    schannels.insert(channel[0], new_vertex)
+                else:
+                    schannels.append(new_vertex)
+                legs = vertex.get('legs')
+                # Insert the new s-channel into vertex
+                legs.insert(0, copy.copy(popped_legs[-1]))
+                # Renumber resulting leg according to minimum leg number
+                legs[-1].set('number', min([l.get('number') for l in legs[:-1]]))
+
         # Finally go through all vertices, sort the legs and replace
         # leg number with propagator number -1, -2, ...
         number_dict = {}
@@ -1772,6 +1795,15 @@
         else:
             return ()
 
+    def get_vertex_leg_numbers(self):
+        """Get a list of the number of legs in vertices in this diagram"""
+
+        vertex_leg_numbers = [len(self.get('mothers'))]
+        for mother in self.get('mothers'):
+            vertex_leg_numbers.extend(mother.get_vertex_leg_numbers())
+
+        return vertex_leg_numbers
+
     def set_coupling_color_factor(self):
         """Check if there is a mismatch between order of fermions
         w.r.t. color"""
@@ -1889,6 +1921,11 @@
 
         return coupling_orders
 
+    def get_vertex_leg_numbers(self):
+        """Get a list of the number of legs in vertices in this diagram"""
+
+        return self.get('amplitudes')[0].get_vertex_leg_numbers()
+
 
 #===============================================================================
 # HelasDiagramList
@@ -2236,31 +2273,37 @@
             helas_diagram.set('number', diagram_number)
             for number_wf_dict, color_list in zip(number_to_wavefunctions,
                                                   color_lists):
+
+                # Now generate HelasAmplitudes from the last vertex.
+                if lastvx.get('id'):
+                    inter = model.get_interaction(lastvx.get('id'))
+                    keys = sorted(inter.get('couplings').keys())
+                    pdg_codes = [p.get_pdg_code() for p in \
+                                 inter.get('particles')]
+                else:
+                    # Special case for decay chain - amplitude is just a
+                    # placeholder for replaced wavefunction
+                    inter = None
+                    keys = [(0, 0)]
+                    pdg_codes = None
+
                 # Find mothers for the amplitude
                 legs = lastvx.get('legs')
                 mothers = self.getmothers(legs, number_wf_dict,
                                           external_wavefunctions,
                                           wavefunctions,
-                                          diagram_wavefunctions)
+                                          diagram_wavefunctions).\
+                                             sort_by_pdg_codes(pdg_codes, 0)[0]
                 # Need to check for clashing fermion flow due to
                 # Majorana fermions, and modify if necessary
                 wf_number = mothers.check_and_fix_fermion_flow(wavefunctions,
                                               diagram_wavefunctions,
                                               external_wavefunctions,
-                                              "Nostate",
+                                              None,
                                               wf_number,
                                               False,
                                               number_to_wavefunctions)
 
-                # Now generate HelasAmplitudes from the last vertex.
-                if lastvx.get('id'):
-                    inter = model.get('interaction_dict')[lastvx.get('id')]
-                    keys = sorted(inter.get('couplings').keys())
-                else:
-                    # Special case for decay chain - amplitude is just a
-                    # placeholder for replaced wavefunction
-                    inter = None
-                    keys = [(0, 0)]
                 for i, coupl_key in enumerate(keys):
                     amp = HelasAmplitude(lastvx, model)
                     if inter:
@@ -2340,7 +2383,8 @@
         for wf in self.get_all_wavefunctions() + \
             sum([d.get_all_wavefunctions() for d in \
                  decay_dict.values()], []):
-            if wf.get('self_antipart') and wf.is_fermion():
+            if wf.get('fermionflow') < 0 or \
+                   wf.get('self_antipart') and wf.is_fermion():
                 got_majoranas = True
 
         # Now insert decays for all legs that have decays
@@ -2533,13 +2577,6 @@
                 numbers[0] = numbers[0] + 1
                 wf.set('number', numbers[0])
 
-            # External wavefunction offset for new wfs
-            incr_new = number_external - \
-                       decay_wfs[0].get('number_external')
-            for wf in decay_wfs:
-                wf.set('number_external', wf.get('number_external') + incr_new)
-
-
         # External wavefunction offset for old wfs, only the first
         # time this external wavefunction is replaced
         (nex, nin) = decay.get_nexternal_ninitial()
@@ -2752,7 +2789,7 @@
                     for wf in final_decay_wfs:
                         wf.set('onshell', True)
 
-                    if len_decay == 1:
+                    if len_decay == 1 and len(final_decay_wfs) == 1:
                         # Can use simplified treatment, by just modifying old_wf
                         self.replace_single_wavefunction(old_wf,
                                                          final_decay_wfs[0])
@@ -3378,7 +3415,7 @@
     def sorted_mothers(arg):
         """Gives a list of mother wavefunctions sorted according to
         1. The order of the particles in the interaction
-        2. Cyclic reordering of particles in same spin group (if boson)
+        2. Cyclic reordering of particles in same spin group
         3. Fermions ordered IOIOIO... according to the pairs in
            the interaction."""
 
@@ -3388,7 +3425,6 @@
         if not arg.get('interaction_id'):
             return arg.get('mothers')
 
-        sorted_mothers = copy.copy(arg.get('mothers'))
         my_pdg_code = 0
         my_spin = 0
         if isinstance(arg, HelasWavefunction):
@@ -3398,32 +3434,9 @@
         sorted_mothers, my_index = arg.get('mothers').sort_by_pdg_codes(\
             arg.get('pdg_codes'), my_pdg_code)
 
-        same_spin_mothers = []
-        if isinstance(arg, HelasWavefunction) and \
-               my_spin % 2 == 1:
-            # For bosons with same spin as this wf, need special treatment
-            same_spin_index = -1
-            i=0
-            while i < len(sorted_mothers):
-                if sorted_mothers[i].get_spin_state_number() == my_spin:
-                    if same_spin_index < 0:
-                        # Remember starting index for same spin states
-                        same_spin_index = i
-                    same_spin_mothers.append(sorted_mothers.pop(i))
-                else:
-                    i += 1
-
-        # Make cyclic reordering of mothers with same spin as this wf
-        if same_spin_mothers:
-            same_spin_mothers = same_spin_mothers[my_index - same_spin_index:] \
-                                + same_spin_mothers[:my_index - same_spin_index]
-
-            # Insert same_spin_mothers in sorted_mothers
-            sorted_mothers = sorted_mothers[:same_spin_index] + \
-                              same_spin_mothers + sorted_mothers[same_spin_index:]
         # If fermion, partner is the corresponding fermion flow partner
         partner = None
-        if isinstance(arg, HelasWavefunction) and my_spin % 2 == 0:
+        if isinstance(arg, HelasWavefunction) and arg.is_fermion():
             # Fermion case, just pick out the fermion flow partner
             if my_index % 2 == 0:
                 # partner is after arg
@@ -3432,10 +3445,15 @@
                 # partner is before arg
                 partner_index = my_index - 1
             partner = sorted_mothers.pop(partner_index)
+            # If partner is incoming, move to before arg
+            if partner.get_spin_state_number() > 0:
+                my_index = partner_index
+            else:
+                my_index = partner_index + 1
 
         # Reorder fermions pairwise according to incoming/outgoing
         for i in range(0, len(sorted_mothers), 2):
-            if sorted_mothers[i].get_spin_state_number() % 2 == 0:
+            if sorted_mothers[i].is_fermion():
                 # This is a fermion, order between this fermion and its brother
                 if sorted_mothers[i].get_spin_state_number() > 0 and \
                    sorted_mothers[i + 1].get_spin_state_number() < 0:
@@ -3447,12 +3465,6 @@
                    sorted_mothers[i + 1].get_spin_state_number() > 0:
                     # This is the right order
                     pass
-                else:
-                    # Two incoming or two outgoing in a row - not good!
-                    raise base_objects.PhysicsObject.PhysicsObjectError, \
-                    "Two incoming or outgoing fermions in a row: %i, %i" % \
-                    (sorted_mothers[i].get_spin_state_number(),
-                     sorted_mothers[i+1].get_spin_state_number())
             else:
                 # No more fermions in sorted_mothers
                 break
@@ -3652,8 +3664,11 @@
                                              me.get('processes')[0].\
                                              get_initial_ids()[0] == fs_id,
                                              decay_elements[index]))
-                else:
-                    # All decays for this particle type are used
+
+                if len(fs_legs) != len(decay_elements) or not chains[0]:
+                    # In second case, or no chains are found
+                    # (e.g. because the order of decays is reversed),
+                    # all decays for this particle type are used
                     chain = sum([filter(lambda me: \
                                         me.get('processes')[0].\
                                         get_initial_ids()[0] == fs_id,

=== modified file 'madgraph/interface/cmd_interface.py'
--- madgraph/interface/cmd_interface.py	2010-12-18 22:39:25 +0000
+++ madgraph/interface/cmd_interface.py	2011-01-27 02:53:44 +0000
@@ -197,6 +197,13 @@
         if '#' in line:
             line = line.split('#')[0]
 
+        # Deal with line splitting
+        if ';' in line:
+            for subline in line.split(';'):
+                stop = self.onecmd(subline)
+                stop = self.postcmd(stop, subline)
+            return ''
+        
         # execute the line command
         return line
 
@@ -1304,7 +1311,8 @@
 
     # Options and formats available
     _display_opts = ['particles', 'interactions', 'processes', 'diagrams', 
-                     'multiparticles', 'couplings', 'lorentz', 'checks']
+                     'multiparticles', 'couplings', 'lorentz', 'checks',
+                     'parameters']
     _add_opts = ['process']
     _save_opts = ['model', 'processes']
     _tutorial_opts = ['start', 'stop']
@@ -1510,7 +1518,26 @@
                     print "Interactions %s has the following property:" % arg
                     print self._curr_model['interactions'][int(arg)-1]
 
-
+        elif args[0] == 'parameters' and len(args) == 1:
+            text = "Current model contains %i parameters\n" % \
+                    sum([len(part) for part in 
+                                       self._curr_model['parameters'].values()])
+            
+            for key, item in self._curr_model['parameters'].items():
+                text += '\nparameter type: %s\n' % str(key)
+                for value in item:
+                    if hasattr(value, 'expr'):
+                        if value.value is not None:
+                            text+= '        %s = %s = %s\n' % (value.name, value.expr ,value.value)
+                        else:
+                            text+= '        %s = %s\n' % (value.name, value.expr)
+                    else:
+                        if value.value is not None:
+                            text+= '        %s = %s\n' % (value.name, value.value)
+                        else:
+                            text+= '        %s \n' % (value.name)
+            pydoc.pager(text)
+            
         elif args[0] == 'processes':
             for amp in self._curr_amps:
                 print amp.nice_string_processes()

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2010-12-20 05:04:27 +0000
+++ madgraph/iolibs/export_v4.py	2011-01-27 02:53:44 +0000
@@ -197,6 +197,11 @@
     """Finalize ME v4 directory by creating jpeg diagrams, html
     pages,proc_card_mg5.dat and madevent.tar.gz."""
 
+    if not misc.which('f77'):
+        logger.info('Change makefiles to use gfortran')
+        subprocess.call(['python','./bin/Passto_gfortran.py'], cwd=dir_path, \
+                        stdout = open(os.devnull, 'w')) 
+    
     old_pos = os.getcwd()
     os.chdir(os.path.join(dir_path, 'SubProcesses'))
     P_dir_list = [proc for proc in os.listdir('.') if os.path.isdir(proc) and \
@@ -408,6 +413,10 @@
                 matrix_element)
     replace_dict['helas_calls'] = "\n".join(helas_calls)
 
+    # Extract amp2 lines
+    amp2_lines = get_amp2_lines(matrix_element)
+    replace_dict['amp2_lines'] = '\n'.join(amp2_lines)
+
     # Extract JAMP lines
     jamp_lines = get_JAMP_lines(matrix_element)
     replace_dict['jamp_lines'] = '\n'.join(jamp_lines)
@@ -495,18 +504,22 @@
 
     s_and_t_channels = []
 
-    for idiag, diag in enumerate(matrix_element.get('base_amplitude').\
-                                                get('diagrams')):
-        if any([len(vert.get('legs')) > 3 for vert in diag.get('vertices')]):
+    base_diagrams = matrix_element.get('base_amplitude').get('diagrams')
+    minvert = min([max([len(vert.get('legs')) for vert in \
+                        diag.get('vertices')]) for diag in base_diagrams])
+
+    for idiag, diag in enumerate(base_diagrams):
+        if any([len(vert.get('legs')) > minvert for vert in
+                diag.get('vertices')]):
             # Only 3-vertices allowed in configs.inc
             continue
         iconfig = iconfig + 1
         helas_diag = matrix_element.get('diagrams')[idiag]
-        amp_number = helas_diag.get('amplitudes')[0].get('number')
-        lines.append("# Diagram %d, Amplitude %d" % \
-                     (helas_diag.get('number'), amp_number))
+        lines.append("# Diagram %d" % \
+                     (helas_diag.get('number')))
         # Correspondance between the config and the amplitudes
-        lines.append("data mapconfig(%d)/%d/" % (iconfig, amp_number))
+        lines.append("data mapconfig(%d)/%d/" % (iconfig,
+                                                 helas_diag.get('number')))
 
         # Need to reorganize the topology so that we start with all
         # final state external particles and work our way inwards
@@ -635,7 +648,8 @@
                 repr_dict = {}
                 for l in legs:
                     repr_dict[l.get('number')] = \
-                        proc.get('model').get_particle(l.get('id')).get_color()
+                        proc.get('model').get_particle(l.get('id')).get_color()\
+                        * (-1)**(1+l.get('state'))
                 # Get the list of color flows
                 color_flow_list = \
                     matrix_element.get('color_basis').color_flow_decomposition(repr_dict,
@@ -804,19 +818,25 @@
     for iconf, configs in enumerate(s_and_t_channels):
         for vertex in configs[0] + configs[1][:-1]:
             leg = vertex.get('legs')[-1]
-            particle = particle_dict[leg.get('id')]
-            # Get mass
-            if particle.get('mass').lower() == 'zero':
-                mass = particle.get('mass')
-            else:
-                mass = "abs(%s)" % particle.get('mass')
-            # Get width
-            if particle.get('width').lower() == 'zero':
-                width = particle.get('width')
-            else:
-                width = "abs(%s)" % particle.get('width')
+            if leg.get('id') == 21 and 21 not in particle_dict:
+                # Fake propagator used in multiparticle vertices
+                mass = 'zero'
+                width = 'zero'
+                pow_part = 0
+            else:
+                particle = particle_dict[leg.get('id')]
+                # Get mass
+                if particle.get('mass').lower() == 'zero':
+                    mass = particle.get('mass')
+                else:
+                    mass = "abs(%s)" % particle.get('mass')
+                # Get width
+                if particle.get('width').lower() == 'zero':
+                    width = particle.get('width')
+                else:
+                    width = "abs(%s)" % particle.get('width')
 
-            pow_part = 1 + int(particle.is_boson())
+                pow_part = 1 + int(particle.is_boson())
 
             lines.append("pmass(%d,%d)  = %s" % \
                          (leg.get('number'), iconf + 1, mass))
@@ -951,6 +971,18 @@
                        fortran_model,
                         len(matrix_element.get_all_amplitudes()))
 
+    # Generate diagrams
+    filename = "matrix.ps"
+    plot = draw.MultiEpsDiagramDrawer(matrix_element.get('base_amplitude').\
+                                         get('diagrams'),
+                                      filename,
+                                      model=matrix_element.get('processes')[0].\
+                                         get('model'),
+                                      amplitude='')
+    logger.info("Generating Feynman diagrams for " + \
+                 matrix_element.get('processes')[0].nice_string())
+    plot.draw()
+
     linkfiles = ['check_sa.f', 'coupl.inc', 'makefile']
 
     
@@ -1263,6 +1295,29 @@
 
     return ret_list
 
+def get_amp2_lines(matrix_element):
+    """Return the amp2(i) = sum(amp for diag(i))^2 lines"""
+
+    nexternal, ninitial = matrix_element.get_nexternal_ninitial()
+    # Get minimum legs in a vertex
+    minvert = min([max(diag.get_vertex_leg_numbers()) for diag in \
+                   matrix_element.get('diagrams')])
+
+    ret_lines = []
+    for idiag, diag in enumerate(matrix_element.get('diagrams')):
+        # Ignore any diagrams with 4-particle vertices.
+        if max(diag.get_vertex_leg_numbers()) > minvert:
+            continue
+        # Now write out the expression for AMP2, meaning the sum of
+        # squared amplitudes belonging to the same diagram
+        line = "AMP2(%(num)d)=AMP2(%(num)d)+" % {"num": (idiag + 1)}
+        line += "+".join(["AMP(%(num)d)*dconjg(AMP(%(num)d))" % \
+                          {"num": a.get('number')} for a in \
+                          diag.get('amplitudes')])
+        ret_lines.append(line)
+    
+    return ret_lines
+
 def get_JAMP_lines(matrix_element):
     """Return the JAMP = sum(fermionfactor * AMP(i)) lines"""
 
@@ -1586,6 +1641,7 @@
             one_mass = particle.get('mass')
             if one_mass.lower() != 'zero':
                 masses.add(one_mass)
+                
             # find width
             one_width = particle.get('width')
             if one_width.lower() != 'zero':

=== modified file 'madgraph/iolibs/helas_call_writers.py'
--- madgraph/iolibs/helas_call_writers.py	2010-12-15 21:03:07 +0000
+++ madgraph/iolibs/helas_call_writers.py	2011-01-27 02:53:44 +0000
@@ -227,8 +227,8 @@
         key = ((3, 3, 5, 3), 'A')
         call = lambda wf: \
                "CALL UVVAXX(W(1,%d),W(1,%d),%s,zero,zero,zero,W(1,%d))" % \
-               (wf.get('mothers')[0].get('number'),
-                wf.get('mothers')[1].get('number'),
+               (FortranHelasCallWriter.sorted_mothers(wf)[0].get('number'),
+                FortranHelasCallWriter.sorted_mothers(wf)[1].get('number'),
                 wf.get('coupling'),
                 wf.get('number'))
         self.add_wavefunction(key, call)
@@ -236,8 +236,8 @@
         key = ((3, 5, 3, 1), 'A')
         call = lambda wf: \
                "CALL JVTAXX(W(1,%d),W(1,%d),%s,zero,zero,W(1,%d))" % \
-               (wf.get('mothers')[0].get('number'),
-                wf.get('mothers')[1].get('number'),
+               (FortranHelasCallWriter.sorted_mothers(wf)[0].get('number'),
+                FortranHelasCallWriter.sorted_mothers(wf)[1].get('number'),
                 wf.get('coupling'),
                 wf.get('number'))
         self.add_wavefunction(key, call)
@@ -245,9 +245,9 @@
         key = ((3, 3, 5), 'A')
         call = lambda amp: \
                "CALL VVTAXX(W(1,%d),W(1,%d),W(1,%d),%s,zero,AMP(%d))" % \
-               (amp.get('mothers')[0].get('number'),
-                amp.get('mothers')[1].get('number'),
-                amp.get('mothers')[2].get('number'),
+               (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('number'))
         self.add_amplitude(key, call)
@@ -257,56 +257,93 @@
         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))" % \
-               (wf.get('mothers')[1].get('number'),
-                wf.get('mothers')[0].get('number'),
-                wf.get('mothers')[2].get('number'),
+               (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('number'))
         self.add_wavefunction(key, call)
         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))" % \
-               (amp.get('mothers')[0].get('number'),
-                amp.get('mothers')[1].get('number'),
-                amp.get('mothers')[2].get('number'),
-                amp.get('mothers')[3].get('number'),
+               (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('number'))
         self.add_amplitude(key, call)
         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))" % \
-               (wf.get('mothers')[0].get('number'),
-                wf.get('mothers')[2].get('number'),
-                wf.get('mothers')[1].get('number'),
+               (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('number'))
         self.add_wavefunction(key, call)
         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))" % \
-               (amp.get('mothers')[2].get('number'),
-                amp.get('mothers')[0].get('number'),
-                amp.get('mothers')[1].get('number'),
-                amp.get('mothers')[3].get('number'),
+               (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('number'))
         self.add_amplitude(key, call)
         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))" % \
-               (wf.get('mothers')[2].get('number'),
-                wf.get('mothers')[1].get('number'),
-                wf.get('mothers')[0].get('number'),
+               (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('number'))
         self.add_wavefunction(key, call)
         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))" % \
-               (amp.get('mothers')[1].get('number'),
+               (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')[0].get('number'),
                 amp.get('mothers')[3].get('number'),
                 amp.get('coupling'),
                 amp.get('number'))
@@ -502,7 +539,7 @@
             if argument.needs_hermitian_conjugate():
                 call = call + 'C'
 
-            assert len(call) < 11, "Call to Helas routine %s should be maximum 6 chars" \
+            assert len(call) < 12, "Call to Helas routine %s should be maximum 6 chars" \
                       % call[5:]
 
             # Fill out with X up to 6 positions
@@ -552,9 +589,9 @@
                 # Create call for wavefunction
                 if len(argument.get('mothers')) == 2:
                     call_function = lambda wf: call % \
-                                    (wf.get('mothers')[0].\
+                                    (FortranHelasCallWriter.sorted_mothers(wf)[0].\
                                      get('number'),
-                                     wf.get('mothers')[1].\
+                                     FortranHelasCallWriter.sorted_mothers(wf)[1].\
                                      get('number'),
                                      wf.get_with_flow('coupling'),
                                      wf.get('mass'),
@@ -562,11 +599,11 @@
                                      wf.get('number'))
                 else:
                     call_function = lambda wf: call % \
-                                    (wf.get('mothers')[0].\
-                                     get('number'),
-                                     wf.get('mothers')[1].\
-                                     get('number'),
-                                     wf.get('mothers')[2].\
+                                    (FortranHelasCallWriter.sorted_mothers(wf)[0].\
+                                     get('number'),
+                                     FortranHelasCallWriter.sorted_mothers(wf)[1].\
+                                     get('number'),
+                                     FortranHelasCallWriter.sorted_mothers(wf)[2].\
                                      get('number'),
                                      wf.get_with_flow('coupling'),
                                      wf.get('mass'),
@@ -576,24 +613,24 @@
                 # Create call for amplitude
                 if len(argument.get('mothers')) == 3:
                     call_function = lambda amp: call % \
-                                    (amp.get('mothers')[0].\
-                                     get('number'),
-                                     amp.get('mothers')[1].\
-                                     get('number'),
-                                     amp.get('mothers')[2].\
+                                    (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('number'))
                 else:
                     call_function = lambda amp: call % \
-                                    (amp.get('mothers')[0].\
-                                     get('number'),
-                                     amp.get('mothers')[1].\
-                                     get('number'),
-                                     amp.get('mothers')[2].\
-                                     get('number'),
-                                     amp.get('mothers')[3].\
+                                    (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('number'))
@@ -625,6 +662,95 @@
                           FortranHelasCallWriter.sort_amp[l2] - \
                           FortranHelasCallWriter.sort_amp[l1]))
 
+    @staticmethod
+    def sorted_mothers(arg):
+        """Gives a list of mother wavefunctions sorted according to
+        1. The order of the particles in the interaction
+        2. Cyclic reordering of particles in same spin group
+        3. Fermions ordered IOIOIO... according to the pairs in
+           the interaction."""
+
+        assert isinstance(arg, (helas_objects.HelasWavefunction, helas_objects.HelasAmplitude)), \
+            "%s is not a valid HelasWavefunction or HelasAmplitude" % repr(arg)
+
+        if not arg.get('interaction_id'):
+            return arg.get('mothers')
+        my_pdg_code = 0
+        my_spin = 0
+        if isinstance(arg, helas_objects.HelasWavefunction):
+            my_pdg_code = arg.get_anti_pdg_code()
+            my_spin = arg.get_spin_state_number()
+
+        sorted_mothers, my_index = arg.get('mothers').sort_by_pdg_codes(\
+            arg.get('pdg_codes'), my_pdg_code)
+
+        # If fermion, partner is the corresponding fermion flow partner
+        partner = None
+        if isinstance(arg, helas_objects.HelasWavefunction) and arg.is_fermion():
+            # Fermion case, just pick out the fermion flow partner
+            if my_index % 2 == 0:
+                # partner is after arg
+                partner_index = my_index
+            else:
+                # partner is before arg
+                partner_index = my_index - 1
+            partner = sorted_mothers.pop(partner_index)
+            # If partner is incoming, move to before arg
+            if partner.get_spin_state_number() > 0:
+                my_index = partner_index
+            else:
+                my_index = partner_index + 1
+
+        # Reorder fermions pairwise according to incoming/outgoing
+        for i in range(0, len(sorted_mothers), 2):
+            if sorted_mothers[i].is_fermion():
+                # This is a fermion, order between this fermion and its brother
+                if sorted_mothers[i].get_spin_state_number() > 0 and \
+                   sorted_mothers[i + 1].get_spin_state_number() < 0:
+                    # Switch places between outgoing and incoming
+                    sorted_mothers = sorted_mothers[:i] + \
+                                      [sorted_mothers[i+1], sorted_mothers[i]] + \
+                                      sorted_mothers[i+2:]
+                elif sorted_mothers[i].get_spin_state_number() < 0 and \
+                   sorted_mothers[i + 1].get_spin_state_number() > 0:
+                    # This is the right order
+                    pass
+            else:
+                # No more fermions in sorted_mothers
+                break
+            
+        # Put back partner into sorted_mothers
+        if partner:
+            sorted_mothers.insert(partner_index, partner)
+
+        same_spin_mothers = []
+        if isinstance(arg, helas_objects.HelasWavefunction):
+            # Pick out mothers with same spin, for cyclic reordering
+            same_spin_index = -1
+            i=0
+            while i < len(sorted_mothers):
+                if abs(sorted_mothers[i].get_spin_state_number()) == \
+                       abs(my_spin):
+                    if same_spin_index < 0:
+                        # Remember starting index for same spin states
+                        same_spin_index = i
+                    same_spin_mothers.append(sorted_mothers.pop(i))
+                else:
+                    i += 1
+
+        # Make cyclic reordering of mothers with same spin as this wf
+        if same_spin_mothers:
+            same_spin_mothers = same_spin_mothers[my_index - same_spin_index:] \
+                                + same_spin_mothers[:my_index - same_spin_index]
+
+            # Insert same_spin_mothers in sorted_mothers
+            sorted_mothers = sorted_mothers[:same_spin_index] + \
+                              same_spin_mothers + sorted_mothers[same_spin_index:]
+
+        # Next sort according to spin_state_number
+        return helas_objects.HelasWavefunctionList(sorted_mothers)
+
+
 #===============================================================================
 # UFOHelasCallWriter
 #===============================================================================
@@ -754,7 +880,8 @@
             # Check if we need to append a charge conjugation flag
             c_flag = '' 
             if argument.needs_hermitian_conjugate():
-                c_flag = 'C1' # MG5 not configure for 4F vertex
+                c_flag = "".join(['C%d' % i for i in \
+                                  argument.get_conjugate_index()])
 
             call = 'CALL %s%s_%s' % (argument.get('lorentz'), c_flag, outgoing) 
 
@@ -888,7 +1015,8 @@
             # Check if we need to append a charge conjugation flag
             c_flag = '' 
             if argument.needs_hermitian_conjugate():
-                c_flag = 'C1' # MG5 not configure for 4F vertex
+                c_flag = "".join(['C%d' % i for i in \
+                                  argument.get_conjugate_index()])
 
             call = '%s%s_%s' % (argument.get('lorentz'), c_flag, outgoing)
 
@@ -1055,7 +1183,8 @@
             # Check if we need to append a charge conjugation flag
             c_flag = '' 
             if argument.needs_hermitian_conjugate():
-                c_flag = 'C1' # MG5 not configure for 4F vertex
+                c_flag = "".join(['C%d' % i for i in \
+                                  argument.get_conjugate_index()])
 
             if isinstance(argument, helas_objects.HelasWavefunction):
                 call = 'w[%d] = '

=== modified file 'madgraph/iolibs/import_v4.py'
--- madgraph/iolibs/import_v4.py	2010-12-18 22:39:25 +0000
+++ madgraph/iolibs/import_v4.py	2011-01-27 02:53:44 +0000
@@ -80,7 +80,12 @@
 
     # treat simple case (model_path is a valid path/ mgme_dir doesn't exist)
     if os.path.isdir(model_path):
-        return model_path
+        pass
+    elif mgme_dir and os.path.isdir(os.path.join(mgme_dir, 'models',
+                                                 model_path + "_v4")):
+        model_path = os.path.join(mgme_dir, 'models', model_path + "_v4")
+    elif mgme_dir and os.path.isdir(os.path.join(mgme_dir, 'Models', model_path)):
+        model_path = os.path.join(mgme_dir, 'Models', model_path)
     elif not mgme_dir:
         error_text = "Path %s is not a valid pathname\n" % model_path
         error_text += "and no MG_ME installation detected in order to search in Models"

=== modified file 'madgraph/iolibs/template_files/matrix_madevent_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_v4.inc	2010-10-18 21:37:36 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_v4.inc	2011-01-27 02:53:44 +0000
@@ -198,9 +198,7 @@
         ENDDO
         MATRIX =MATRIX+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
     ENDDO
-    Do I = 1, NGRAPHS
-        amp2(i)=amp2(i)+amp(i)*dconjg(amp(i))
-    Enddo
+%(amp2_lines)s
     Do I = 1, NCOLOR
         Jamp2(i)=Jamp2(i)+Jamp(i)*dconjg(Jamp(i))
     Enddo

=== modified file 'madgraph/various/process_checks.py'
--- madgraph/various/process_checks.py	2010-12-15 21:03:07 +0000
+++ madgraph/various/process_checks.py	2011-01-27 02:53:44 +0000
@@ -41,12 +41,12 @@
 
 import madgraph.various.rambo as rambo
 
-from madgraph import MG5DIR
+from madgraph import MG5DIR, MadGraph5Error
 
 import models.model_reader as model_reader
 import aloha.template_files.wavefunctions as wavefunctions
 from aloha.template_files.wavefunctions import \
-     ixxxxx, oxxxxx, vxxxxx, sxxxxx
+     ixxxxx, oxxxxx, vxxxxx, sxxxxx, txxxxx
 
 #===============================================================================
 # Logger for process_checks
@@ -110,7 +110,7 @@
         return p, w_rambo
 
     if nincoming != 2:
-        raise RAMBOError('Need 1 or 2 incoming particles')
+        raise rambo.RAMBOError('Need 1 or 2 incoming particles')
 
     if nfinal == 1:
         energy = masses[0]
@@ -1070,4 +1070,3 @@
         return res_str        
     else: 
         return fail_proc
-    return res_str        

=== modified file 'models/template_files/fortran/formats.inc'
--- models/template_files/fortran/formats.inc	2010-06-02 23:59:09 +0000
+++ models/template_files/fortran/formats.inc	2011-01-27 02:53:44 +0000
@@ -11,7 +11,7 @@
 c Simple real
  1    format(1x,a15,f11.5)
 c Simple Complex
- 2    format(1x,a15,f11.5,1x,f11.5)
+ 2    format(1x,a15,e13.5,1x,e13.5)
 c Real with mass dimension
  3    format(1x,a15,f11.5,' GeV')
 c Chiral couplings

=== added directory 'models/uutt_sch_4fermion'
=== added file 'models/uutt_sch_4fermion/__init__.py'
--- models/uutt_sch_4fermion/__init__.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/__init__.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,15 @@
+
+import particles
+import couplings
+import lorentz
+import parameters
+import vertices
+import write_param_card
+
+
+all_particles = particles.all_particles
+all_vertices = vertices.all_vertices
+all_couplings = couplings.all_couplings
+all_lorentz = lorentz.all_lorentz
+all_parameters = parameters.all_parameters
+all_functions = function_library.all_functions

=== added file 'models/uutt_sch_4fermion/couplings.py'
--- models/uutt_sch_4fermion/couplings.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/couplings.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,119 @@
+# This file was automatically created by FeynRules $Revision: 302 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Tue 31 Aug 2010 16:54:46
+
+
+from object_library import all_couplings, Coupling
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+
+
+GC_1 = Coupling(name = 'GC_1',
+                value = '-(ee*complex(0,1))/3.',
+                order = {'QED':1})
+
+GC_2 = Coupling(name = 'GC_2',
+                value = '(2*ee*complex(0,1))/3.',
+                order = {'QED':1})
+
+GC_3 = Coupling(name = 'GC_3',
+                value = '-(ee*complex(0,1))',
+                order = {'QED':1})
+
+GC_4 = Coupling(name = 'GC_4',
+                value = '-G',
+                order = {'QCD':1})
+
+GC_5 = Coupling(name = 'GC_5',
+                value = 'complex(0,1)*G',
+                order = {'QCD':1})
+
+GC_6 = Coupling(name = 'GC_6',
+                value = 'complex(0,1)*G**2',
+                order = {'QCD':2})
+
+GC_7 = Coupling(name = 'GC_7',
+                value = 'cw*complex(0,1)*gw',
+                order = {'QED':1})
+
+GC_8 = Coupling(name = 'GC_8',
+                value = '-(complex(0,1)*gw**2)',
+                order = {'QED':2})
+
+GC_9 = Coupling(name = 'GC_9',
+                value = 'cw**2*complex(0,1)*gw**2',
+                order = {'QED':2})
+
+GC_10 = Coupling(name = 'GC_10',
+                 value = '(ee**2*complex(0,1))/(2.*sw**2)',
+                 order = {'QED':2})
+
+GC_11 = Coupling(name = 'GC_11',
+                 value = '(ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_12 = Coupling(name = 'GC_12',
+                 value = '-(cw*ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_13 = Coupling(name = 'GC_13',
+                 value = '(cw*ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_14 = Coupling(name = 'GC_14',
+                 value = '-(ee*complex(0,1)*sw)/(6.*cw)',
+                 order = {'QED':1})
+
+GC_15 = Coupling(name = 'GC_15',
+                 value = '(ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_16 = Coupling(name = 'GC_16',
+                 value = 'complex(0,1)*gw*sw',
+                 order = {'QED':1})
+
+GC_17 = Coupling(name = 'GC_17',
+                 value = '-2*cw*complex(0,1)*gw**2*sw',
+                 order = {'QED':2})
+
+GC_18 = Coupling(name = 'GC_18',
+                 value = 'complex(0,1)*gw**2*sw**2',
+                 order = {'QED':2})
+
+GC_19 = Coupling(name = 'GC_19',
+                 value = '(cw*ee*complex(0,1))/(2.*sw) + (ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_20 = Coupling(name = 'GC_20',
+                 value = 'ee**2*complex(0,1) + (cw**2*ee**2*complex(0,1))/(2.*sw**2) + (ee**2*complex(0,1)*sw**2)/(2.*cw**2)',
+                 order = {'QED':2})
+
+GC_21 = Coupling(name = 'GC_21',
+                 value = '-6*complex(0,1)*lam*v',
+                 order = {'QED':1})
+
+GC_22 = Coupling(name = 'GC_22',
+                 value = '(ee**2*complex(0,1)*v)/(2.*sw**2)',
+                 order = {'QED':1})
+
+GC_23 = Coupling(name = 'GC_23',
+                 value = 'ee**2*complex(0,1)*v + (cw**2*ee**2*complex(0,1)*v)/(2.*sw**2) + (ee**2*complex(0,1)*sw**2*v)/(2.*cw**2)',
+                 order = {'QED':1})
+
+GC_24 = Coupling(name = 'GC_24',
+                 value = '-((complex(0,1)*yb)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_25 = Coupling(name = 'GC_25',
+                 value = '-((complex(0,1)*yt)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_26 = Coupling(name = 'GC_26',
+                 value = '-((complex(0,1)*ytau)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+# 4-fermion coupling is 2*g^2*prop for t-channel (due to FS symmetry)
+GC_27 = Coupling(name = 'GC_27',
+                 value = '-uusix*ttsix/Msix**2',
+                 order = {'NP':2})

=== added file 'models/uutt_sch_4fermion/function_library.py'
--- models/uutt_sch_4fermion/function_library.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/function_library.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,54 @@
+# This file is part of the UFO.
+#
+# This file contains definitions for functions that
+# are extensions of the cmath library, and correspond
+# either to functions that are in cmath, but inconvenient
+# to access from there (e.g. z.conjugate()),
+# or functions that are simply not defined.
+#
+#
+
+__date__ = "22 July 2010"
+__author__ = "[e-mail]"
+
+import cmath
+from object_library import all_functions, Function
+
+#
+# shortcuts for functions from cmath
+#
+
+complexconjugate = Function(name = 'complexconjugate',
+                            arguments = ('z',),
+                            expression = 'z.conjugate()')
+
+
+re = Function(name = 're',
+              arguments = ('z',),
+              expression = 'z.real')
+
+im = Function(name = 'im',
+              arguments = ('z',),
+              expression = 'z.imag')
+
+# New functions (trigonometric)
+
+sec = Function(name = 'sec',
+             arguments = ('z',),
+             expression = '1./cmath.cos(z)')
+
+asec = Function(name = 'asec',
+             arguments = ('z',),
+             expression = 'cmath.acos(1./z)')
+
+csc = Function(name = 'csc',
+             arguments = ('z',),
+             expression = '1./cmath.sin(z)')
+
+acsc = Function(name = 'acsc',
+             arguments = ('z',),
+             expression = 'cmath.asin(1./z)')
+
+
+
+

=== added file 'models/uutt_sch_4fermion/lorentz.py'
--- models/uutt_sch_4fermion/lorentz.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/lorentz.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,75 @@
+# This file was automatically created by FeynRules $Revision: 302 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Tue 31 Aug 2010 16:54:46
+
+
+from object_library import all_lorentz, Lorentz
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+
+
+SSS1 = Lorentz(name = 'SSS1',
+               spins = [ 1, 1, 1 ],
+               structure = '1')
+
+FFS1 = Lorentz(name = 'FFS1',
+               spins = [ 2, 2, 1 ],
+               structure = 'Identity(1,2)')
+
+FFV1 = Lorentz(name = 'FFV1',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,1)')
+
+FFV2 = Lorentz(name = 'FFV2',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,\'s1\')*ProjM(\'s1\',1)')
+
+FFV3 = Lorentz(name = 'FFV3',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,\'s1\')*ProjM(\'s1\',1) - 2*Gamma(3,2,\'s1\')*ProjP(\'s1\',1)')
+
+FFV4 = Lorentz(name = 'FFV4',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,\'s1\')*ProjM(\'s1\',1) + 2*Gamma(3,2,\'s1\')*ProjP(\'s1\',1)')
+
+FFV5 = Lorentz(name = 'FFV5',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,\'s1\')*ProjM(\'s1\',1) + 4*Gamma(3,2,\'s1\')*ProjP(\'s1\',1)')
+
+VVS1 = Lorentz(name = 'VVS1',
+               spins = [ 3, 3, 1 ],
+               structure = 'Metric(1,2)')
+
+VVV1 = Lorentz(name = 'VVV1',
+               spins = [ 3, 3, 3 ],
+               structure = 'P(3,1)*Metric(1,2) - P(3,2)*Metric(1,2) - P(2,1)*Metric(1,3) + P(2,3)*Metric(1,3) + P(1,2)*Metric(2,3) - P(1,3)*Metric(2,3)')
+
+VVSS1 = Lorentz(name = 'VVSS1',
+                spins = [ 3, 3, 1, 1 ],
+                structure = 'Metric(1,2)')
+
+VVVV1 = Lorentz(name = 'VVVV1',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - Metric(1,3)*Metric(2,4)')
+
+VVVV2 = Lorentz(name = 'VVVV2',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) + Metric(1,3)*Metric(2,4) - 2*Metric(1,2)*Metric(3,4)')
+
+VVVV3 = Lorentz(name = 'VVVV3',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - Metric(1,2)*Metric(3,4)')
+
+VVVV4 = Lorentz(name = 'VVVV4',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,3)*Metric(2,4) - Metric(1,2)*Metric(3,4)')
+
+VVVV5 = Lorentz(name = 'VVVV5',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - (Metric(1,3)*Metric(2,4))/2. - (Metric(1,2)*Metric(3,4))/2.')
+
+FFFF1 = Lorentz(name = 'FFFF1',
+                spins = [ 2, 2, 2, 2 ],
+                structure = 'Identity(1,2)*Identity(4,3)')
+

=== added file 'models/uutt_sch_4fermion/object_library.py'
--- models/uutt_sch_4fermion/object_library.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/object_library.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,191 @@
+##
+##
+## Feynrules Header
+##
+##
+##
+##
+##
+
+import cmath
+
+
+class FRBaseClass(object):
+    """The class from which all FeynRules classes are derived."""
+
+    require_args = []
+
+    def __init__(self, *args, **options):
+        assert(len(self.require_args) == len (args))
+    
+        for i, name in enumerate(self.require_args):
+            setattr(self, name, args[i])
+    
+        for (option, value) in options.items():
+            setattr(self, option, value)
+
+    def get(self, name):
+        return getattr(self, name)
+    
+    def set(self, name, value):
+        setattr(self, name, value)
+        
+    def get_all(self):
+        """Return a dictionary containing all the information of the object"""
+        return self.__dict__
+
+    def __str__(self):
+        return self.name
+
+
+    def __repr__(self):
+        replacements = [
+            ('+','__plus__'),
+            ('-','__minus__'),
+            ('@','__at__'),
+            ('!','__exclam__'),
+            ('?','__quest__'),
+            ('*','__star__'),
+            ('~','__tilde__')
+            ]
+        text = self.name
+        for orig,sub in replacements:
+            text = text.replace(orig,sub)
+        return text
+
+
+
+all_particles = []
+
+    
+
+class Particle(FRBaseClass):
+    """A standard Particle"""
+
+    require_args=['pdg_code', 'name', 'antiname', 'spin', 'color', 'mass', 'width', 'texname', 'antitexname', 'line', 'charge']
+
+    require_args_all = ['pdg_code', 'name', 'antiname', 'spin', 'color', 'mass', 'width', 'texname', 'antitexname', 'line', 'charge', 'propagating', 'goldstoneboson']
+
+    def __init__(self, pdg_code, name, antiname, spin, color, mass, width, texname,
+                 antitexname, line, charge , propagating=True, goldstoneboson=False, **options):
+
+        args= (pdg_code, name, antiname, spin, color, mass, width, texname,
+                 antitexname, line, float(charge))
+
+        FRBaseClass.__init__(self, *args,  **options)
+
+        global all_particles
+        all_particles.append(self)
+
+        self.propagating = propagating
+        self.goldstoneboson= goldstoneboson
+
+        self.selfconjugate = (name == antiname)
+
+    def anti(self):
+        if self.selfconjugate:
+            raise Exception('%s has no anti particle.' % self.name) 
+        outdic = {}
+        for k,v in self.__dict__.iteritems():
+            if k not in self.require_args_all:                
+                outdic[k] = -v
+        if self.color in [1,8]:
+            newcolor = self.color
+        else:
+            newcolor = -self.color
+                
+        return Particle(-self.pdg_code, self.antiname, self.name, self.spin, newcolor, self.mass, self.width,
+                        self.antitexname, self.texname, self.line, -self.charge, self.propagating, self.goldstoneboson, **outdic)
+
+
+
+all_parameters = []
+
+class Parameter(FRBaseClass):
+
+    require_args=['name', 'nature', 'type', 'value', 'texname']
+
+    def __init__(self, name, nature, type, value, texname, lhablock=None, lhacode=None):
+
+        args = (name,nature,type,value,texname)
+
+        FRBaseClass.__init__(self, *args)
+
+        args=(name,nature,type,value,texname)
+
+        global all_parameters
+        all_parameters.append(self)
+
+        if (lhablock is None or lhacode is None)  and nature == 'external':
+            raise Exception('Need LHA information for external parameter "%s".' % name)
+        self.lhablock = lhablock
+        self.lhacode = lhacode
+
+all_vertices = []
+
+class Vertex(FRBaseClass):
+
+    require_args=['name', 'particles', 'color', 'lorentz', 'couplings']
+
+    def __init__(self, name, particles, color, lorentz, couplings, **opt):
+ 
+        args = (name, particles, color, lorentz, couplings)
+
+        FRBaseClass.__init__(self, *args, **opt)
+
+        args=(particles,color,lorentz,couplings)
+
+        global all_vertices
+        all_vertices.append(self)
+
+all_couplings = []
+
+class Coupling(FRBaseClass):
+
+    require_args=['name', 'value', 'order']
+
+    def __init__(self, name, value, order, **opt):
+
+        args =(name, value, order)	
+        FRBaseClass.__init__(self, *args, **opt)
+        global all_couplings
+        all_couplings.append(self)
+  
+
+
+all_lorentz = []
+
+class Lorentz(FRBaseClass):
+
+    require_args=['name','spins','structure']
+    
+    def __init__(self, name, spins, structure, **opt):
+        args = (name, spins, structure)
+        FRBaseClass.__init__(self, *args, **opt)
+
+        global all_lorentz
+        all_lorentz.append(self)
+
+
+all_functions = []
+
+class Function(object):
+
+    def __init__(self, name, arguments, expression):
+
+        global all_functions
+        all_functions.append(self)
+
+        self.name = name
+        self.arguments = arguments
+        self.expr = expression
+    
+    def __call__(self, *opt):
+
+        for i, arg in enumerate(self.arguments):
+            exec('%s = %s' % (arg, opt[i] ))
+
+        return eval(self.expr)
+
+
+        

=== added file 'models/uutt_sch_4fermion/parameters.py'
--- models/uutt_sch_4fermion/parameters.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/parameters.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,245 @@
+# This file was automatically created by FeynRules $Revision: 302 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Tue 31 Aug 2010 16:54:46
+
+
+
+from object_library import all_parameters, Parameter
+
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+aEWM1 = Parameter(name = 'aEWM1',
+                  nature = 'external',
+                  type = 'real',
+                  value = 132.50698,
+                  texname = '\\text{aEWM1}',
+                  lhablock = 'SMINPUTS',
+                  lhacode = [ 1 ])
+
+Gf = Parameter(name = 'Gf',
+               nature = 'external',
+               type = 'real',
+               value = 0.0000116639,
+               texname = 'G_f',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 2 ])
+
+aS = Parameter(name = 'aS',
+               nature = 'external',
+               type = 'real',
+               value = 0.118,
+               texname = '\\text{aS}',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 3 ])
+
+ymb = Parameter(name = 'ymb',
+                nature = 'external',
+                type = 'real',
+                value = 4.2,
+                texname = '\\text{ymb}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 5 ])
+
+ymt = Parameter(name = 'ymt',
+                nature = 'external',
+                type = 'real',
+                value = 164.5,
+                texname = '\\text{ymt}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 6 ])
+
+ymtau = Parameter(name = 'ymtau',
+                  nature = 'external',
+                  type = 'real',
+                  value = 1.777,
+                  texname = '\\text{ymtau}',
+                  lhablock = 'YUKAWA',
+                  lhacode = [ 15 ])
+
+MTA = Parameter(name = 'MTA',
+                nature = 'external',
+                type = 'real',
+                value = 1.777,
+                texname = '\\text{MTA}',
+                lhablock = 'MASS',
+                lhacode = [ 15 ])
+
+MT = Parameter(name = 'MT',
+               nature = 'external',
+               type = 'real',
+               value = 174.3,
+               texname = '\\text{MT}',
+               lhablock = 'MASS',
+               lhacode = [ 6 ])
+
+MB = Parameter(name = 'MB',
+               nature = 'external',
+               type = 'real',
+               value = 4.7,
+               texname = '\\text{MB}',
+               lhablock = 'MASS',
+               lhacode = [ 5 ])
+
+MZ = Parameter(name = 'MZ',
+               nature = 'external',
+               type = 'real',
+               value = 91.188,
+               texname = '\\text{MZ}',
+               lhablock = 'MASS',
+               lhacode = [ 23 ])
+
+MH = Parameter(name = 'MH',
+               nature = 'external',
+               type = 'real',
+               value = 120.,
+               texname = '\\text{MH}',
+               lhablock = 'MASS',
+               lhacode = [ 25 ])
+
+WT = Parameter(name = 'WT',
+               nature = 'external',
+               type = 'real',
+               value = 1.50833649,
+               texname = '\\text{WT}',
+               lhablock = 'DECAY',
+               lhacode = [ 6 ])
+
+WZ = Parameter(name = 'WZ',
+               nature = 'external',
+               type = 'real',
+               value = 2.44140351,
+               texname = '\\text{WZ}',
+               lhablock = 'DECAY',
+               lhacode = [ 23 ])
+
+WW = Parameter(name = 'WW',
+               nature = 'external',
+               type = 'real',
+               value = 2.04759951,
+               texname = '\\text{WW}',
+               lhablock = 'DECAY',
+               lhacode = [ 24 ])
+
+WH = Parameter(name = 'WH',
+               nature = 'external',
+               type = 'real',
+               value = 0.00575308848,
+               texname = '\\text{WH}',
+               lhablock = 'DECAY',
+               lhacode = [ 25 ])
+
+aEW = Parameter(name = 'aEW',
+                nature = 'internal',
+                type = 'real',
+                value = '1/aEWM1',
+                texname = '\\text{aEW}')
+
+G = Parameter(name = 'G',
+              nature = 'internal',
+              type = 'real',
+              value = '2*cmath.sqrt(aS)*cmath.sqrt(cmath.pi)',
+              texname = 'G')
+
+MW = Parameter(name = 'MW',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(MZ**2/2. + cmath.sqrt(MZ**4/4. - (aEW*cmath.pi*MZ**2)/(Gf*cmath.sqrt(2))))',
+               texname = 'M_W')
+
+ee = Parameter(name = 'ee',
+               nature = 'internal',
+               type = 'real',
+               value = '2*cmath.sqrt(aEW)*cmath.sqrt(cmath.pi)',
+               texname = 'e')
+
+sw2 = Parameter(name = 'sw2',
+                nature = 'internal',
+                type = 'real',
+                value = '1 - MW**2/MZ**2',
+                texname = '\\text{sw2}')
+
+cw = Parameter(name = 'cw',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(1 - sw2)',
+               texname = 'c_w')
+
+sw = Parameter(name = 'sw',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(sw2)',
+               texname = 's_w')
+
+g1 = Parameter(name = 'g1',
+               nature = 'internal',
+               type = 'real',
+               value = 'ee/cw',
+               texname = 'g_1')
+
+gw = Parameter(name = 'gw',
+               nature = 'internal',
+               type = 'real',
+               value = 'ee/sw',
+               texname = 'g_w')
+
+v = Parameter(name = 'v',
+              nature = 'internal',
+              type = 'real',
+              value = '(2*MW*sw)/ee',
+              texname = 'v')
+
+lam = Parameter(name = 'lam',
+                nature = 'internal',
+                type = 'real',
+                value = 'MH**2/(2.*v**2)',
+                texname = '\\text{lam}')
+
+yb = Parameter(name = 'yb',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymb*cmath.sqrt(2))/v',
+               texname = '\\text{yb}')
+
+yt = Parameter(name = 'yt',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymt*cmath.sqrt(2))/v',
+               texname = '\\text{yt}')
+
+ytau = Parameter(name = 'ytau',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '(ymtau*cmath.sqrt(2))/v',
+                 texname = '\\text{ytau}')
+
+muH = Parameter(name = 'muH',
+                nature = 'internal',
+                type = 'real',
+                value = 'cmath.sqrt(lam*v**2)',
+                texname = '\\mu ')
+
+ttsix = Parameter(name = 'ttsix',
+                  nature = 'external',
+                  type = 'real',
+                  value = 1,
+                  texname = '\\text{ttsix}',
+                  lhablock = 'BSMINPUTS',
+                  lhacode = [ 1 ])
+
+uusix = Parameter(name = 'uusix',
+                  nature = 'external',
+                  type = 'real',
+                  value = 0.1,
+                  texname = '\\text{ttsix}',
+                  lhablock = 'BSMINPUTS',
+                  lhacode = [ 2 ])
+
+Msix = Parameter(name = 'Msix',
+                nature = 'external',
+                type = 'real',
+                value = 1.e5,
+                texname = '\\text{Mut}',
+                lhablock = 'BSMINPUTS',
+                lhacode = [ 3 ])
+

=== added file 'models/uutt_sch_4fermion/particles.py'
--- models/uutt_sch_4fermion/particles.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/particles.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,384 @@
+# This file was automatically created by FeynRules $Revision: 302 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Tue 31 Aug 2010 16:54:46
+
+
+from __future__ import division
+from object_library import all_particles, Particle
+
+ve = Particle(pdg_code = 12,
+              name = 've',
+              antiname = 've~',
+              spin = 2,
+              color = 1,
+              mass = 'ZERO',
+              width = 'ZERO',
+              texname = 've',
+              antitexname = 've',
+              line = 'straight',
+              charge = 0,
+              LeptonNumber = 1,
+              GhostNumber = 0)
+
+ve__tilde__ = ve.anti()
+
+vm = Particle(pdg_code = 14,
+              name = 'vm',
+              antiname = 'vm~',
+              spin = 2,
+              color = 1,
+              mass = 'ZERO',
+              width = 'ZERO',
+              texname = 'vm',
+              antitexname = 'vm',
+              line = 'straight',
+              charge = 0,
+              LeptonNumber = 1,
+              GhostNumber = 0)
+
+vm__tilde__ = vm.anti()
+
+vt = Particle(pdg_code = 16,
+              name = 'vt',
+              antiname = 'vt~',
+              spin = 2,
+              color = 1,
+              mass = 'ZERO',
+              width = 'ZERO',
+              texname = 'vt',
+              antitexname = 'vt',
+              line = 'straight',
+              charge = 0,
+              LeptonNumber = 1,
+              GhostNumber = 0)
+
+vt__tilde__ = vt.anti()
+
+e__minus__ = Particle(pdg_code = 11,
+                      name = 'e-',
+                      antiname = 'e+',
+                      spin = 2,
+                      color = 1,
+                      mass = 'ZERO',
+                      width = 'ZERO',
+                      texname = 'e-',
+                      antitexname = 'e-',
+                      line = 'straight',
+                      charge = -1,
+                      LeptonNumber = 1,
+                      GhostNumber = 0)
+
+e__plus__ = e__minus__.anti()
+
+m__minus__ = Particle(pdg_code = 13,
+                      name = 'm-',
+                      antiname = 'm+',
+                      spin = 2,
+                      color = 1,
+                      mass = 'ZERO',
+                      width = 'ZERO',
+                      texname = 'm-',
+                      antitexname = 'm-',
+                      line = 'straight',
+                      charge = -1,
+                      LeptonNumber = 1,
+                      GhostNumber = 0)
+
+m__plus__ = m__minus__.anti()
+
+tt__minus__ = Particle(pdg_code = 15,
+                       name = 'tt-',
+                       antiname = 'tt+',
+                       spin = 2,
+                       color = 1,
+                       mass = 'MTA',
+                       width = 'ZERO',
+                       texname = 'tt-',
+                       antitexname = 'tt-',
+                       line = 'straight',
+                       charge = -1,
+                       LeptonNumber = 1,
+                       GhostNumber = 0)
+
+tt__plus__ = tt__minus__.anti()
+
+u = Particle(pdg_code = 2,
+             name = 'u',
+             antiname = 'u~',
+             spin = 2,
+             color = 3,
+             mass = 'ZERO',
+             width = 'ZERO',
+             texname = 'u',
+             antitexname = 'u',
+             line = 'straight',
+             charge = 2/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+u__tilde__ = u.anti()
+
+c = Particle(pdg_code = 4,
+             name = 'c',
+             antiname = 'c~',
+             spin = 2,
+             color = 3,
+             mass = 'ZERO',
+             width = 'ZERO',
+             texname = 'c',
+             antitexname = 'c',
+             line = 'straight',
+             charge = 2/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+c__tilde__ = c.anti()
+
+t = Particle(pdg_code = 6,
+             name = 't',
+             antiname = 't~',
+             spin = 2,
+             color = 3,
+             mass = 'MT',
+             width = 'WT',
+             texname = 't',
+             antitexname = 't',
+             line = 'straight',
+             charge = 2/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+t__tilde__ = t.anti()
+
+d = Particle(pdg_code = 1,
+             name = 'd',
+             antiname = 'd~',
+             spin = 2,
+             color = 3,
+             mass = 'ZERO',
+             width = 'ZERO',
+             texname = 'd',
+             antitexname = 'd',
+             line = 'straight',
+             charge = -1/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+d__tilde__ = d.anti()
+
+s = Particle(pdg_code = 3,
+             name = 's',
+             antiname = 's~',
+             spin = 2,
+             color = 3,
+             mass = 'ZERO',
+             width = 'ZERO',
+             texname = 's',
+             antitexname = 's',
+             line = 'straight',
+             charge = -1/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+s__tilde__ = s.anti()
+
+b = Particle(pdg_code = 5,
+             name = 'b',
+             antiname = 'b~',
+             spin = 2,
+             color = 3,
+             mass = 'MB',
+             width = 'ZERO',
+             texname = 'b',
+             antitexname = 'b',
+             line = 'straight',
+             charge = -1/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+b__tilde__ = b.anti()
+
+ghA = Particle(pdg_code = 9000001,
+               name = 'ghA',
+               antiname = 'ghA~',
+               spin = -1,
+               color = 1,
+               mass = 'ZERO',
+               width = 'WghA',
+               texname = 'ghA',
+               antitexname = 'ghA',
+               line = 'dotted',
+               charge = 0,
+               LeptonNumber = 0,
+               GhostNumber = 1)
+
+ghA__tilde__ = ghA.anti()
+
+ghZ = Particle(pdg_code = 9000002,
+               name = 'ghZ',
+               antiname = 'ghZ~',
+               spin = -1,
+               color = 1,
+               mass = 'MZ',
+               width = 'WghZ',
+               texname = 'ghZ',
+               antitexname = 'ghZ',
+               line = 'dotted',
+               charge = 0,
+               LeptonNumber = 0,
+               GhostNumber = 1)
+
+ghZ__tilde__ = ghZ.anti()
+
+ghWp = Particle(pdg_code = 9000003,
+                name = 'ghWp',
+                antiname = 'ghWp~',
+                spin = -1,
+                color = 1,
+                mass = 'MW',
+                width = 'WghWp',
+                texname = 'ghWp',
+                antitexname = 'ghWp',
+                line = 'dotted',
+                charge = 1,
+                LeptonNumber = 0,
+                GhostNumber = 1)
+
+ghWp__tilde__ = ghWp.anti()
+
+ghWm = Particle(pdg_code = 9000004,
+                name = 'ghWm',
+                antiname = 'ghWm~',
+                spin = -1,
+                color = 1,
+                mass = 'MW',
+                width = 'WghWm',
+                texname = 'ghWm',
+                antitexname = 'ghWm',
+                line = 'dotted',
+                charge = -1,
+                LeptonNumber = 0,
+                GhostNumber = 1)
+
+ghWm__tilde__ = ghWm.anti()
+
+ghG = Particle(pdg_code = 9000005,
+               name = 'ghG',
+               antiname = 'ghG~',
+               spin = -1,
+               color = 8,
+               mass = 'ZERO',
+               width = 'WghG',
+               texname = 'ghG',
+               antitexname = 'ghG',
+               line = 'dotted',
+               charge = 0,
+               LeptonNumber = 0,
+               GhostNumber = 1)
+
+ghG__tilde__ = ghG.anti()
+
+A = Particle(pdg_code = 22,
+             name = 'A',
+             antiname = 'A',
+             spin = 3,
+             color = 1,
+             mass = 'ZERO',
+             width = 'ZERO',
+             texname = 'A',
+             antitexname = 'A',
+             line = 'wavy',
+             charge = 0,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+Z = Particle(pdg_code = 23,
+             name = 'Z',
+             antiname = 'Z',
+             spin = 3,
+             color = 1,
+             mass = 'MZ',
+             width = 'WZ',
+             texname = 'Z',
+             antitexname = 'Z',
+             line = 'wavy',
+             charge = 0,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+W__plus__ = Particle(pdg_code = 24,
+                     name = 'W+',
+                     antiname = 'W-',
+                     spin = 3,
+                     color = 1,
+                     mass = 'MW',
+                     width = 'WW',
+                     texname = 'W+',
+                     antitexname = 'W+',
+                     line = 'wavy',
+                     charge = 1,
+                     LeptonNumber = 0,
+                     GhostNumber = 0)
+
+W__minus__ = W__plus__.anti()
+
+G = Particle(pdg_code = 21,
+             name = 'G',
+             antiname = 'G',
+             spin = 3,
+             color = 8,
+             mass = 'ZERO',
+             width = 'ZERO',
+             texname = 'G',
+             antitexname = 'G',
+             line = 'curly',
+             charge = 0,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+H = Particle(pdg_code = 25,
+             name = 'H',
+             antiname = 'H',
+             spin = 1,
+             color = 1,
+             mass = 'MH',
+             width = 'WH',
+             texname = '\\phi',
+             antitexname = '\\phi',
+             line = 'dashed',
+             charge = 0,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+phi0 = Particle(pdg_code = 250,
+                name = 'phi0',
+                antiname = 'phi0',
+                spin = 1,
+                color = 1,
+                mass = 'MZ',
+                width = 'Wphi',
+                texname = 'phi0',
+                antitexname = 'phi0',
+                line = 'dashed',
+                GoldstoneBoson = True,
+                charge = 0,
+                LeptonNumber = 0,
+                GhostNumber = 0)
+
+phi__plus__ = Particle(pdg_code = 251,
+                       name = 'phi+',
+                       antiname = 'phi-',
+                       spin = 1,
+                       color = 1,
+                       mass = 'MW',
+                       width = 'Wphi2',
+                       texname = '\\phi^+',
+                       antitexname = '\\phi^+',
+                       line = 'dashed',
+                       GoldstoneBoson = True,
+                       charge = 1,
+                       LeptonNumber = 0,
+                       GhostNumber = 0)
+
+phi__minus__ = phi__plus__.anti()
+

=== added file 'models/uutt_sch_4fermion/vertices.py'
--- models/uutt_sch_4fermion/vertices.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/vertices.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,347 @@
+# This file was automatically created by FeynRules $Revision: 302 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Tue 31 Aug 2010 16:54:46
+
+
+from object_library import all_vertices, Vertex
+import particles as P
+import couplings as C
+import lorentz as L
+
+
+V_1 = Vertex(name = 'V_1',
+             particles = [ P.H, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSS1 ],
+             couplings = {(0,0):C.GC_21})
+
+V_2 = Vertex(name = 'V_2',
+             particles = [ P.G, P.G, P.G ],
+             color = [ 'f(1,2,3)' ],
+             lorentz = [ L.VVV1 ],
+             couplings = {(0,0):C.GC_4})
+
+V_3 = Vertex(name = 'V_3',
+             particles = [ P.G, P.G, P.G, P.G ],
+             color = [ 'f(2,3,\'a1\')*f(\'a1\',1,4)', 'f(2,4,\'a1\')*f(\'a1\',1,3)', 'f(3,4,\'a1\')*f(\'a1\',1,2)' ],
+             lorentz = [ L.VVVV1, L.VVVV3, L.VVVV4 ],
+             couplings = {(1,1):C.GC_6,(2,0):C.GC_6,(0,2):C.GC_6})
+
+V_4 = Vertex(name = 'V_4',
+             particles = [ P.A, P.W__minus__, P.W__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.VVV1 ],
+             couplings = {(0,0):C.GC_16})
+
+V_5 = Vertex(name = 'V_5',
+             particles = [ P.W__minus__, P.W__plus__, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.VVSS1 ],
+             couplings = {(0,0):C.GC_10})
+
+V_6 = Vertex(name = 'V_6',
+             particles = [ P.W__minus__, P.W__plus__, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.VVS1 ],
+             couplings = {(0,0):C.GC_22})
+
+V_7 = Vertex(name = 'V_7',
+             particles = [ P.A, P.A, P.W__minus__, P.W__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.VVVV2 ],
+             couplings = {(0,0):C.GC_18})
+
+V_8 = Vertex(name = 'V_8',
+             particles = [ P.W__minus__, P.W__plus__, P.Z ],
+             color = [ '1' ],
+             lorentz = [ L.VVV1 ],
+             couplings = {(0,0):C.GC_7})
+
+V_9 = Vertex(name = 'V_9',
+             particles = [ P.W__minus__, P.W__minus__, P.W__plus__, P.W__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.VVVV2 ],
+             couplings = {(0,0):C.GC_8})
+
+V_10 = Vertex(name = 'V_10',
+              particles = [ P.A, P.W__minus__, P.W__plus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV5 ],
+              couplings = {(0,0):C.GC_17})
+
+V_11 = Vertex(name = 'V_11',
+              particles = [ P.Z, P.Z, P.H, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_20})
+
+V_12 = Vertex(name = 'V_12',
+              particles = [ P.Z, P.Z, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_23})
+
+V_13 = Vertex(name = 'V_13',
+              particles = [ P.W__minus__, P.W__plus__, P.Z, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV2 ],
+              couplings = {(0,0):C.GC_9})
+
+V_14 = Vertex(name = 'V_14',
+              particles = [ P.d__tilde__, P.d, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_1})
+
+V_15 = Vertex(name = 'V_15',
+              particles = [ P.s__tilde__, P.s, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_1})
+
+V_16 = Vertex(name = 'V_16',
+              particles = [ P.b__tilde__, P.b, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_1})
+
+V_17 = Vertex(name = 'V_17',
+              particles = [ P.e__plus__, P.e__minus__, P.A ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_18 = Vertex(name = 'V_18',
+              particles = [ P.m__plus__, P.m__minus__, P.A ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_19 = Vertex(name = 'V_19',
+              particles = [ P.tt__plus__, P.tt__minus__, P.A ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_20 = Vertex(name = 'V_20',
+              particles = [ P.u__tilde__, P.u, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_21 = Vertex(name = 'V_21',
+              particles = [ P.c__tilde__, P.c, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_22 = Vertex(name = 'V_22',
+              particles = [ P.t__tilde__, P.t, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_23 = Vertex(name = 'V_23',
+              particles = [ P.d__tilde__, P.d, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_24 = Vertex(name = 'V_24',
+              particles = [ P.s__tilde__, P.s, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_25 = Vertex(name = 'V_25',
+              particles = [ P.b__tilde__, P.b, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_26 = Vertex(name = 'V_26',
+              particles = [ P.b__tilde__, P.b, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_24})
+
+V_27 = Vertex(name = 'V_27',
+              particles = [ P.d__tilde__, P.d, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              couplings = {(0,0):C.GC_12,(0,1):C.GC_14})
+
+V_28 = Vertex(name = 'V_28',
+              particles = [ P.s__tilde__, P.s, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              couplings = {(0,0):C.GC_12,(0,1):C.GC_14})
+
+V_29 = Vertex(name = 'V_29',
+              particles = [ P.b__tilde__, P.b, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              couplings = {(0,0):C.GC_12,(0,1):C.GC_14})
+
+V_30 = Vertex(name = 'V_30',
+              particles = [ P.d__tilde__, P.u, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_31 = Vertex(name = 'V_31',
+              particles = [ P.s__tilde__, P.c, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_32 = Vertex(name = 'V_32',
+              particles = [ P.b__tilde__, P.t, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_33 = Vertex(name = 'V_33',
+              particles = [ P.u__tilde__, P.d, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_34 = Vertex(name = 'V_34',
+              particles = [ P.c__tilde__, P.s, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_35 = Vertex(name = 'V_35',
+              particles = [ P.t__tilde__, P.b, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_36 = Vertex(name = 'V_36',
+              particles = [ P.u__tilde__, P.u, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_37 = Vertex(name = 'V_37',
+              particles = [ P.c__tilde__, P.c, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_38 = Vertex(name = 'V_38',
+              particles = [ P.t__tilde__, P.t, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_39 = Vertex(name = 'V_39',
+              particles = [ P.tt__plus__, P.tt__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_26})
+
+V_40 = Vertex(name = 'V_40',
+              particles = [ P.t__tilde__, P.t, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_25})
+
+V_41 = Vertex(name = 'V_41',
+              particles = [ P.e__plus__, P.e__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_12,(0,1):C.GC_15})
+
+V_42 = Vertex(name = 'V_42',
+              particles = [ P.m__plus__, P.m__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_12,(0,1):C.GC_15})
+
+V_43 = Vertex(name = 'V_43',
+              particles = [ P.tt__plus__, P.tt__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_12,(0,1):C.GC_15})
+
+V_44 = Vertex(name = 'V_44',
+              particles = [ P.e__plus__, P.ve, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_45 = Vertex(name = 'V_45',
+              particles = [ P.m__plus__, P.vm, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_46 = Vertex(name = 'V_46',
+              particles = [ P.tt__plus__, P.vt, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_47 = Vertex(name = 'V_47',
+              particles = [ P.ve__tilde__, P.e__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_48 = Vertex(name = 'V_48',
+              particles = [ P.vm__tilde__, P.m__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_49 = Vertex(name = 'V_49',
+              particles = [ P.vt__tilde__, P.tt__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_50 = Vertex(name = 'V_50',
+              particles = [ P.u__tilde__, P.u, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              couplings = {(0,0):C.GC_13,(0,1):C.GC_14})
+
+V_51 = Vertex(name = 'V_51',
+              particles = [ P.c__tilde__, P.c, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              couplings = {(0,0):C.GC_13,(0,1):C.GC_14})
+
+V_52 = Vertex(name = 'V_52',
+              particles = [ P.t__tilde__, P.t, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              couplings = {(0,0):C.GC_13,(0,1):C.GC_14})
+
+V_53 = Vertex(name = 'V_53',
+              particles = [ P.ve__tilde__, P.ve, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_19})
+
+V_54 = Vertex(name = 'V_54',
+              particles = [ P.vm__tilde__, P.vm, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_19})
+
+V_55 = Vertex(name = 'V_55',
+              particles = [ P.vt__tilde__, P.vt, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_19})
+
+V_56 = Vertex(name = 'V_56',
+              particles = [ P.u__tilde__, P.u__tilde__, P.t, P.t ],
+              color = [ 'K6(-1,1,2)*K6Bar(-1,3,4)' ],
+              lorentz = [ L.FFFF1 ],
+              couplings = {(0,0):C.GC_27})
+

=== added file 'models/uutt_sch_4fermion/write_param_card.py'
--- models/uutt_sch_4fermion/write_param_card.py	1970-01-01 00:00:00 +0000
+++ models/uutt_sch_4fermion/write_param_card.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,65 @@
+
+__date__ = "3 june 2010"
+__author__ = '[e-mail]'
+
+class ParamCardWriter(object):
+    
+    header = \
+    """######################################################################\n""" + \
+    """## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################\n""" + \
+    """######################################################################\n"""   
+    
+    def __init__(self, filename, list_of_parameters=None):
+        """write a valid param_card.dat"""
+        
+        if not list_of_parameters:
+            from parameters import all_parameters
+            list_of_parameters = [param for param in all_parameters if \
+                                                       param.nature=='external']
+        
+        self.fsock = open(filename, 'w')
+        self.fsock.write(self.header)
+        
+        self.write_card(list_of_parameters)
+    
+    
+    def write_card(self, all_ext_param):
+        """ """
+        
+      
+        
+        # list all lhablock
+        all_lhablock = set([param.lhablock for param in all_ext_param])
+        
+        # ordonate lhablock alphabeticaly
+        list(all_lhablock).sort()
+        
+        for lhablock in all_lhablock:
+            self.write_block(lhablock)
+            [self.write_param(param, lhablock) for param in all_ext_param if \
+                                                     param.lhablock == lhablock]
+    def write_block(self, name):
+        """ write a comment for a block"""
+        
+        self.fsock.writelines(
+        """\n###################################""" + \
+        """\n## INFORMATION FOR %s""" % name.upper() +\
+        """\n###################################\n"""
+         )
+        if name!='DECAY':
+            self.fsock.write("""Block %s \n""" % name)
+
+    def write_param(self, param, lhablock):
+        
+        lhacode=' '.join(['%3s' % key for key in param.lhacode])
+        if lhablock != 'DECAY':
+            text = """  %s %e # %s \n""" % (lhacode, param.value, param.name ) 
+        else:
+            text = '''DECAY %s %e \n''' % (lhacode, param.value)
+        self.fsock.write(text) 
+            
+            
+if '__main__' == __name__:
+    ParamCardWriter('./param_card.dat')
+    print 'done'
+    

=== added directory 'models/uutt_tch_4fermion'
=== added file 'models/uutt_tch_4fermion/__init__.py'
--- models/uutt_tch_4fermion/__init__.py	1970-01-01 00:00:00 +0000
+++ models/uutt_tch_4fermion/__init__.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,15 @@
+
+import particles
+import couplings
+import lorentz
+import parameters
+import vertices
+import write_param_card
+
+
+all_particles = particles.all_particles
+all_vertices = vertices.all_vertices
+all_couplings = couplings.all_couplings
+all_lorentz = lorentz.all_lorentz
+all_parameters = parameters.all_parameters
+all_functions = function_library.all_functions

=== added file 'models/uutt_tch_4fermion/couplings.py'
--- models/uutt_tch_4fermion/couplings.py	1970-01-01 00:00:00 +0000
+++ models/uutt_tch_4fermion/couplings.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,122 @@
+# This file was automatically created by FeynRules $Revision: 302 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Tue 31 Aug 2010 16:54:46
+
+
+from object_library import all_couplings, Coupling
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+
+
+GC_1 = Coupling(name = 'GC_1',
+                value = '-(ee*complex(0,1))/3.',
+                order = {'QED':1})
+
+GC_2 = Coupling(name = 'GC_2',
+                value = '(2*ee*complex(0,1))/3.',
+                order = {'QED':1})
+
+GC_3 = Coupling(name = 'GC_3',
+                value = '-(ee*complex(0,1))',
+                order = {'QED':1})
+
+GC_4 = Coupling(name = 'GC_4',
+                value = '-G',
+                order = {'QCD':1})
+
+GC_5 = Coupling(name = 'GC_5',
+                value = 'complex(0,1)*G',
+                order = {'QCD':1})
+
+GC_6 = Coupling(name = 'GC_6',
+                value = 'complex(0,1)*G**2',
+                order = {'QCD':2})
+
+GC_7 = Coupling(name = 'GC_7',
+                value = 'cw*complex(0,1)*gw',
+                order = {'QED':1})
+
+GC_8 = Coupling(name = 'GC_8',
+                value = '-(complex(0,1)*gw**2)',
+                order = {'QED':2})
+
+GC_9 = Coupling(name = 'GC_9',
+                value = 'cw**2*complex(0,1)*gw**2',
+                order = {'QED':2})
+
+GC_10 = Coupling(name = 'GC_10',
+                 value = '(ee**2*complex(0,1))/(2.*sw**2)',
+                 order = {'QED':2})
+
+GC_11 = Coupling(name = 'GC_11',
+                 value = '(ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_12 = Coupling(name = 'GC_12',
+                 value = '-(cw*ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_13 = Coupling(name = 'GC_13',
+                 value = '(cw*ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_14 = Coupling(name = 'GC_14',
+                 value = '-(ee*complex(0,1)*sw)/(6.*cw)',
+                 order = {'QED':1})
+
+GC_15 = Coupling(name = 'GC_15',
+                 value = '(ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_16 = Coupling(name = 'GC_16',
+                 value = 'complex(0,1)*gw*sw',
+                 order = {'QED':1})
+
+GC_17 = Coupling(name = 'GC_17',
+                 value = '-2*cw*complex(0,1)*gw**2*sw',
+                 order = {'QED':2})
+
+GC_18 = Coupling(name = 'GC_18',
+                 value = 'complex(0,1)*gw**2*sw**2',
+                 order = {'QED':2})
+
+GC_19 = Coupling(name = 'GC_19',
+                 value = '(cw*ee*complex(0,1))/(2.*sw) + (ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_20 = Coupling(name = 'GC_20',
+                 value = 'ee**2*complex(0,1) + (cw**2*ee**2*complex(0,1))/(2.*sw**2) + (ee**2*complex(0,1)*sw**2)/(2.*cw**2)',
+                 order = {'QED':2})
+
+GC_21 = Coupling(name = 'GC_21',
+                 value = '-6*complex(0,1)*lam*v',
+                 order = {'QED':1})
+
+GC_22 = Coupling(name = 'GC_22',
+                 value = '(ee**2*complex(0,1)*v)/(2.*sw**2)',
+                 order = {'QED':1})
+
+GC_23 = Coupling(name = 'GC_23',
+                 value = 'ee**2*complex(0,1)*v + (cw**2*ee**2*complex(0,1)*v)/(2.*sw**2) + (ee**2*complex(0,1)*sw**2*v)/(2.*cw**2)',
+                 order = {'QED':1})
+
+GC_24 = Coupling(name = 'GC_24',
+                 value = '-((complex(0,1)*yb)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_25 = Coupling(name = 'GC_25',
+                 value = '-((complex(0,1)*yt)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_26 = Coupling(name = 'GC_26',
+                 value = '-((complex(0,1)*ytau)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+# 4-fermion coupling is 2*g^2*prop for t-channel (due to FS symmetry)
+GC_27 = Coupling(name = 'GC_27',
+                 value = 'complex(0,1)*GSUT**2/MSut**2',
+                 order = {'NP':2})
+GC_28 = Coupling(name = 'GC_28',
+                 value = '-complex(0,1)*GSUT**2/MSut**2',
+                 order = {'NP':2})

=== added file 'models/uutt_tch_4fermion/function_library.py'
--- models/uutt_tch_4fermion/function_library.py	1970-01-01 00:00:00 +0000
+++ models/uutt_tch_4fermion/function_library.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,54 @@
+# This file is part of the UFO.
+#
+# This file contains definitions for functions that
+# are extensions of the cmath library, and correspond
+# either to functions that are in cmath, but inconvenient
+# to access from there (e.g. z.conjugate()),
+# or functions that are simply not defined.
+#
+#
+
+__date__ = "22 July 2010"
+__author__ = "[e-mail]"
+
+import cmath
+from object_library import all_functions, Function
+
+#
+# shortcuts for functions from cmath
+#
+
+complexconjugate = Function(name = 'complexconjugate',
+                            arguments = ('z',),
+                            expression = 'z.conjugate()')
+
+
+re = Function(name = 're',
+              arguments = ('z',),
+              expression = 'z.real')
+
+im = Function(name = 'im',
+              arguments = ('z',),
+              expression = 'z.imag')
+
+# New functions (trigonometric)
+
+sec = Function(name = 'sec',
+             arguments = ('z',),
+             expression = '1./cmath.cos(z)')
+
+asec = Function(name = 'asec',
+             arguments = ('z',),
+             expression = 'cmath.acos(1./z)')
+
+csc = Function(name = 'csc',
+             arguments = ('z',),
+             expression = '1./cmath.sin(z)')
+
+acsc = Function(name = 'acsc',
+             arguments = ('z',),
+             expression = 'cmath.asin(1./z)')
+
+
+
+

=== added file 'models/uutt_tch_4fermion/lorentz.py'
--- models/uutt_tch_4fermion/lorentz.py	1970-01-01 00:00:00 +0000
+++ models/uutt_tch_4fermion/lorentz.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,79 @@
+# This file was automatically created by FeynRules $Revision: 302 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Tue 31 Aug 2010 16:54:46
+
+
+from object_library import all_lorentz, Lorentz
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+
+
+SSS1 = Lorentz(name = 'SSS1',
+               spins = [ 1, 1, 1 ],
+               structure = '1')
+
+FFS1 = Lorentz(name = 'FFS1',
+               spins = [ 2, 2, 1 ],
+               structure = 'Identity(1,2)')
+
+FFV1 = Lorentz(name = 'FFV1',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,1)')
+
+FFV2 = Lorentz(name = 'FFV2',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,\'s1\')*ProjM(\'s1\',1)')
+
+FFV3 = Lorentz(name = 'FFV3',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,\'s1\')*ProjM(\'s1\',1) - 2*Gamma(3,2,\'s1\')*ProjP(\'s1\',1)')
+
+FFV4 = Lorentz(name = 'FFV4',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,\'s1\')*ProjM(\'s1\',1) + 2*Gamma(3,2,\'s1\')*ProjP(\'s1\',1)')
+
+FFV5 = Lorentz(name = 'FFV5',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,\'s1\')*ProjM(\'s1\',1) + 4*Gamma(3,2,\'s1\')*ProjP(\'s1\',1)')
+
+VVS1 = Lorentz(name = 'VVS1',
+               spins = [ 3, 3, 1 ],
+               structure = 'Metric(1,2)')
+
+VVV1 = Lorentz(name = 'VVV1',
+               spins = [ 3, 3, 3 ],
+               structure = 'P(3,1)*Metric(1,2) - P(3,2)*Metric(1,2) - P(2,1)*Metric(1,3) + P(2,3)*Metric(1,3) + P(1,2)*Metric(2,3) - P(1,3)*Metric(2,3)')
+
+VVSS1 = Lorentz(name = 'VVSS1',
+                spins = [ 3, 3, 1, 1 ],
+                structure = 'Metric(1,2)')
+
+VVVV1 = Lorentz(name = 'VVVV1',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - Metric(1,3)*Metric(2,4)')
+
+VVVV2 = Lorentz(name = 'VVVV2',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) + Metric(1,3)*Metric(2,4) - 2*Metric(1,2)*Metric(3,4)')
+
+VVVV3 = Lorentz(name = 'VVVV3',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - Metric(1,2)*Metric(3,4)')
+
+VVVV4 = Lorentz(name = 'VVVV4',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,3)*Metric(2,4) - Metric(1,2)*Metric(3,4)')
+
+VVVV5 = Lorentz(name = 'VVVV5',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - (Metric(1,3)*Metric(2,4))/2. - (Metric(1,2)*Metric(3,4))/2.')
+
+FFFF1 = Lorentz(name = 'FFFF1',
+               spins = [ 2, 2, 2, 2 ],
+                structure = 'Identity(1,2)*Identity(3,4)')
+
+FFFF2 = Lorentz(name = 'FFFF2',
+                spins = [ 2, 2, 2, 2 ],
+                structure = 'Identity(1,4)*Identity(3,2)')
+

=== added file 'models/uutt_tch_4fermion/object_library.py'
--- models/uutt_tch_4fermion/object_library.py	1970-01-01 00:00:00 +0000
+++ models/uutt_tch_4fermion/object_library.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,191 @@
+##
+##
+## Feynrules Header
+##
+##
+##
+##
+##
+
+import cmath
+
+
+class FRBaseClass(object):
+    """The class from which all FeynRules classes are derived."""
+
+    require_args = []
+
+    def __init__(self, *args, **options):
+        assert(len(self.require_args) == len (args))
+    
+        for i, name in enumerate(self.require_args):
+            setattr(self, name, args[i])
+    
+        for (option, value) in options.items():
+            setattr(self, option, value)
+
+    def get(self, name):
+        return getattr(self, name)
+    
+    def set(self, name, value):
+        setattr(self, name, value)
+        
+    def get_all(self):
+        """Return a dictionary containing all the information of the object"""
+        return self.__dict__
+
+    def __str__(self):
+        return self.name
+
+
+    def __repr__(self):
+        replacements = [
+            ('+','__plus__'),
+            ('-','__minus__'),
+            ('@','__at__'),
+            ('!','__exclam__'),
+            ('?','__quest__'),
+            ('*','__star__'),
+            ('~','__tilde__')
+            ]
+        text = self.name
+        for orig,sub in replacements:
+            text = text.replace(orig,sub)
+        return text
+
+
+
+all_particles = []
+
+    
+
+class Particle(FRBaseClass):
+    """A standard Particle"""
+
+    require_args=['pdg_code', 'name', 'antiname', 'spin', 'color', 'mass', 'width', 'texname', 'antitexname', 'line', 'charge']
+
+    require_args_all = ['pdg_code', 'name', 'antiname', 'spin', 'color', 'mass', 'width', 'texname', 'antitexname', 'line', 'charge', 'propagating', 'goldstoneboson']
+
+    def __init__(self, pdg_code, name, antiname, spin, color, mass, width, texname,
+                 antitexname, line, charge , propagating=True, goldstoneboson=False, **options):
+
+        args= (pdg_code, name, antiname, spin, color, mass, width, texname,
+                 antitexname, line, float(charge))
+
+        FRBaseClass.__init__(self, *args,  **options)
+
+        global all_particles
+        all_particles.append(self)
+
+        self.propagating = propagating
+        self.goldstoneboson= goldstoneboson
+
+        self.selfconjugate = (name == antiname)
+
+    def anti(self):
+        if self.selfconjugate:
+            raise Exception('%s has no anti particle.' % self.name) 
+        outdic = {}
+        for k,v in self.__dict__.iteritems():
+            if k not in self.require_args_all:                
+                outdic[k] = -v
+        if self.color in [1,8]:
+            newcolor = self.color
+        else:
+            newcolor = -self.color
+                
+        return Particle(-self.pdg_code, self.antiname, self.name, self.spin, newcolor, self.mass, self.width,
+                        self.antitexname, self.texname, self.line, -self.charge, self.propagating, self.goldstoneboson, **outdic)
+
+
+
+all_parameters = []
+
+class Parameter(FRBaseClass):
+
+    require_args=['name', 'nature', 'type', 'value', 'texname']
+
+    def __init__(self, name, nature, type, value, texname, lhablock=None, lhacode=None):
+
+        args = (name,nature,type,value,texname)
+
+        FRBaseClass.__init__(self, *args)
+
+        args=(name,nature,type,value,texname)
+
+        global all_parameters
+        all_parameters.append(self)
+
+        if (lhablock is None or lhacode is None)  and nature == 'external':
+            raise Exception('Need LHA information for external parameter "%s".' % name)
+        self.lhablock = lhablock
+        self.lhacode = lhacode
+
+all_vertices = []
+
+class Vertex(FRBaseClass):
+
+    require_args=['name', 'particles', 'color', 'lorentz', 'couplings']
+
+    def __init__(self, name, particles, color, lorentz, couplings, **opt):
+ 
+        args = (name, particles, color, lorentz, couplings)
+
+        FRBaseClass.__init__(self, *args, **opt)
+
+        args=(particles,color,lorentz,couplings)
+
+        global all_vertices
+        all_vertices.append(self)
+
+all_couplings = []
+
+class Coupling(FRBaseClass):
+
+    require_args=['name', 'value', 'order']
+
+    def __init__(self, name, value, order, **opt):
+
+        args =(name, value, order)	
+        FRBaseClass.__init__(self, *args, **opt)
+        global all_couplings
+        all_couplings.append(self)
+  
+
+
+all_lorentz = []
+
+class Lorentz(FRBaseClass):
+
+    require_args=['name','spins','structure']
+    
+    def __init__(self, name, spins, structure, **opt):
+        args = (name, spins, structure)
+        FRBaseClass.__init__(self, *args, **opt)
+
+        global all_lorentz
+        all_lorentz.append(self)
+
+
+all_functions = []
+
+class Function(object):
+
+    def __init__(self, name, arguments, expression):
+
+        global all_functions
+        all_functions.append(self)
+
+        self.name = name
+        self.arguments = arguments
+        self.expr = expression
+    
+    def __call__(self, *opt):
+
+        for i, arg in enumerate(self.arguments):
+            exec('%s = %s' % (arg, opt[i] ))
+
+        return eval(self.expr)
+
+
+        

=== added file 'models/uutt_tch_4fermion/parameters.py'
--- models/uutt_tch_4fermion/parameters.py	1970-01-01 00:00:00 +0000
+++ models/uutt_tch_4fermion/parameters.py	2011-01-27 02:53:44 +0000
@@ -0,0 +1,237 @@
+# This file was automatically created by FeynRules $Revision: 302 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Tue 31 Aug 2010 16:54:46
+
+
+
+from object_library import all_parameters, Parameter
+
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+aEWM1 = Parameter(name = 'aEWM1',
+                  nature = 'external',
+                  type = 'real',
+                  value = 132.50698,
+                  texname = '\\text{aEWM1}',
+                  lhablock = 'SMINPUTS',
+                  lhacode = [ 1 ])
+
+Gf = Parameter(name = 'Gf',
+               nature = 'external',
+               type = 'real',
+               value = 0.0000116639,
+               texname = 'G_f',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 2 ])
+
+aS = Parameter(name = 'aS',
+               nature = 'external',
+               type = 'real',
+               value = 0.118,
+               texname = '\\text{aS}',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 3 ])
+
+ymb = Parameter(name = 'ymb',
+                nature = 'external',
+                type = 'real',
+                value = 4.2,
+                texname = '\\text{ymb}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 5 ])
+
+ymt = Parameter(name = 'ymt',
+                nature = 'external',
+                type = 'real',
+                value = 164.5,
+                texname = '\\text{ymt}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 6 ])
+
+ymtau = Parameter(name = 'ymtau',
+                  nature = 'external',
+                  type = 'real',
+                  value = 1.777,
+                  texname = '\\text{ymtau}',
+                  lhablock = 'YUKAWA',
+                  lhacode = [ 15 ])
+
+MTA = Parameter(name = 'MTA',
+                nature = 'external',
+                type = 'real',
+                value = 1.777,
+                texname = '\\text{MTA}',
+                lhablock = 'MASS',
+                lhacode = [ 15 ])
+
+MT = Parameter(name = 'MT',
+               nature = 'external',
+               type = 'real',
+               value = 174.3,
+               texname = '\\text{MT}',
+               lhablock = 'MASS',
+               lhacode = [ 6 ])
+
+MB = Parameter(name = 'MB',
+               nature = 'external',
+               type = 'real',
+               value = 4.7,
+               texname = '\\text{MB}',
+               lhablock = 'MASS',
+               lhacode = [ 5 ])
+
+MZ = Parameter(name = 'MZ',
+               nature = 'external',
+               type = 'real',
+               value = 91.188,
+               texname = '\\text{MZ}',
+               lhablock = 'MASS',
+               lhacode = [ 23 ])
+
+MH = Parameter(name = 'MH',
+               nature = 'external',
+               type = 'real',
+               value = 120.,
+               texname = '\\text{MH}',
+               lhablock = 'MASS',
+               lhacode = [ 25 ])
+
+WT = Parameter(name = 'WT',
+               nature = 'external',
+               type = 'real',
+               value = 1.50833649,
+               texname = '\\text{WT}',
+               lhablock = 'DECAY',
+               lhacode = [ 6 ])
+
+WZ = Parameter(name = 'WZ',
+               nature = 'external',
+               type = 'real',
+               value = 2.44140351,
+               texname = '\\text{WZ}',
+               lhablock = 'DECAY',
+               lhacode = [ 23 ])
+
+WW = Parameter(name = 'WW',
+               nature = 'external',
+               type = 'real',
+               value = 2.04759951,
+               texname = '\\text{WW}',
+               lhablock = 'DECAY',
+               lhacode = [ 24 ])
+
+WH = Parameter(name = 'WH',
+               nature = 'external',
+               type = 'real',
+               value = 0.00575308848,
+               texname = '\\text{WH}',
+               lhablock = 'DECAY',
+               lhacode = [ 25 ])
+
+aEW = Parameter(name = 'aEW',
+                nature = 'internal',
+                type = 'real',
+                value = '1/aEWM1',
+                texname = '\\text{aEW}')
+
+G = Parameter(name = 'G',
+              nature = 'internal',
+              type = 'real',
+              value = '2*cmath.sqrt(aS)*cmath.sqrt(cmath.pi)',
+              texname = 'G')
+
+MW = Parameter(name = 'MW',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(MZ**2/2. + cmath.sqrt(MZ**4/4. - (aEW*cmath.pi*MZ**2)/(Gf*cmath.sqrt(2))))',
+               texname = 'M_W')
+
+ee = Parameter(name = 'ee',
+               nature = 'internal',
+               type = 'real',
+               value = '2*cmath.sqrt(aEW)*cmath.sqrt(cmath.pi)',
+               texname = 'e')
+
+sw2 = Parameter(name = 'sw2',
+                nature = 'internal',
+                type = 'real',
+                value = '1 - MW**2/MZ**2',
+                texname = '\\text{sw2}')
+
+cw = Parameter(name = 'cw',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(1 - sw2)',
+               texname = 'c_w')
+
+sw = Parameter(name = 'sw',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(sw2)',
+               texname = 's_w')
+
+g1 = Parameter(name = 'g1',
+               nature = 'internal',
+               type = 'real',
+               value = 'ee/cw',
+               texname = 'g_1')
+
+gw = Parameter(name = 'gw',
+               nature = 'internal',
+               type = 'real',
+               value = 'ee/sw',
+               texname = 'g_w')
+
+v = Parameter(name = 'v',
+              nature = 'internal',
+              type = 'real',
+              value = '(2*MW*sw)/ee',
+              texname = 'v')
+
+lam = Parameter(name = 'lam',
+                nature = 'internal',
+                type = 'real',
+                value = 'MH**2/(2.*v**2)',
+                texname = '\\text{lam}')
+
+yb = Parameter(name = 'yb',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymb*cmath.sqrt(2))/v',
+               texname = '\\text{yb}')
+
+yt = Parameter(name = 'yt',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymt*cmath.sqrt(2))/v',
+               texname = '\\text{yt}')
+
+ytau = Parameter(name = 'ytau',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '(ymtau*cmath.sqrt(2))/v',
+                 texname = '\\text{ytau}')
+
+muH = Parameter(name = 'muH',
+                nature = 'internal',
+                type = 'real',
+                value = 'cmath.sqrt(lam*v**2)',
+                texname = '\\mu ')
+
+GSUT = Parameter(name = 'GSUT',
+                  nature = 'external',
+                  type = 'real',
+                  value = 0.1,
+                  texname = '\\text{GSUT}',
+                  lhablock = 'BSMINPUTS',
+                  lhacode = [ 1 ])
+
+MSut = Parameter(name = 'MSut',
+                nature = 'external',
+                type = 'real',
+                value = 1.e5,
+                texname = '\\text{Mut}',
+                lhablock = 'BSMINPUTS',
+                lhacode = [ 2 ])