Merge lp:~madteam/mg5amcnlo/madevent_output into lp:~madteam/mg5amcnlo/trunk

[Johan Alwall]

* Decay chain implementation (by Johan).
* MadEvent output (by Johan and Michel).
* A number of bugs and minor problems fixed, both relating to Helas objects and color.

I would recommend that you divide the revision work, e.g. Tim and Michel focus on the decay chains (mainly in helas_objects.py) while Fabio and Olivier take care of the MadEvent output (mainly in export_v4.py).

[lp:~madteam/mg5amcnlo/madevent_output]

54. By Johan Alwall <alwall@alwall-laptop> on 2010-04-19

Added possibility for comments starting with # in the command interface

55. By Johan Alwall <alwall@alwall-laptop> on 2010-04-19

Added possibility to add process to existing processes. Note that for decay chain processes, can not do export twice.

56. By Johan Alwall <alwall@alwall-laptop> on 2010-04-19

Fixed problem with adding processes after decay chain export

57. By Johan Alwall <alwall@alwall-laptop> on 2010-04-20

Fixed error discovered by Olivier when there are no valid decays

58. By Johan Alwall <alwall@alwall-laptop> on 2010-04-20

Fixed path for generate_subprocess_directory_v4_madevent, added writeout of maxamps.inc file

59. By Johan Alwall <alwall@alwall-laptop> on 2010-04-20

Added help line for export v4madevent

60. By Johan Alwall <alwall@alwall-laptop> on 2010-04-21

Updated MadEvent matrix.f output according to new standard in MG/ME v. 4.4.36

61. By Johan Alwall <alwall@alwall-laptop> on 2010-04-21

Updated MadEvent auto_dsig.f output according to new standard in MG/ME v. 4.4.36

[Michel Herquet (herquet)]

Honestly I looked at it for quite some time, and I can only say… fantastic work Johan! Every single line of code is very well documented, and even though there are many tricky parts (some of them I am not sure I understand, but it would require more time), I think the code is quite DRY. I haven't yet ran manual test, but I already give my formal approval for the quality of the code. Once again, great job...

review: Approve

[lp:~madteam/mg5amcnlo/madevent_output]

62. By Johan Alwall <alwall@alwall-laptop> on 2010-04-21

Fixed bug in propagator in configs.inc for t-channel case

63. By Johan Alwall <alwall@alwall-laptop> on 2010-04-21

Fixed bug and test

64. By Johan Alwall <alwall@alwall-laptop> on 2010-04-21

Merged with Oliviers changes

65. By Johan Alwall <alwall@alwall-laptop> on 2010-04-22

Fixed same error as rev 62 but for wavefunction - should work but cannot find example to test

66. By Johan Alwall <alwall@alwall-laptop> on 2010-04-23

Added emptyline() so empty line does not repeat previous line

67. By Michel Herquet on 2010-04-23

fixed ordering of element from the color basis dictionary when creating the color matrix and outputting JAMPs

68. By Johan Alwall <alwall@alwall-laptop> on 2010-04-24

First pass at initial state mirror process optimization

69. By Johan Alwall <alwall@alwall-laptop> on 2010-04-24

Reverted to rev. 67, since mirror process optimization turned out to not give optimal subproc combination

70. By Michel Herquet on 2010-04-26

Fix ordering of color matrix (2d attempt)

71. By Johan Alwall <alwall@alwall-laptop> on 2010-04-26

Fixed the tests after Michels fix

72. By Michel Herquet on 2010-04-26

Trying to fix col matrix ordering, 3rd attempt...

73. By Michel Herquet on 2010-04-26

Fixed VERY strange bug where a color basis element was shared between two different instances (not reset in __init__)

74. By Johan Alwall <alwall@alwall-laptop> on 2010-04-27

First pass on >50% memory optimization. Still need to fix tests.

75. By Johan Alwall <alwall@alwall-laptop> on 2010-04-27

Second pass on >50% memory optimization, still need to fix tests.

76. By Johan Alwall <alwall@alwall-laptop> on 2010-04-28

All tests fixed for memory optimization, except one drawing test.

77. By Johan Alwall <alwall@alwall-laptop> on 2010-04-28

Oops, fixed slight bug

78. By Johan Alwall <alwall@alwall-laptop> on 2010-04-28

Fixed memory for decay chains, in particular problems with deepcopied model

79. By Johan Alwall <alwall@alwall-laptop> on 2010-04-28

Merged with new version info

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

I've still not work intensively on the reviewing. But this is a first comment.

In the files export_v4, you have define a lot of long string.
I think that it will be more clear if you put those files content in some files in the
Template_file directory.

I suggest to create one files for each of those long string (with the associate name of mg4)
and then load them via the string.Template class.

Olivier

PS: Off course this is only a suggestion.

In this sense this breaks some rules of how to

Why do not use the Templates

review: Needs Information (suggestion)

[Johan Alwall (johan-alwall)]

Hi Olivier,

I just hadn't thought about the possibility to use templates. I'll try it out and see how it works. The only worry is that it might make the code slightly less clear, but on the other hand more concise so probably the pros are larger than the cons.

Johan

[lp:~madteam/mg5amcnlo/madevent_output]

80. By Johan Alwall <alwall@alwall-laptop> on 2010-04-28

Fix cmd_interface error

81. By Johan Alwall <alwall@alwall-laptop> on 2010-04-30

Did Oliviers suggested move of long file strings in export_v4.py to template_files/. Improved memory optimization yet a couple of notches.

[Johan Alwall (johan-alwall)]

Ok I did as Olivier suggested and moved the long file strings (matrix.f and auto_dsig.f) to template files in the iolibs/template_files directory. Also improved memory optimization another couple of notches (by reducing the number of Legs and Vertices kept in memory).

[lp:~madteam/mg5amcnlo/madevent_output]

82. By Olivier Mattelaer on 2010-05-01

correct the test problem with decay chains. and update the script creating the decay chains

83. By Olivier Mattelaer on 2010-05-01

review of Olivier:
minor modification (mainly remove unused variable)

[Olivier Mattelaer (olivier-mattelaer)]

Thanks Johan for your modification,

I've made a more deep review of the export_v4.py.
And I should say that everything is perfect. Thanks a lot for this beautifull work!!!

I've suppress/modify some lines, mainly because they were simply not use. This doesn't change anything deep. And the test passes afterwards.

Also what do you think to put the class FortranWriter outside export_v4.py.
I think this will be more clear if this class is not in export_v4. but in some specific module
file_writer.py (where we will insert the writer for C++ files for examples).

Thanks a lot,

Olivier

review: Approve

[lp:~madteam/mg5amcnlo/madevent_output]

84. By Olivier Mattelaer on 2010-05-01

minor modification

[Olivier Mattelaer (olivier-mattelaer)]

Hi Johan,

Could you also look at this:
generate u u~ > t t~ / a $ a QED=1 , (t > w+ b / a $ a , w+ > l+ vl / a $ a  ) , (t~ > w- b~ / a $ a , w- > l- vl~ / a $ a  )
INFO:diagram_generation:Trying process: u u~ > t t~ / a a QED=1
...

It seems that the $ request is not take correctly take into account (at least for the information message)

review: Needs Fixing

[Tim Stelzer (tstelzer)]

Very impressive. If I understand things correctly, this works VERY differently from V4. You are separately calculating the production M.E. and the decay M.E. and then "attaching" them. That will certainly have great speed advantages.... but required a lot of very nice coding.

It is very well documented, and what I followed in (helas_objects.py) looked good. But, I have not had the opportunity to run a battery of tests. (Have you already done that? is it part of the test suite?)

That is my only remaining activity before approving.

review: Needs Information

[lp:~madteam/mg5amcnlo/madevent_output]

85. By Johan Alwall <alwall@alwall-laptop> on 2010-05-05

Fixed broken small MadEvent files, added tests

[Johan Alwall (johan-alwall)]

> Could you also look at this:
> generate u u~ > t t~ / a $ a QED=1 , (t > w+ b / a $ a , w+ > l+ vl / a $ a  )
> , (t~ > w- b~ / a $ a , w- > l- vl~ / a $ a  )
> INFO:diagram_generation:Trying process: u u~ > t t~ / a a QED=1

Hi Olivier,

In fact the order is fixed to be $ before /, as indicated in the help_generate. Please go ahead and change this is you want - at present the order of all features is complately fixed.

Best,
Johan

[Olivier Mattelaer (olivier-mattelaer)]

Ok Thanks,

I'll check that the converter mg4-> mg5 respect this convention.

review: Approve

[lp:~madteam/mg5amcnlo/madevent_output]

86. By Johan Alwall <alwall@alwall-laptop> on 2010-05-06

Changed line in auto_dsig from core process to complete decay process

87. By Johan Alwall <alwall@alwall-laptop> on 2010-05-07

Allow @ flag anywhere, allow arbitrary order of $ and /

88. By Johan Alwall <alwall@alwall-laptop> on 2010-05-07

Woops, fixed bug when @N had no characters after it

89. By Johan Alwall <alwall@alwall-laptop> on 2010-05-11

Added some better error messages to cmd_interface

=== modified file '.pydevproject'
--- .pydevproject	2010-02-12 13:09:25 +0000
+++ .pydevproject	2010-05-07 07:20:46 +0000
@@ -1,9 +1,10 @@
-<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<?xml version="1.0" encoding="UTF-8"?>
 <?eclipse-pydev version="1.0"?>
+
 <pydev_project>
 <pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.6</pydev_property>
 <pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Default</pydev_property>
 <pydev_pathproperty name="org.python.pydev.PROJECT_SOURCE_PATH">
-<path>/MG-trunk</path>
+<path>/madevent_output</path>
 </pydev_pathproperty>
 </pydev_project>

=== modified file 'madgraph/VERSION'
--- madgraph/VERSION	2010-02-12 13:07:52 +0000
+++ madgraph/VERSION	2010-05-07 07:20:46 +0000
@@ -1,2 +1,2 @@
-version = 0.3.0
-date = 2010-02-12
+version = 0.4.0
+date = 2010-05-01

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2010-02-12 12:53:59 +0000
+++ madgraph/core/base_objects.py	2010-05-07 07:20:46 +0000
@@ -59,7 +59,7 @@
 
         try:
             return dict.__getitem__(self, name)
-        except:
+        except KeyError:
             self.is_valid_prop(name) #raise the correct error
 
 
@@ -81,16 +81,9 @@
 
         return True
 
-    def __getitem__(self, name):
-        """ Get the value of the property name."""
-
-        if self.is_valid_prop(name):
-            return dict.__getitem__(self, name)
-
     def get(self, name):
         """Get the value of the property name."""
 
-        #if self.is_valid_prop(name): #done automaticaly in __getitem__
         return self[name]
 
     def set(self, name, value):
@@ -332,7 +325,16 @@
         """Return the color code with a correct minus sign"""
 
         if not self['is_part'] and self['color'] in [3, 6]:
-            return -self['color']
+            return - self['color']
+        else:
+            return self['color']
+
+    def get_anti_color(self):
+        """Return the color code of the antiparticle with a correct minus sign
+        """
+
+        if self['is_part'] and self['color'] in [3, 6]:
+            return - self['color']
         else:
             return self['color']
 
@@ -369,6 +371,11 @@
 
         return self['spin'] % 2 == 0
 
+    def is_boson(self):
+        """Returns True if this is a boson, False if fermion"""
+
+        return self['spin'] % 2 == 1
+
 #===============================================================================
 # ParticleList
 #===============================================================================
@@ -799,7 +806,8 @@
 
         self['id'] = 0
         self['number'] = 0
-        self['state'] = 'final'
+        # True = final, False = initial (boolean to save memory)
+        self['state'] = True
         self['from_group'] = True
 
     def filter(self, name, value):
@@ -811,17 +819,13 @@
                         "%s is not a valid integer for leg id" % str(value)
 
         if name == 'state':
-            if not isinstance(value, str):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid string for leg state" % \
-                                                                    str(value)
-            if value not in ['initial', 'final']:
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid leg state (initial|final)" % \
+            if not isinstance(value, bool):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid leg state (True|False)" % \
                                                                     str(value)
 
         if name == 'from_group':
-            if not isinstance(value, bool):
+            if not isinstance(value, bool) and value != None:
                 raise self.PhysicsObjectError, \
                         "%s is not a valid boolean for leg flag from_group" % \
                                                                     str(value)
@@ -853,8 +857,8 @@
 
         part = model.get('particle_dict')[self['id']]
         return part.is_fermion() and \
-               (self.get('state') == 'initial' and part.get('is_part') or \
-                self.get('state') == 'final' and not part.get('is_part'))
+               (self.get('state') == False and part.get('is_part') or \
+                self.get('state') == True and not part.get('is_part'))
 
     def is_outgoing_fermion(self, model):
         """Returns True if leg is an outgoing fermion, i.e., initial
@@ -866,8 +870,8 @@
 
         part = model.get('particle_dict')[self['id']]
         return part.is_fermion() and \
-               (self.get('state') == 'final' and part.get('is_part') or \
-                self.get('state') == 'initial' and not part.get('is_part'))
+               (self.get('state') == True and part.get('is_part') or \
+                self.get('state') == False and not part.get('is_part'))
 
 #===============================================================================
 # LegList
@@ -906,10 +910,24 @@
         else:
             return False
 
-    def can_combine_to_0(self, ref_dict_to0):
+    def can_combine_to_0(self, ref_dict_to0, is_decay_chain=False):
         """If has at least two 'from_group' True and in ref_dict_to0,
-           return the vertex (with id from ref_dict_to0), otherwise return None
-           """
+        
+        return the vertex (with id from ref_dict_to0), otherwise return None
+
+        If is_decay_chain = True, we only allow clustering of the
+        initial leg, since we want this to be the last wavefunction to
+        be evaluated.
+        """
+        if is_decay_chain:
+            # Special treatment - here we only allow combination to 0
+            # if the initial leg (marked by from_group = None) is
+            # unclustered, since we want this to stay until the very
+            # end.
+            return any(leg.get('from_group') == None for leg in self) and \
+                   ref_dict_to0.has_key(tuple(sorted([leg.get('id') \
+                                                      for leg in self])))
+
         if self.minimum_two_from_group():
             return ref_dict_to0.has_key(tuple(sorted([leg.get('id') for leg in self])))
         else:
@@ -926,7 +944,7 @@
             return res
 
         for leg in self:
-            if leg.get('state') == 'initial':
+            if leg.get('state') == False:
                 res.append(model.get('particle_dict')[leg.get('id')].get_anti_pdg_code())
             else:
                 res.append(leg.get('id'))
@@ -944,7 +962,7 @@
         """Default values for all properties"""
 
         self['ids'] = []
-        self['state'] = 'final'
+        self['state'] = True
 
     def filter(self, name, value):
         """Filter for valid multileg property values."""
@@ -959,11 +977,7 @@
                           "%s is not a valid list of integers" % str(value)
 
         if name == 'state':
-            if not isinstance(value, str):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid string for leg state" % \
-                                                                    str(value)
-            if value not in ['initial', 'final']:
+            if not isinstance(value, bool):
                 raise self.PhysicsObjectError, \
                         "%s is not a valid leg state (initial|final)" % \
                                                                     str(value)
@@ -1030,8 +1044,8 @@
 
         if ninitial == 1:
             # For one initial particle, all legs are s-channel
-            # Only need to flip particle id if state is 'initial'
-            if leg.get('state') == 'final':
+            # Only need to flip particle id if state is False
+            if leg.get('state') == True:
                 return leg.get('id')
             else:
                 return model.get('particle_dict')[leg.get('id')].\
@@ -1039,7 +1053,7 @@
 
         # Number of initial particles is at least 2
         if self.get('id') == 0 or \
-           leg.get('state') == 'initial':
+           leg.get('state') == False:
             # identity vertex or t-channel particle
             return 0
 
@@ -1053,7 +1067,7 @@
 
         ## Check if the other legs are initial or final.
         ## If the latter, return leg id, if the former, return -leg id
-        #if self.get('legs')[0].get('state') == 'final':
+        #if self.get('legs')[0].get('state') == True:
         #    return leg.get('id')
         #else:
         #    return model.get('particle_dict')[leg.get('id')].\
@@ -1146,11 +1160,12 @@
 
         return isinstance(obj, Diagram)
 
-    def nice_string(self):
+    def nice_string(self, indent=0):
         """Returns a nicely formatted string"""
-        mystr = str(len(self)) + ' diagrams:\n'
-        for diag in self:
-            mystr = mystr + "  " + diag.nice_string() + '\n'
+        mystr = " " * indent + str(len(self)) + ' diagrams:\n'
+        for i, diag in enumerate(self):
+            mystr = mystr + " " * indent + str(i+1) + "  " + \
+                    diag.nice_string() + '\n'
         return mystr[:-1]
 
 
@@ -1175,6 +1190,9 @@
         self['required_s_channels'] = []
         self['forbidden_s_channels'] = []
         self['forbidden_particles'] = []
+        self['is_decay_chain'] = False
+        # Decay chain processes associated with this process
+        self['decay_chains'] = ProcessList()
 
     def filter(self, name, value):
         """Filter for valid process property values."""
@@ -1220,25 +1238,43 @@
                     raise self.PhysicsObjectError, \
                       "Forbidden particles should have a positive PDG code" % str(value)
 
+        if name == 'is_decay_chain':
+            if not isinstance(value, bool):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid bool" % str(value)
+
+        if name == 'decay_chains':
+            if not isinstance(value, ProcessList):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid ProcessList" % str(value)
+
         return True
 
+    def set(self, name, value):
+        """Special set for forbidden particles - set to abs value."""
+
+        if name == 'forbidden_particles':
+            value = [abs(i) for i in value]
+
+        return super(Process, self).set(name, value) # call the mother routine
+
     def get_sorted_keys(self):
         """Return process property names as a nicely sorted list."""
 
         return ['legs', 'orders', 'model', 'id',
                 'required_s_channels', 'forbidden_s_channels',
-                'forbidden_particles']
+                'forbidden_particles', 'is_decay_chain', 'decay_chains']
 
-    def nice_string(self):
+    def nice_string(self, indent=0):
         """Returns a nicely formated string about current process
         content"""
 
-        mystr = "Process: "
+        mystr = " " * indent + "Process: "
         prevleg = None
         for leg in self['legs']:
             mypart = self['model'].get('particle_dict')[leg['id']]
-            if prevleg and prevleg['state'] == 'initial' \
-                   and leg['state'] == 'final':
+            if prevleg and prevleg['state'] == False \
+                   and leg['state'] == True:
                 # Separate initial and final legs by ">"
                 mystr = mystr + '> '
                 # Add required s-channels
@@ -1267,10 +1303,35 @@
                 mystr = mystr + forbpart.get_name() + ' '
 
         if self['orders']:
-            mystr = mystr[:-1] + "\n"
-            mystr = mystr + 'Orders: '
-            mystr = mystr + ", ".join([key + '=' + repr(self['orders'][key]) \
+            mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
                        for key in self['orders']]) + ' '
+
+        # Remove last space
+        mystr = mystr[:-1]
+
+        if not self.get('decay_chains'):
+            return mystr
+
+        for decay in self['decay_chains']:
+            mystr = mystr + '\n' + \
+                    decay.nice_string(indent + 2).replace('Process', 'Decay')
+
+        return mystr
+
+    def base_string(self):
+        """Returns a string containing only the basic process (w/ decays)."""
+
+        mystr = ""
+        prevleg = None
+        for leg in self.get_legs_with_decays():
+            mypart = self['model'].get('particle_dict')[leg['id']]
+            if prevleg and prevleg['state'] == False \
+                   and leg['state'] == True:
+                # Separate initial and final legs by ">"
+                mystr = mystr + '> '
+            mystr = mystr + mypart.get_name() + ' '
+            prevleg = leg
+
         # Remove last space
         return mystr[:-1]
 
@@ -1283,8 +1344,8 @@
         prevleg = None
         for leg in self['legs']:
             mypart = self['model'].get('particle_dict')[leg['id']]
-            if prevleg and prevleg['state'] == 'initial' \
-                   and leg['state'] == 'final':
+            if prevleg and prevleg['state'] == False \
+                   and leg['state'] == True:
                 # Separate initial and final legs by ">"
                 mystr = mystr + '_'
                 # Add required s-channels
@@ -1311,6 +1372,9 @@
         # Just to be safe, remove all spaces
         mystr = mystr.replace(' ', '')
 
+        for decay in self.get('decay_chains'):
+            mystr = mystr + decay.shell_string().replace("%d_" % decay.get('id'),
+                                                        "_", 1)
         return mystr
 
     def shell_string_v4(self):
@@ -1319,10 +1383,10 @@
 
         mystr = "%d_" % self['id']
         prevleg = None
-        for leg in self['legs']:
+        for leg in self.get_legs_with_decays():
             mypart = self['model'].get('particle_dict')[leg['id']]
-            if prevleg and prevleg['state'] == 'initial' \
-                   and leg['state'] == 'final':
+            if prevleg and prevleg['state'] == False \
+                   and leg['state'] == True:
                 # Separate initial and final legs by ">"
                 mystr = mystr + '_'
             if mypart['is_part']:
@@ -1338,6 +1402,54 @@
 
         return mystr
 
+    # Helper functions
+
+    def get_ninitial(self):
+        """Gives number of initial state particles"""
+
+        return len(filter(lambda leg: leg.get('state') == False,
+                           self.get('legs')))
+
+    def get_initial_ids(self):
+        """Gives the pdg codes for initial state particles"""
+
+        return [leg.get('id') for leg in \
+                filter(lambda leg: leg.get('state') == False,
+                       self.get('legs'))]
+
+    def get_initial_pdg(self, number):
+        """Return the pdg codes for initial state particles for beam number"""
+
+        return filter(lambda leg: leg.get('state') == False and\
+                       leg.get('number') == number,
+                       self.get('legs'))[0].get('id')
+
+    def get_final_legs(self):
+        """Gives the pdg codes for initial state particles"""
+
+        return filter(lambda leg: leg.get('state') == True,
+                       self.get('legs'))
+    
+    def get_legs_with_decays(self):
+        """Return process with all decay chains substituted in."""
+
+        legs = copy.deepcopy(self.get('legs'))
+        if self.get('is_decay_chain'):
+            legs.pop(0)
+        ileg = 0
+        for decay in self.get('decay_chains'):
+            while legs[ileg].get('state') == False or \
+                      legs[ileg].get('id') != decay.get('legs')[0].get('id'):
+                ileg = ileg + 1
+            decay_legs = decay.get_legs_with_decays()
+            legs = legs[:ileg] + decay_legs + legs[ileg+1:]
+            ileg = ileg + len(decay_legs)
+
+        for ileg, leg in enumerate(legs):
+            leg.set('number', ileg + 1)
+            
+        return LegList(legs)
+
 #===============================================================================
 # ProcessList
 #===============================================================================
@@ -1353,7 +1465,7 @@
 #===============================================================================
 # ProcessDefinition
 #===============================================================================
-class ProcessDefinition(PhysicsObject):
+class ProcessDefinition(Process):
     """ProcessDefinition: list of multilegs (ordered)
                           dictionary of orders
                           model
@@ -1363,13 +1475,11 @@
     def default_setup(self):
         """Default values for all properties"""
 
+        super(ProcessDefinition, self).default_setup()
+
         self['legs'] = MultiLegList()
-        self['orders'] = {}
-        self['model'] = Model()
-        self['id'] = 0
-        self['required_s_channels'] = []
-        self['forbidden_s_channels'] = []
-        self['forbidden_particles'] = []
+        # Decay chain processes associated with this process
+        self['decay_chains'] = ProcessDefinitionList()
 
     def filter(self, name, value):
         """Filter for valid process property values."""
@@ -1378,51 +1488,20 @@
             if not isinstance(value, MultiLegList):
                 raise self.PhysicsObjectError, \
                         "%s is not a valid MultiLegList object" % str(value)
-        if name == 'orders':
-            Interaction.filter(Interaction(), 'orders', value)
-
-        if name == 'model':
-            if not isinstance(value, Model):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid Model object" % str(value)
-        if name is 'id':
-            if not isinstance(value, int):
-                raise self.PhysicsObjectError, \
-                    "Process id %s is not an integer" % repr(value)
-
-        if name in ['required_s_channels',
-                    'forbidden_s_channels']:
-            if not isinstance(value, list):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid list" % str(value)
-            for i in value:
-                if not isinstance(i, int):
-                    raise self.PhysicsObjectError, \
-                          "%s is not a valid list of integers" % str(value)
-                if i == 0:
-                    raise self.PhysicsObjectError, \
-                      "Not valid PDG code %d for s-channel particle" % str(value)
-
-        if name == 'forbidden_particles':
-            if not isinstance(value, list):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid list" % str(value)
-            for i in value:
-                if not isinstance(i, int):
-                    raise self.PhysicsObjectError, \
-                          "%s is not a valid list of integers" % str(value)
-                if i <= 0:
-                    raise self.PhysicsObjectError, \
-                      "Forbidden particles should have a positive PDG code" % str(value)
+        elif name == 'decay_chains':
+            if not isinstance(value, ProcessDefinitionList):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid ProcessDefinitionList" % str(value)
+
+        else:
+            return super(ProcessDefinition, self).filter(name, value)
 
         return True
 
     def get_sorted_keys(self):
         """Return process property names as a nicely sorted list."""
 
-        return ['legs', 'orders', 'model', 'id',
-                'required_s_channels', 'forbidden_s_channels',
-                'forbidden_particles']
+        return super(ProcessDefinition, self).get_sorted_keys()
 
 #===============================================================================
 # ProcessDefinitionList

=== modified file 'madgraph/core/color_algebra.py'
--- madgraph/core/color_algebra.py	2010-02-11 13:27:47 +0000
+++ madgraph/core/color_algebra.py	2010-05-07 07:20:46 +0000
@@ -78,6 +78,7 @@
 
     __copy__ = create_copy
 
+
 #===============================================================================
 # Tr
 #===============================================================================
@@ -201,11 +202,10 @@
 
         return None
 
-    def pair_simplify(self, col_obj, simplify_T_product=False):
+    def pair_simplify(self, col_obj):
         """Implement T(a,...,i,j)T(b,...,j,k) = T(a,...,b,...,i,k)
         and T(a,x,b,i,j)T(c,x,d,k,l) = 1/2(T(a,d,i,l)T(c,b,k,j)    
-                                        -1/Nc T(a,b,i,j)T(c,d,k,l))
-        but only if the simplify_T_product tag is True."""
+                                        -1/Nc T(a,b,i,j)T(c,d,k,l))."""
 
         if isinstance(col_obj, T):
             ij1 = self[-2:]
@@ -219,30 +219,29 @@
                                                                ij2[1]])))])])
             # T(a,x,b,i,j)T(c,x,d,k,l) = 1/2(T(a,d,i,l)T(c,b,k,j)    
             #                          -1/Nc T(a,b,i,j)T(c,d,k,l))
-            if simplify_T_product:
-                for i1, index1 in enumerate(self[:-2]):
-                    for i2, index2 in enumerate(col_obj[:-2]):
-                        if index1 == index2:
-                            a = self[:i1]
-                            b = self[i1 + 1:-2]
-                            c = col_obj[:i2]
-                            d = col_obj[i2 + 1:-2]
-                            col_str1 = ColorString([T(*(a + d + \
-                                                       array.array('i',
-                                                        [ij1[0], ij2[1]]))),
-                                                   T(*(c + b + \
-                                                       array.array('i',
-                                                        [ij2[0], ij1[1]])))])
-                            col_str2 = ColorString([T(*(a + b + \
-                                                       array.array('i',
-                                                        [ij1[0], ij1[1]]))),
-                                                   T(*(c + d + \
-                                                       array.array('i',
-                                                        [ij2[0], ij2[1]])))])
-                            col_str1.coeff = fractions.Fraction(1, 2)
-                            col_str2.coeff = fractions.Fraction(-1, 2)
-                            col_str2.Nc_power = -1
-                            return ColorFactor([col_str1, col_str2])
+            for i1, index1 in enumerate(self[:-2]):
+                for i2, index2 in enumerate(col_obj[:-2]):
+                    if index1 == index2:
+                        a = self[:i1]
+                        b = self[i1 + 1:-2]
+                        c = col_obj[:i2]
+                        d = col_obj[i2 + 1:-2]
+                        col_str1 = ColorString([T(*(a + d + \
+                                                   array.array('i',
+                                                    [ij1[0], ij2[1]]))),
+                                               T(*(c + b + \
+                                                   array.array('i',
+                                                    [ij2[0], ij1[1]])))])
+                        col_str2 = ColorString([T(*(a + b + \
+                                                   array.array('i',
+                                                    [ij1[0], ij1[1]]))),
+                                               T(*(c + d + \
+                                                   array.array('i',
+                                                    [ij2[0], ij2[1]])))])
+                        col_str1.coeff = fractions.Fraction(1, 2)
+                        col_str2.coeff = fractions.Fraction(-1, 2)
+                        col_str2.Nc_power = -1
+                        return ColorFactor([col_str1, col_str2])
 
     def complex_conjugate(self):
         """Complex conjugation. Overwritten here because the two last indices
@@ -631,3 +630,4 @@
 
 
 
+

=== modified file 'madgraph/core/color_amp.py'
--- madgraph/core/color_amp.py	2010-02-12 09:52:24 +0000
+++ madgraph/core/color_amp.py	2010-05-07 07:20:46 +0000
@@ -43,6 +43,11 @@
     # Dictionary store the raw colorize information
     _list_color_dict = []
 
+    class ColorBasisError(Exception):
+        """Exception raised if an error occurs in the definition
+        or the execution of a color basis object."""
+        pass
+
     def colorize(self, diagram, model):
         """Takes a diagram and a model and outputs a dictionary with keys being
         color coefficient index tuples and values a color string (before 
@@ -76,6 +81,7 @@
                         for cs in res_dict.values()]):
                     res_dict = {}
                 # Return since this must be the last vertex
+
                 return res_dict
 
         # NORMAL VERTICES WITH ID != 0 -----------------------------------------
@@ -86,6 +92,7 @@
         if all([cs == color_algebra.ColorString() \
                         for cs in res_dict.values()]):
             res_dict = {}
+
         return res_dict
 
     def add_vertex(self, vertex, diagram, model,
@@ -95,16 +102,14 @@
         If the id0_rep list is not None, perform the requested replacement on the
         last leg number before going further."""
 
-        # Create a list of pdg codes entering the vertex ordered as in
-        # interactions.py
-        list_pdg = [part.get_pdg_code() for part in \
-               model.get_interaction(vertex.get('id')).get('particles')]
-        # Create a dictionary pdg code --> leg(s)
-        dict_pdg_leg = {}
+        # Create a list of (color,leg number) pairs for the vertex, where color
+        # can be negative for anti particles
+
+        color_num_pairs = []
 
         for index, leg in enumerate(vertex.get('legs')):
             curr_num = leg.get('number')
-            curr_pdg = leg.get('id')
+            curr_color = model.get('particle_dict')[leg.get('id')].get_color()
 
             # If this is the next-to-last vertex and the last vertex is
             # the special identity id=0, start by applying the replacement rule
@@ -114,12 +119,12 @@
                     curr_num = id0_rep[id0_rep.index(curr_num) - 1]
 
             # If this is the last leg and not the last vertex 
-            # flip part/antipart. If it is not the last, AND not the next-to-last
+            # flip color. If it is not the last, AND not the next-to-last
             # before an id=0 vertex, replace last index by a new summed index.
             if index == len(vertex.get('legs')) - 1 and \
                 vertex != diagram.get('vertices')[-1]:
-                curr_pdg = \
-                    model.get('particle_dict')[curr_pdg].get_anti_pdg_code()
+                curr_color = \
+                    model.get('particle_dict')[leg.get('id')].get_anti_color()
                 if not id0_rep:
                     repl_dict[curr_num] = min_index
                     min_index = min_index - 1
@@ -130,21 +135,23 @@
             except KeyError:
                 pass
 
-            try:
-                dict_pdg_leg[curr_pdg].append(curr_num)
-            except KeyError:
-                dict_pdg_leg[curr_pdg] = [curr_num]
+            # Discard color singlets
+            if curr_color != 1:
+                color_num_pairs.append((curr_color, curr_num))
+
+        # Order the color/number pairs according to increasing color (assumed
+        # to be the ordering chose in interactions.py). For identical colors,
+        # keep the normal leg ordering.
+        color_num_pairs = sorted(color_num_pairs, lambda p1, p2:p1[0] - p2[0])
 
         # Create a list of associated leg number following the same order
-        list_numbers = []
-        for pdg_code in list_pdg:
-            list_numbers.append(dict_pdg_leg[pdg_code].pop(0))
+        list_numbers = [p[1] for p in color_num_pairs]
+
         # ... and the associated dictionary for replacement
         match_dict = dict(enumerate(list_numbers))
 
         # Update the result dict using the current vertex ColorString object
         # If more than one, create different entries
-
         inter_color = model.get_interaction(vertex['id'])['color']
 
         # For colorless vertices, return a copy of res_dict
@@ -285,6 +292,14 @@
         or 1 arguments). If one arguments is given, it's interpreted as 
         an amplitude."""
 
+        dict.__init__(self)
+
+        # Dictionary to save simplifications already done in a canonical form
+        self._canonical_dict = {}
+
+        # Dictionary store the raw colorize information
+        self._list_color_dict = []
+
         if len(args) not in (0, 1):
             raise ValueError, \
                 "Object ColorBasis must be initialized with 0 or 1 arguments"
@@ -318,6 +333,120 @@
 
         return dict([v, k] for k, v in mydict.items())
 
+    @staticmethod
+    def get_color_flow_string(my_color_string, octet_indices):
+        """Return the color_flow_string (i.e., composed only of T's with 2 
+        indices) associated to my_color_string. Take a list of the external leg
+        color octet state indices as an input. Returns only the leading N 
+        contribution!"""
+
+        # Create a new color factor to allow for simplification
+        my_cf = color_algebra.ColorFactor([my_color_string])
+
+        # Add one T per external octet
+        for indices in octet_indices:
+            my_cf[0].append(color_algebra.T(indices[0],
+                                            indices[1],
+                                            indices[2]))
+
+        # Simplify the whole thing
+        my_cf = my_cf.full_simplify()
+
+        # Return the string with the highest N coefficient 
+        # (leading N decomposition), and the value of this coeff
+        max_coeff = max([cs.Nc_power for cs in my_cf])
+
+        res_cs = [cs for cs in my_cf if cs.Nc_power == max_coeff]
+
+        # If more than one string at leading N...
+        if len(res_cs) > 1:
+            raise ColorBasis.ColorBasisError, \
+             "More than one color string with leading N coeff: %s" % str(res_cs)
+
+        res_cs = res_cs[0]
+
+        # If the result string does not contain only T's with two indices
+        for col_obj in res_cs:
+            if col_obj.__class__.__name__ != 'T':
+                raise ColorBasis.ColorBasisError, \
+                  "Color flow decomposition %s contains non T elements" % \
+                                                                    str(res_cs)
+            if len(col_obj) != 2:
+                raise ColorBasis.ColorBasisError, \
+                  "Color flow decomposition %s contains T's w/o 2 indices" % \
+                                                                    str(res_cs)
+
+        return res_cs
+
+    def color_flow_decomposition(self, repr_dict, ninitial):
+        """Returns the color flow decomposition of the current basis, i.e. a 
+        list of dictionaries (one per color basis entry) with keys corresponding
+        to external leg numbers and values tuples containing two color indices
+        ( (0,0) for singlets, (X,0) for triplet, (0,X) for antitriplet and 
+        (X,Y) for octets). Other color representations are not yet supported 
+        here (an error is raised). Needs a dictionary with keys being external
+        leg numbers, and value the corresponding color representation."""
+
+        # Offsets used to introduce fake quark indices for gluons
+        offset1 = 1000
+        offset2 = 2000
+
+        res = []
+
+        for col_basis_entry in sorted(self.keys()):
+
+            res_dict = {}
+            fake_repl = []
+
+            # Rebuild a color string from a CB entry
+            col_str = color_algebra.ColorString()
+            col_str.from_immutable(col_basis_entry)
+
+            for (leg_num, leg_repr) in repr_dict.items():
+                # By default, assign a (0,0) color flow
+                res_dict[leg_num] = [0, 0]
+
+                # Raise an error if external legs contain non supported repr
+                if leg_repr not in [1, 3, -3, 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:
+                    fake_repl.append((leg_num,
+                                      offset1 + leg_num,
+                                      offset2 + leg_num))
+
+            # Get the actual color flow
+            col_str_flow = self.get_color_flow_string(col_str, fake_repl)
+
+            # Offset for color flow
+            offset = 501
+
+            for col_obj in col_str_flow:
+                for i, index in enumerate(col_obj):
+                    if index < offset1:
+                        res_dict[index][i] = offset
+                    elif index > offset1 and index < offset2:
+                        res_dict[index - offset1][0] = offset
+                    elif index > offset2:
+                        res_dict[index - offset2][1] = offset
+                offset = offset + 1
+
+            # Reverse ordering for initial state to stick to the (weird)
+            # les houches convention
+
+            for key in res_dict.keys():
+                if key <= ninitial:
+                    res_dict[key].reverse()
+
+            res.append(res_dict)
+
+        return res
+
+
+
+
 #===============================================================================
 # ColorMatrix
 #===============================================================================
@@ -364,11 +493,10 @@
         to be symmetric."""
 
         canonical_dict = {}
-
         for i1, struct1 in \
-                    enumerate(self._col_basis1.keys()):
+                    enumerate(sorted(self._col_basis1.keys())):
             for i2, struct2 in \
-                    enumerate(self._col_basis2.keys()):
+                    enumerate(sorted(self._col_basis2.keys())):
 
                 # Only scan upper right triangle if symmetric
                 if is_symmetric and i2 < i1:
@@ -528,3 +656,4 @@
     def lcmm(*args):
         """Return lcm of args."""
         return reduce(ColorMatrix.lcm, args)
+

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2010-02-12 12:53:59 +0000
+++ madgraph/core/diagram_generation.py	2010-05-07 07:20:46 +0000
@@ -83,6 +83,25 @@
 
         return ['process', 'diagrams']
 
+    def get_number_of_diagrams(self):
+        """Returns number of diagrams for this amplitude"""
+        return len(self.get('diagrams'))
+
+    def get_amplitudes(self):
+        """Return an AmplitudeList with just this amplitude.
+        Needed for DecayChainAmplitude."""
+
+        return AmplitudeList([self])
+
+    def nice_string(self, indent=0):
+        """Returns a nicely formatted string of the amplitude content."""
+        return self.get('process').nice_string(indent) + "\n" + \
+               self.get('diagrams').nice_string(indent)
+
+    def nice_string_processes(self, indent=0):
+        """Returns a nicely formatted string of the amplitude process."""
+        return self.get('process').nice_string(indent)
+
     def generate_diagrams(self):
         """Generate diagrams. Algorithm:
 
@@ -116,6 +135,11 @@
 
         Be aware that the resulting vertices have all particles outgoing,
         so need to flip for incoming particles when used.
+
+        SPECIAL CASE: For A>BC... processes which are legs in decay
+        chains, we need to ensure that BC... combine first, giving A=A
+        as a final vertex. This case is defined by the Process
+        property is_decay_chain = True.
         """
 
         model = self['process'].get('model')
@@ -164,7 +188,7 @@
 
             # Need to flip part-antipart for incoming particles, 
             # so they are all outgoing
-            if leg.get('state') == 'initial':
+            if leg.get('state') == False:
                 part = model.get('particle_dict')[leg.get('id')]
                 leg.set('id', part.get_anti_pdg_code())
 
@@ -177,9 +201,30 @@
         # Reduce the leg list and return the corresponding
         # list of vertices
 
-        reduced_leglist = self.reduce_leglist(leglist,
-                                              max_multi_to1,
-                                              self.get('process').get('orders'))
+        # Set is_decay_chain to True if this is a 1->N decay process
+        # as part of a decay chain.
+        is_decay_chain = self['process'].get('is_decay_chain')
+        if is_decay_chain:
+            part = model.get('particle_dict')[leglist[0].get('id')]
+            # For decay chain legs, we want everything to combine to
+            # the initial leg. This is done by only allowing the
+            # initial leg to combine as a final identity.
+            ref_dict_to0 = {(part.get_pdg_code(),part.get_anti_pdg_code()):0,
+                            (part.get_anti_pdg_code(),part.get_pdg_code()):0}
+            # Need to set initial leg from_group to None, to make sure
+            # it can only be combined at the end.
+            leglist[0].set('from_group', None)
+            reduced_leglist = self.reduce_leglist(leglist,
+                                                  max_multi_to1,
+                                                  ref_dict_to0,
+                                                  is_decay_chain,
+                                                  self.get('process').get('orders'))
+        else:
+            reduced_leglist = self.reduce_leglist(leglist,
+                                                  max_multi_to1,
+                                                  model.get('ref_dict_to0'),
+                                                  is_decay_chain,
+                                                  self.get('process').get('orders'))
 
         for vertex_list in reduced_leglist:
             res.append(base_objects.Diagram(
@@ -194,7 +239,7 @@
         # Note that we shouldn't look at the last vertex in each
         # diagram, since that is the n->0 vertex
         if self['process'].get('required_s_channels'):
-            ninitial = len(filter(lambda leg: leg.get('state') == 'initial',
+            ninitial = len(filter(lambda leg: leg.get('state') == False,
                               self['process'].get('legs')))
             res = base_objects.DiagramList(\
                 filter(lambda diagram: \
@@ -210,7 +255,7 @@
         # Note that we shouldn't look at the last vertex in each
         # diagram, since that is the n->0 vertex
         if self['process'].get('forbidden_s_channels'):
-            ninitial = len(filter(lambda leg: leg.get('state') == 'initial',
+            ninitial = len(filter(lambda leg: leg.get('state') == False,
                               self['process'].get('legs')))
             res = base_objects.DiagramList(\
                 filter(lambda diagram: \
@@ -223,13 +268,16 @@
         # Set diagrams to res
         self['diagrams'] = res
 
+        # Trim down number of legs and vertices used to save memory
+        self.trim_diagrams()
+
         if res:
             logger.info("Process has %d diagrams" % len(res))
 
         return not failed_crossing
 
-    def reduce_leglist(self, curr_leglist, max_multi_to1,
-                       coupling_orders=None):
+    def reduce_leglist(self, curr_leglist, max_multi_to1, ref_dict_to0,
+                       is_decay_chain = False, coupling_orders = None):
         """Recursive function to reduce N LegList to N-1
            For algorithm, see doc for generate_diagrams.
         """
@@ -245,13 +293,13 @@
 
         # Extract ref dict information
         model = self['process'].get('model')
-        ref_dict_to0 = self['process'].get('model').get('ref_dict_to0')
         ref_dict_to1 = self['process'].get('model').get('ref_dict_to1')
 
 
         # If all legs can be combined in one single vertex, add this
-        # vertex to res and continue
-        if curr_leglist.can_combine_to_0(ref_dict_to0):
+        # vertex to res and continue.
+        # Special treatment for decay chain legs
+        if curr_leglist.can_combine_to_0(ref_dict_to0, is_decay_chain):
             # Extract the interaction id associated to the vertex 
             vertex_id = ref_dict_to0[tuple(sorted([leg.get('id') for \
                                                    leg in curr_leglist]))]
@@ -300,6 +348,8 @@
             # First, reduce again the leg part
             reduced_diagram = self.reduce_leglist(leg_vertex_tuple[0],
                                                   max_multi_to1,
+                                                  ref_dict_to0,
+                                                  is_decay_chain,
                                                   new_coupling_orders)
             # If there is a reduced diagram
             if reduced_diagram:
@@ -444,11 +494,11 @@
                     number = min([leg.get('number') for leg in entry])
                     # 3) state is final, unless there is exactly one initial 
                     # state particle involved in the combination -> t-channel
-                    if len(filter(lambda leg: leg.get('state') == 'initial',
+                    if len(filter(lambda leg: leg.get('state') == False,
                                   entry)) == 1:
-                        state = 'initial'
+                        state = False
                     else:
-                        state = 'final'
+                        state = True
                     # 4) from_group is True, by definition
 
                     # Create and add the object
@@ -484,7 +534,11 @@
                 # and add it
                 else:
                     cp_entry = copy.copy(entry)
-                    cp_entry.set('from_group', False)
+                    # Need special case for from_group == None; this
+                    # is for initial state leg of decay chain process
+                    # (see Leg.can_combine_to_0)
+                    if cp_entry.get('from_group') != None:
+                        cp_entry.set('from_group', False)
                     reduced_list.append(cp_entry)
 
             # Flatten the obtained leg and vertex lists
@@ -498,6 +552,28 @@
 
         return res
 
+    def trim_diagrams(self):
+        """Reduce the number of legs and vertices used in memory."""
+
+        legs = []
+        vertices = []
+
+        for diagram in self.get('diagrams'):
+            for ivx, vertex in enumerate(diagram.get('vertices')):
+                for ileg, leg in enumerate(vertex.get('legs')):
+                    leg.set('from_group', False)
+                    try:
+                        index = legs.index(leg)
+                        vertex.get('legs')[ileg] = legs[index]
+                    except ValueError:
+                        legs.append(leg)
+                try:
+                    index = vertices.index(vertex)
+                    diagram.get('vertices')[ivx] = vertices[index]
+                except ValueError:
+                    vertices.append(vertex)
+        
+
 #===============================================================================
 # AmplitudeList
 #===============================================================================
@@ -511,6 +587,143 @@
         return isinstance(obj, Amplitude)
 
 #===============================================================================
+# DecayChainAmplitude
+#===============================================================================
+class DecayChainAmplitude(Amplitude):
+    """A list of amplitudes + a list of decay chain amplitude lists;
+    corresponding to a ProcessDefinition with a list of decay chains
+    """
+
+    def default_setup(self):
+        """Default values for all properties"""
+
+        self['amplitudes'] = AmplitudeList()
+        self['decay_chains'] = DecayChainAmplitudeList()
+
+    def __init__(self, argument=None):
+        """Allow initialization with Process and with ProcessDefinition"""
+
+        if isinstance(argument, base_objects.Process):
+            super(DecayChainAmplitude, self).__init__()
+            if isinstance(argument, base_objects.ProcessDefinition):
+                self['amplitudes'].extend(\
+                MultiProcess.generate_multi_amplitudes(argument))
+            else:
+                self['amplitudes'].append(Amplitude(argument))
+                # Clean decay chains from process, since we haven't
+                # combined processes with decay chains yet
+                process = copy.copy(self['amplitudes'][0].get('process'))
+                process.set('decay_chains', base_objects.ProcessList())
+                self['amplitudes'][0].set('process', process)
+            for process in argument.get('decay_chains'):
+                if not process.get('is_decay_chain'):
+                    process.set('is_decay_chain',True)
+                if not process.get_ninitial() == 1:
+                    raise self.PhysicsObjectError,\
+                          "Decay chain process must have exactly one" + \
+                          " incoming particle"
+                self['decay_chains'].append(\
+                    DecayChainAmplitude(process))
+        elif argument != None:
+            # call the mother routine
+            super(DecayChainAmplitude, self).__init__(argument)
+        else:
+            # call the mother routine
+            super(DecayChainAmplitude, self).__init__()
+
+    def filter(self, name, value):
+        """Filter for valid amplitude property values."""
+
+        if name == 'amplitudes':
+            if not isinstance(value, AmplitudeList):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid AmplitudeList" % str(value)
+        if name == 'decay_chains':
+            if not isinstance(value, DecayChainAmplitudeList):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid DecayChainAmplitudeList object" % \
+                        str(value)
+        return True
+
+    def get_sorted_keys(self):
+        """Return diagram property names as a nicely sorted list."""
+
+        return ['amplitudes', 'decay_chains']
+
+    # Helper functions
+
+    def get_number_of_diagrams(self):
+        """Returns number of diagrams for this amplitude"""
+        return sum(len(a.get('diagrams')) for a in self.get('amplitudes')) \
+               + sum(d.get_number_of_diagrams() for d in \
+                                        self.get('decay_chains'))
+
+    def nice_string(self, indent = 0):
+        """Returns a nicely formatted string of the amplitude content."""
+        mystr = ""
+        for amplitude in self.get('amplitudes'):
+            mystr = mystr + amplitude.nice_string(indent) + "\n"
+
+        if self.get('decay_chains'):
+            mystr = mystr + " " * indent + "Decays:\n"
+        for dec in self.get('decay_chains'):
+            mystr = mystr + dec.nice_string(indent + 2) + "\n"
+
+        return  mystr[:-1]
+
+    def nice_string_processes(self, indent = 0):
+        """Returns a nicely formatted string of the amplitude processes."""
+        mystr = ""
+        for amplitude in self.get('amplitudes'):
+            mystr = mystr + amplitude.nice_string_processes(indent) + "\n"
+
+        if self.get('decay_chains'):
+            mystr = mystr + " " * indent + "Decays:\n"
+        for dec in self.get('decay_chains'):
+            mystr = mystr + dec.nice_string_processes(indent + 2) + "\n"
+
+        return  mystr[:-1]
+
+    def get_decay_ids(self):
+        """Returns a set of all particle ids for which a decay is defined"""
+
+        decay_ids = []
+
+        # Get all amplitudes for the decay processes
+        for amp in sum([dc.get('amplitudes') for dc \
+                        in self['decay_chains']], []):
+            # For each amplitude, find the initial state leg
+            decay_ids.append(amp.get('process').get_initial_ids()[0])
+            
+        # Return a list with unique ids
+        return list(set(decay_ids))
+    
+    def get_amplitudes(self):
+        """Recursive function to extract all amplitudes for this process"""
+
+        amplitudes = AmplitudeList()
+
+        amplitudes.extend(self.get('amplitudes'))
+        for decay in self.get('decay_chains'):
+            amplitudes.extend(decay.get_amplitudes())
+
+        return amplitudes
+            
+
+#===============================================================================
+# DecayChainAmplitudeList
+#===============================================================================
+class DecayChainAmplitudeList(base_objects.PhysicsObjectList):
+    """List of DecayChainAmplitude objects
+    """
+
+    def is_valid_element(self, obj):
+        """Test if object obj is a valid DecayChainAmplitude for the list."""
+
+        return isinstance(obj, DecayChainAmplitude)
+
+    
+#===============================================================================
 # MultiProcess
 #===============================================================================
 class MultiProcess(base_objects.PhysicsObject):
@@ -523,8 +736,30 @@
         """Default values for all properties"""
 
         self['process_definitions'] = base_objects.ProcessDefinitionList()
+        # self['amplitudes'] can be an AmplitudeList or a
+        # DecayChainAmplitudeList, depending on whether there are
+        # decay chains in the process definitions or not.
         self['amplitudes'] = AmplitudeList()
 
+    def __init__(self, argument=None):
+        """Allow initialization with ProcessDefinition or
+        ProcessDefinitionList"""
+
+        if isinstance(argument, base_objects.ProcessDefinition):
+            super(MultiProcess, self).__init__()
+            self['process_definitions'].append(argument)
+            self.get('amplitudes')
+        elif isinstance(argument, base_objects.ProcessDefinitionList):
+            super(MultiProcess, self).__init__()
+            self['process_definitions'] = argument
+            self.get('amplitudes')
+        elif argument != None:
+            # call the mother routine
+            super(MultiProcess, self).__init__(argument)
+        else:
+            # call the mother routine
+            super(MultiProcess, self).__init__()
+
     def filter(self, name, value):
         """Filter for valid process property values."""
 
@@ -544,8 +779,15 @@
         """Get the value of the property name."""
 
         if (name == 'amplitudes') and not self[name]:
-            if self.get('process_definitions'):
-                self.generate_amplitudes()
+            for process_def in self.get('process_definitions'):
+                if process_def.get('decay_chains'):
+                    # This is a decay chain process
+                    # Store amplitude(s) as DecayChainAmplitude
+                    self['amplitudes'].append(\
+                        DecayChainAmplitude(process_def))
+                else:
+                    self['amplitudes'].extend(\
+                        MultiProcess.generate_multi_amplitudes(process_def))
 
         return MultiProcess.__bases__[0].get(self, name) # call the mother routine
 
@@ -554,7 +796,8 @@
 
         return ['process_definitions', 'amplitudes']
 
-    def generate_amplitudes(self):
+    @staticmethod
+    def generate_multi_amplitudes(process_definition):
         """Generate amplitudes in a semi-efficient way.
         Make use of crossing symmetry for processes that fail diagram
         generation, but not for processes that succeed diagram
@@ -562,79 +805,86 @@
         identify processes with identical amplitudes.
         """
 
-        # Loop over all process definitions and multilegs
+        if not isinstance(process_definition, base_objects.ProcessDefinition):
+            raise base_objects.PhysicsObjectError,\
+                  "%s not valid ProcessDefinition object" % \
+                  repr(process_definition)
+
         processes = base_objects.ProcessList()
-
-        for process_def in self['process_definitions']:
-
-            # failed_procs are processes that have already failed
-            # based on crossing symmetry
-            failed_procs = []
-                    
-            model = process_def['model']
-
-            isids = [leg['ids'] for leg in \
-                     filter(lambda leg: leg['state'] == 'initial', process_def['legs'])]
-            fsids = [leg['ids'] for leg in \
-                     filter(lambda leg: leg['state'] == 'final', process_def['legs'])]
-
-            # Generate all combinations for the initial state
-
-            for prod in apply(itertools.product, isids):
-                islegs = [\
-                    base_objects.Leg({'id':id, 'state': 'initial'}) \
+        amplitudes = AmplitudeList()
+
+        # failed_procs are processes that have already failed
+        # based on crossing symmetry
+        failed_procs = []
+        
+        model = process_definition['model']
+        
+        isids = [leg['ids'] for leg in \
+                 filter(lambda leg: leg['state'] == False, process_definition['legs'])]
+        fsids = [leg['ids'] for leg in \
+                 filter(lambda leg: leg['state'] == True, process_definition['legs'])]
+
+        # Generate all combinations for the initial state
+        
+        for prod in apply(itertools.product, isids):
+            islegs = [\
+                    base_objects.Leg({'id':id, 'state': False}) \
                     for id in prod]
 
-                # Generate all combinations for the final state, and make
-                # sure to remove double counting
-
-                red_fsidlist = []
-
-                for prod in apply(itertools.product, fsids):
-
-                    # Remove double counting between final states
-                    if tuple(sorted(prod)) in red_fsidlist:
-                        continue
-
-                    red_fsidlist.append(tuple(sorted(prod)));
-
-                    # Generate leg list for process
-                    leg_list = [copy.copy(leg) for leg in islegs]
-
-                    leg_list.extend([\
-                        base_objects.Leg({'id':id, 'state': 'final'}) \
+            # Generate all combinations for the final state, and make
+            # sure to remove double counting
+
+            red_fsidlist = []
+
+            for prod in apply(itertools.product, fsids):
+
+                # Remove double counting between final states
+                if tuple(sorted(prod)) in red_fsidlist:
+                    continue
+                
+                red_fsidlist.append(tuple(sorted(prod)));
+                
+                # Generate leg list for process
+                leg_list = [copy.copy(leg) for leg in islegs]
+                
+                leg_list.extend([\
+                        base_objects.Leg({'id':id, 'state': True}) \
                         for id in prod])
-
-                    legs = base_objects.LegList(leg_list)
-
-                    # Setup process
-                    process = base_objects.Process({\
-                                       'legs':legs,
-                                       'model':process_def.get('model'),
-                                       'id': process_def.get('id'),
-                                       'orders': process_def.get('orders'),
-                                       'required_s_channels': \
-                                       process_def.get('required_s_channels'),
-                                       'forbidden_s_channels': \
-                                       process_def.get('forbidden_s_channels'),
-                                       'forbidden_particles': \
-                                       process_def.get('forbidden_particles')})
-
-                    # Check for crossed processes
-                    sorted_legs = sorted(legs.get_outgoing_id_list(model))
-                    # Check if crossed process has already failed
-                    # In that case don't check process
-                    if tuple(sorted_legs) in failed_procs:
-                        continue
-
-                    amplitude = Amplitude({"process": process})
-                    if not amplitude.generate_diagrams():
-                        # Add process to failed_procs
-                        failed_procs.append(tuple(sorted_legs))
-                    if amplitude.get('diagrams'):
-                        self['amplitudes'].append(amplitude)
-                    
-
+                
+                legs = base_objects.LegList(leg_list)
+
+                # Setup process
+                process = base_objects.Process({\
+                              'legs':legs,
+                              'model':process_definition.get('model'),
+                              'id': process_definition.get('id'),
+                              'orders': process_definition.get('orders'),
+                              'required_s_channels': \
+                                 process_definition.get('required_s_channels'),
+                              'forbidden_s_channels': \
+                                 process_definition.get('forbidden_s_channels'),
+                              'forbidden_particles': \
+                                 process_definition.get('forbidden_particles'),
+                              'is_decay_chain': \
+                                 process_definition.get('is_decay_chain')})
+
+                # Check for crossed processes
+                sorted_legs = sorted(legs.get_outgoing_id_list(model))
+                # Check if crossed process has already failed
+                # In that case don't check process
+                if tuple(sorted_legs) in failed_procs:
+                    continue
+
+                amplitude = Amplitude({"process": process})
+                if not amplitude.generate_diagrams():
+                    # Add process to failed_procs
+                    failed_procs.append(tuple(sorted_legs))
+                if amplitude.get('diagrams'):
+                    amplitudes.append(amplitude)
+
+        # Return the produced amplitudes
+        return amplitudes
+            
 #===============================================================================
 # Global helper methods
 #===============================================================================

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2010-02-13 04:20:50 +0000
+++ madgraph/core/helas_objects.py	2010-05-07 07:20:46 +0000
@@ -15,6 +15,7 @@
 
 import array
 import copy
+import gc
 import logging
 import re
 import itertools
@@ -64,15 +65,9 @@
         # width = 'zero'
         # is_part = 'true'    Particle not antiparticle
         # self_antipart='false'   gluon, photo, h, or majorana would be true
-        self['pdg_code'] = 0
-        self['name'] = 'none'
-        self['antiname'] = 'none'
-        self['spin'] = 1
-        self['color'] = 1
-        self['mass'] = 'zero'
-        self['width'] = 'zero'
+        self['particle'] = base_objects.Particle()
+        self['antiparticle'] = base_objects.Particle()
         self['is_part'] = True
-        self['self_antipart'] = False
         # Properties related to the interaction generating the propagator
         # For an e- produced from an e+e-A vertex would have the following
         # proporties:
@@ -92,6 +87,8 @@
         # state = initial/final (for external bosons),
         #         intermediate (for intermediate bosons),
         #         incoming/outgoing (for fermions)
+        # leg_state = initial/final for initial/final legs
+        #             intermediate for non-external wavefunctions
         # number_external = the 'number' property of the corresponding Leg,
         #                   corresponds to the number of the first external
         #                   particle contributing to this leg
@@ -99,10 +96,18 @@
         #                    -1 is used only if there is a fermion flow clash
         #                    due to a Majorana particle 
         self['state'] = 'incoming'
+        self['leg_state'] = True
         self['mothers'] = HelasWavefunctionList()
         self['number_external'] = 0
         self['number'] = 0
         self['fermionflow'] = 1
+        # The decay flag is used in processes with defined decay chains,
+        # to indicate that this wavefunction has a decay defined
+        self['decay'] = False
+        # The onshell flag is used in processes with defined decay
+        # chains, to indicate that this wavefunction is decayed and
+        # should be onshell
+        self['onshell'] = False
 
     # Customized constructor
     def __init__(self, *arguments):
@@ -117,17 +122,30 @@
                 leg = arguments[0]
                 interaction_id = arguments[1]
                 model = arguments[2]
-                self.set('pdg_code', leg.get('id'), model)
+                # decay_ids is the pdg codes for particles with decay
+                # chains defined
+                decay_ids = []
+                if len(arguments) > 3:
+                    decay_ids = arguments[3]
+                self.set('particle', leg.get('id'), model)
                 self.set('number_external', leg.get('number'))
                 self.set('number', leg.get('number'))
-                self.set('state', leg.get('state'))
+                self.set('state', {False: 'initial', True: 'final'}[leg.get('state')])
+                self.set('leg_state', leg.get('state'))
+                # Need to set 'decay' to True for particles which will be
+                # decayed later, in order to not combine such processes
+                # although they might have identical matrix elements before
+                # the decay is applied
+                if self['state'] == 'final' and self.get('pdg_code') in decay_ids:
+                    self.set('decay', True)
+
                 # Set fermion flow state. Initial particle and final
                 # antiparticle are incoming, and vice versa for
                 # outgoing
                 if self.is_fermion():
-                    if leg.get('state') == 'initial' and \
+                    if leg.get('state') == False and \
                            self.get('is_part') or \
-                           leg.get('state') == 'final' and \
+                           leg.get('state') == True and \
                            not self.get('is_part'):
                         self.set('state', 'incoming')
                     else:
@@ -141,45 +159,12 @@
     def filter(self, name, value):
         """Filter for valid wavefunction property values."""
 
-        if name == 'pdg_code':
-            if not isinstance(value, int):
-                raise self.PhysicsObjectError, \
-                      "%s is not a valid pdg_code for wavefunction" % \
-                      str(value)
-
-        if name in ['name', 'antiname']:
-            # Must start with a letter, followed by letters,  digits,
-            # - and + only
-            p = re.compile('\A[a-zA-Z]+[\w]*[\-\+]*~?\Z')
-            if not p.match(value):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid particle name" % value
-
-        if name is 'spin':
-            if not isinstance(value, int):
-                raise self.PhysicsObjectError, \
-                    "Spin %s is not an integer" % repr(value)
-            if value < 1 or value > 5:
-                raise self.PhysicsObjectError, \
-                   "Spin %i is smaller than one" % value
-
-        if name is 'color':
-            if not isinstance(value, int):
-                raise self.PhysicsObjectError, \
-                    "Color %s is not an integer" % repr(value)
-            if value not in [1, 3, 6, 8]:
-                raise self.PhysicsObjectError, \
-                   "Color %i is not valid" % value
-
-        if name in ['mass', 'width']:
-            # Must start with a letter, followed by letters, digits or _
-            p = re.compile('\A[a-zA-Z]+[\w\_]*\Z')
-            if not p.match(value):
-                raise self.PhysicsObjectError, \
-                        "%s is not a valid name for mass/width variable" % \
-                        value
-
-        if name in ['is_part', 'self_antipart']:
+        if name in ['particle', 'antiparticle']:
+            if not isinstance(value, base_objects.Particle):
+                raise self.PhysicsObjectError, \
+                    "%s tag %s is not a particle" % (name, repr(value))            
+
+        if name == 'is_part':
             if not isinstance(value, bool):
                 raise self.PhysicsObjectError, \
                     "%s tag %s is not a boolean" % (name, repr(value))
@@ -220,7 +205,7 @@
                                                                 str(value)
 
         if name == 'coupl_key':
-            if not isinstance(value, tuple):
+            if value and not isinstance(value, tuple):
                 raise self.PhysicsObjectError, \
                       "%s is not a valid tuple" % str(value)
             if len(value) != 2:
@@ -240,6 +225,11 @@
                 raise self.PhysicsObjectError, \
                         "%s is not a valid wavefunction " % str(value) + \
                         "state (incoming|outgoing|intermediate)"
+        if name == 'leg_state':
+            if value not in [False, True]:
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid wavefunction " % str(value) + \
+                        "state (incoming|outgoing|intermediate)"
         if name in ['fermionflow']:
             if not isinstance(value, int):
                 raise self.PhysicsObjectError, \
@@ -260,8 +250,33 @@
                       "%s is not a valid list of mothers for wavefunction" % \
                       str(value)
 
+        if name in ['decay', 'onshell']:
+            if not isinstance(value, bool):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid bool" % str(value) + \
+                        " for decay or onshell"
+
         return True
 
+    # Enhanced get function, where we can directly call the properties of the particle
+    def get(self, name):
+        """When calling any property related to the particle,
+        automatically call the corresponding property of the particle."""
+
+        if name in ['spin', 'mass', 'width', 'self_antipart']:
+            return self['particle'].get(name)
+        elif name == 'pdg_code':
+            return self['particle'].get_pdg_code()
+        elif name == 'color':
+            return self['particle'].get_color()
+        elif name == 'name':
+            return self['particle'].get_name()
+        elif name == 'antiname':
+            return self['particle'].get_anti_name()
+        else:
+            return super(HelasWavefunction, self).get(name)
+        
+
     # Enhanced set function, where we can append a model
 
     def set(self, *arguments):
@@ -296,17 +311,13 @@
                     if inter.get('couplings'):
                         self.set('coupling', inter.get('couplings').values()[0])
                 return True
-            elif name == 'pdg_code':
-                self.set('pdg_code', value)
-                part = model.get('particle_dict')[value]
-                self.set('name', part.get('name'))
-                self.set('antiname', part.get('antiname'))
-                self.set('spin', part.get('spin'))
-                self.set('color', part.get('color'))
-                self.set('mass', part.get('mass'))
-                self.set('width', part.get('width'))
-                self.set('is_part', part.get('is_part'))
-                self.set('self_antipart', part.get('self_antipart'))
+            elif name == 'particle':
+                self.set('particle', model.get('particle_dict')[value])
+                self.set('is_part', self['particle'].get('is_part'))
+                if self['particle'].get('self_antipart'):
+                    self.set('antiparticle', self['particle'])
+                else:
+                    self.set('antiparticle', model.get('particle_dict')[-value])
                 return True
             else:
                 raise self.PhysicsObjectError, \
@@ -317,14 +328,19 @@
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
 
-        return ['pdg_code', 'name', 'antiname', 'spin', 'color',
-                'mass', 'width', 'is_part', 'self_antipart',
+        return ['particle', 'antiparticle', 'is_part',
                 'interaction_id', 'pdg_codes', 'inter_color', 'lorentz',
                 'coupling', 'coupl_key', 'state', 'number_external',
                 'number', 'fermionflow', 'mothers']
 
     # Helper functions
 
+    def flip_part_antipart(self):
+        """Flip between particle and antiparticle."""
+        part = self.get('particle')
+        self.set('particle', self.get('antiparticle'))
+        self.set('antiparticle', part)
+
     def is_fermion(self):
         return self.get('spin') % 2 == 0
 
@@ -343,29 +359,17 @@
         # Lorentz or color structures
         array_rep.extend(list(self['coupl_key']))
         # Finally, the mother numbers
-        array_rep.extend([mother.get('number') for \
+        array_rep.extend([mother['number'] for \
                           mother in self['mothers']])
         return array_rep
 
-    def get_pdg_code_outgoing(self):
+    def get_pdg_code(self):
         """Generate the corresponding pdg_code for an outgoing particle,
         taking into account fermion flow, for mother wavefunctions"""
 
-        if self.get('self_antipart'):
-            #This is its own antiparticle e.g. a gluon
-            return self.get('pdg_code')
-
-        if self.is_boson():
-            # This is a boson
-            return self.get('pdg_code')
-
-        if (self.get('state') == 'incoming' and self.get('is_part') \
-                or self.get('state') == 'outgoing' and not self.get('is_part')):
-            return - self.get('pdg_code')
-        else:
-            return self.get('pdg_code')
-
-    def get_pdg_code_incoming(self):
+        return self.get('pdg_code')
+
+    def get_anti_pdg_code(self):
         """Generate the corresponding pdg_code for an incoming particle,
         taking into account fermion flow, for mother wavefunctions"""
 
@@ -373,15 +377,7 @@
             #This is its own antiparticle e.g. gluon
             return self.get('pdg_code')
 
-        if self.is_boson():
-            # This is a boson
-            return - self.get('pdg_code')
-
-        if (self.get('state') == 'outgoing' and self.get('is_part') \
-                or self.get('state') == 'incoming' and not self.get('is_part')):
-            return - self.get('pdg_code')
-        else:
-            return self.get('pdg_code')
+        return - self.get('pdg_code')
 
     def set_scalar_coupling_sign(self, model):
         """Check if we need to add a minus sign due to non-identical
@@ -408,6 +404,16 @@
             raise self.PhysicsObjectError, \
                   "%s is not a valid model for call to set_state_and_particle" \
                   % repr(model)
+
+        # leg_state is final, unless there is exactly one initial 
+        # state particle involved in the combination -> t-channel
+        if len(filter(lambda mother: mother.get('leg_state') == False,
+                      self.get('mothers'))) == 1:
+            leg_state = False
+        else:
+            leg_state = True
+        self.set('leg_state', leg_state)
+
         # Start by setting the state of the wavefunction
         if self.is_boson():
             # For boson, set state to intermediate
@@ -449,14 +455,6 @@
                 self.set('state', mother.get_with_flow('state'))
                 self.set('is_part', mother.get_with_flow('is_part'))
 
-        # We want the particle created here to go into the next
-        # vertex, so we need to flip identity for incoming
-        # antiparticle and outgoing particle.
-        if not self.get('self_antipart') and \
-               (self.get('state') == 'incoming' and not self.get('is_part') \
-                or self.get('state') == 'outgoing' and self.get('is_part')):
-            self.set('pdg_code', -self.get('pdg_code'), model)
-
         return True
 
     def check_and_fix_fermion_flow(self,
@@ -497,13 +495,14 @@
                                    external_wavefunctions,
                                    self.get_with_flow('state'),
                                    wf_number)
+
         return self, wf_number
 
     def check_majorana_and_flip_flow(self, found_majorana,
                                      wavefunctions,
                                      diagram_wavefunctions,
                                      external_wavefunctions,
-                                     wf_number):
+                                     wf_number, force_flip_flow=False):
         """Recursive function. Check for Majorana fermion. If found,
         continue down to external leg, then flip all the fermion flows
         on the way back up, in the correct way:
@@ -511,6 +510,15 @@
         wavefunctions before the last Majorana fermion, instead flip
         particle identities and state. Return the new (or old)
         wavefunction, and the present wavefunction number.
+
+        Arguments:
+          found_majorana: boolean
+          wavefunctions: HelasWavefunctionList with previously
+                         defined wavefunctions
+          diagram_wavefunctions: HelasWavefunctionList with the wavefunctions
+                         already defined in this diagram
+          external_wavefunctions: dictionary from legnumber to external wf
+          wf_number: The present wavefunction number
         """
 
         if not found_majorana:
@@ -523,7 +531,9 @@
         # Stop recursion at the external leg
         mothers = copy.copy(self.get('mothers'))
         if not mothers:
-            if not self.get('self_antipart'):
+            if force_flip_flow:
+                flip_flow = True
+            elif not self.get('self_antipart'):
                 flip_flow = found_majorana
             else:
                 flip_sign = found_majorana
@@ -538,17 +548,20 @@
                 raise self.PhysicsObjectError, \
                       "6-fermion vertices not yet implemented"
             if len(fermion_mother) == 0:
-                raise self.PhysicsObjectError, \
-                      "Previous unresolved fermion flow in mother chain"
-
-            # Perform recursion by calling on mother
-            new_mother, wf_number = fermion_mother[0].\
-                                    check_majorana_and_flip_flow(\
-                                         found_majorana,
-                                         wavefunctions,
-                                         diagram_wavefunctions,
-                                         external_wavefunctions,
-                                         wf_number)
+                # Fermion flow already resolved for mothers
+                fermion_mother = filter(lambda wf: wf.is_fermion(),
+                                    mothers)
+                new_mother = fermion_mother[0]
+            else:
+                # Perform recursion by calling on mother
+                new_mother, wf_number = fermion_mother[0].\
+                                        check_majorana_and_flip_flow(\
+                                           found_majorana,
+                                           wavefunctions,
+                                           diagram_wavefunctions,
+                                           external_wavefunctions,
+                                           wf_number,
+                                           force_flip_flow)
 
             # If this is Majorana and mother has different fermion
             # flow, it means that we should from now on in the chain
@@ -575,8 +588,29 @@
                 # Update wavefunction number
                 wf_number = wf_number + 1
                 new_wf.set('number', wf_number)
-                new_wf.set('mothers', mothers)
-                diagram_wavefunctions.append(new_wf)
+                try:
+                    # In call from insert_decay, we want to replace
+                    # also identical wavefunctions in the same diagram
+                    old_wf_index = diagram_wavefunctions.index(self)
+                    diagram_wavefunctions[old_wf_index] = new_wf
+                except ValueError:
+                    diagram_wavefunctions.append(new_wf)
+                    # Make sure that new_wf comes before any wavefunction
+                    # which has it as mother
+                    for i, wf in enumerate(diagram_wavefunctions):
+                        if self in wf.get('mothers'):
+                            # Remove new_wf, in order to insert it below
+                            diagram_wavefunctions.pop()
+                            # Update wf numbers
+                            new_wf.set('number', wf.get('number'))
+                            for w in diagram_wavefunctions[i:]:
+                                w.set('number', w.get('number') + 1)
+                            # Insert wavefunction
+                            diagram_wavefunctions.insert(i, new_wf)
+                            break
+
+            # Set new mothers
+            new_wf.set('mothers', mothers)
 
             # Now flip flow or sign
             if flip_flow:
@@ -590,14 +624,17 @@
                                            state != new_wf.get('state'),
                                            ['incoming', 'outgoing'])[0])
                 new_wf.set('is_part', not new_wf.get('is_part'))
-                if not new_wf.get('self_antipart'):
-                    new_wf.set('pdg_code', -new_wf.get('pdg_code'))
 
             try:
                 # Use the copy in wavefunctions instead.
                 # Remove this copy from diagram_wavefunctions
                 new_wf = wavefunctions[wavefunctions.index(new_wf)]
-                diagram_wavefunctions.remove(new_wf)
+                index = diagram_wavefunctions.index(new_wf)
+                diagram_wavefunctions.pop(index)
+                # We need to decrease the wf number for later
+                # diagram wavefunctions
+                for wf in diagram_wavefunctions[index:]:
+                    wf.set('number', wf.get('number') - 1)
                 # Since we reuse the old wavefunction, reset wf_number
                 wf_number = wf_number - 1
             except ValueError:
@@ -659,7 +696,7 @@
             fermion_number_list.extend(in_list[1])
             fermion_number_list.extend(out_list[1])
         elif len(in_fermions) != len(out_fermions):
-            raise self.HelasWavefunctionError, \
+            raise self.PhysicsObjectError, \
                   "Error: %d incoming fermions != %d outgoing fermions" % \
                   (len(in_fermions), len(out_fermions))
 
@@ -718,6 +755,211 @@
 
         return (tuple(res), self.get('lorentz'))
 
+    def get_base_vertices(self, wf_dict = {}, vx_list = [], optimization = 1):
+        """Recursive method to get a base_objects.VertexList
+        corresponding to this wavefunction and its mothers."""
+
+        vertices = base_objects.VertexList()
+
+        if not self.get('mothers'):
+            return vertices
+
+        # Add vertices for all mothers
+        for mother in self.get('mothers'):
+            # This is where recursion happens
+            vertices.extend(mother.get_base_vertices(wf_dict, vx_list,
+                                                     optimization))
+        # Generate last vertex
+        legs = base_objects.LegList()
+
+        # We use the from_group flag to indicate whether this outgoing
+        # leg corresponds to a decaying (onshell) particle or not
+        try:
+            lastleg = wf_dict[self.get('number')]
+        except KeyError:            
+            lastleg = base_objects.Leg({
+                'id': self.get_pdg_code(),
+                'number': self.get('number_external'),
+                'state': self.get('leg_state'),
+                'from_group': self.get('onshell')
+                })
+            if optimization != 0:
+                wf_dict[self.get('number')] = lastleg
+
+        for mother in self.get('mothers'):
+            try:
+                leg = wf_dict[mother.get('number')]
+            except KeyError:
+                leg = base_objects.Leg({
+                    'id': mother.get_pdg_code(),
+                    'number': mother.get('number_external'),
+                    'state': mother.get('leg_state'),
+                    'from_group': mother.get('onshell')
+                    })
+                if optimization != 0:
+                    wf_dict[mother.get('number')] = leg
+            legs.append(leg)
+
+        legs.append(lastleg)
+
+        vertex = base_objects.Vertex({
+            'id': self.get('interaction_id'),
+            'legs': legs})
+
+        try:
+            index = vx_list.index(vertex)
+            vertex = vx_list[index]
+        except ValueError:
+            pass
+        
+        vertices.append(vertex)
+
+        return vertices
+
+    def get_color_indices(self):
+        """Recursive method to get the color indices corresponding to
+        this wavefunction and its mothers."""
+
+        if not self.get('mothers'):
+            return []
+
+        color_indices = []
+
+        # Add color indices for all mothers
+        for mother in self.get('mothers'):
+            # This is where recursion happens
+            color_indices.extend(mother.get_color_indices())
+        # Add this wf's color index
+        color_indices.append(self.get('coupl_key')[0])
+
+        return color_indices
+
+    def get_s_and_t_channels(self, ninitial, mother_leg):
+        """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."""
+
+        schannels = base_objects.VertexList()
+        tchannels = base_objects.VertexList()
+
+        mother_leg = copy.copy(mother_leg)
+
+        # Add vertices for all s-channel mothers
+        final_mothers = filter(lambda wf: wf.get('number_external') > ninitial,
+                               self.get('mothers'))
+
+        for mother in final_mothers:
+            schannels.extend(mother.get_base_vertices(optimization = 0))
+
+        # Extract initial state mothers
+        init_mothers = filter(lambda wf: wf.get('number_external') <= ninitial,
+                              self.get('mothers'))
+
+        if len(init_mothers) > 2:
+            raise self.PhysicsObjectError, \
+                  "get_s_and_t_channels can only handle up to 2 initial states"
+
+        if len(init_mothers) == 1:
+            # This is an s-channel or t-channel leg. 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
+
+            for mother in mothers:
+                legs.append(base_objects.Leg({
+                    'id': mother.get_pdg_code(),
+                    'number': mother.get('number_external'),
+                    'state': mother.get('leg_state'),
+                    'from_group': False
+                    }))
+
+            if ninitial == 1 or (init_mothers[0].get('number_external') == 1 \
+                             and init_mothers[0].get('leg_state') == False):
+                # For decay processes or if the mother is going
+                # towards external leg 1, mother leg is resulting wf
+                legs.append(mother_leg)
+            else:
+                # If this is an s-channel leg or we are going towards
+                # external leg 2, mother leg is one of the mothers
+                legs.insert(-1, mother_leg)
+
+            # Renumber resulting leg according to minimum leg number
+            legs[-1].set('number', min([l.get('number') for l in legs[:-1]]))
+
+            vertex = base_objects.Vertex({
+                'id': self.get('interaction_id'),
+                'legs': legs})
+
+            # Add s- and t-channels from further down
+            mother_s, tchannels = \
+                      init_mothers[0].get_s_and_t_channels(ninitial, legs[-1])
+
+            schannels.extend(mother_s)
+
+            if ninitial == 1 or init_mothers[0].get('leg_state') == True:
+                # This leg 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
+                tchannels.append(vertex)
+            else:
+                # If we are going towards external leg 2, add to
+                # t-channels, at start
+                tchannels.insert(0, vertex)
+
+        elif len(init_mothers) == 2:
+            # This is a t-channel junction. Start with the leg going
+            # towards external particle 1, and then do external
+            # particle 2
+            init_mothers1 = filter(lambda wf: wf.get('number_external') == 1,
+                                   init_mothers)[0]
+            init_mothers2 = filter(lambda wf: wf.get('number_external') == 2,
+                                   init_mothers)[0]
+
+            # Create vertex
+            legs = base_objects.LegList()
+            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
+                    }))
+            legs.insert(-1, mother_leg)
+
+            # Renumber resulting leg according to minimum leg number
+            legs[-1].set('number', min([l.get('number') for l in legs[:-1]]))
+
+            vertex = base_objects.Vertex({
+                'id': self.get('interaction_id'),
+                'legs': legs})
+
+            # Add s- and t-channels going down towards leg 1
+            mother_s, tchannels = \
+                      init_mothers1.get_s_and_t_channels(ninitial, legs[0])
+            schannels.extend(mother_s)
+
+            # Add vertex
+            tchannels.append(vertex)
+
+            # Add s- and t-channels going down towards leg 2
+            mother_s, mother_t = \
+                      init_mothers2.get_s_and_t_channels(ninitial, legs[-1])
+            schannels.extend(mother_s)
+            tchannels.extend(mother_t)
+
+        return schannels, tchannels
+
     # Overloaded operators
 
     def __eq__(self, other):
@@ -733,20 +975,24 @@
 
         # Check relevant directly defined properties
         if self['number_external'] != other['number_external'] or \
-           self['spin'] != other['spin'] or \
-           self['self_antipart'] != other['self_antipart'] or \
            self['fermionflow'] != other['fermionflow'] or \
-           self['mass'] != other['mass'] or \
-           self['width'] != other['width'] or \
-           self['color'] != other['color'] or \
+           self['coupl_key'] != other['coupl_key'] or \
            self['lorentz'] != other['lorentz'] or \
            self['coupling'] != other['coupling'] or \
-           self['state'] != other['state']:
+           self['state'] != other['state'] or \
+           self['onshell'] != other['onshell'] or \
+           self.get('spin') != other.get('spin') or \
+           self.get('self_antipart') != other.get('self_antipart') or \
+           self.get('mass') != other.get('mass') or \
+           self.get('width') != other.get('width') or \
+           self.get('color') != other.get('color') or \
+           self['decay'] != other['decay'] or \
+           self['decay'] and self['particle'] != other['particle']:
             return False
 
         # Check that mothers have the same numbers (only relevant info)
-        return [ mother.get('number') for mother in self['mothers'] ] == \
-               [ mother.get('number') for mother in other['mothers'] ]
+        return sorted([mother['number'] for mother in self['mothers']]) == \
+               sorted([mother['number'] for mother in other['mothers']])
 
     def __ne__(self, other):
         """Overloading the nonequality operator, to make comparison easy"""
@@ -768,12 +1014,16 @@
 
     # Helper functions
 
+    def to_array(self):
+        return array.array('i', [w['number'] for w in self])
+
     def check_and_fix_fermion_flow(self,
                                    wavefunctions,
                                    diagram_wavefunctions,
                                    external_wavefunctions,
                                    my_state,
-                                   wf_number):
+                                   wf_number,
+                                   force_flip_flow=False):
         """Check for clashing fermion flow (N(incoming) !=
         N(outgoing)). If found, we need to trace back through the
         mother structure (only looking at fermions), until we find a
@@ -813,6 +1063,11 @@
                       Please decompose your vertex into 2-fermion
                       vertices to get fermion flow correct."""
 
+        # Note that the order of mothers is given by the external
+        # number meaning that the order will in general be different
+        # for t- and s-channel diagrams, hence giving duplication of
+        # wave functions. I don't think there's anything to do about
+        # this.
         for mother in fermion_mothers:
             if Nincoming > Noutgoing and \
                mother.get_with_flow('state') == 'outgoing' or \
@@ -824,14 +1079,16 @@
 
             # 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)
+                                                wf_number,
+                                                force_flip_flow)
+
             # Replace old mother with new mother
             self[self.index(mother)] = new_mother
 
@@ -849,12 +1106,41 @@
                                        mother_states))
 
         if Nincoming != Noutgoing:
-            raise self.PhysicsObjectListError, \
-                  "Failed to fix fermion flow, %d != %d" % \
-                  (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)
 
         return wf_number
 
+    def insert_own_mothers(self):
+        """Recursively go through a wavefunction list and insert the
+        mothers of all wavefunctions, return the result.
+        Assumes that all wavefunctions have unique numbers."""
+
+        res = copy.copy(self)
+        # Recursively build up res
+        for wf in self:
+            index = res.index(wf)
+            res = res[:index] + wf.get('mothers').insert_own_mothers() \
+                  + res[index:]
+
+        # Make sure no wavefunctions occur twice, by removing doublets
+        # from the back
+        i = len(res) - 1
+        while res[:i]:
+            if res[i].get('number') in [w.get('number') for w in res[:i]]:
+                res.pop(i)
+            i = i - 1
+
+        return res
+
 #===============================================================================
 # HelasAmplitude
 #===============================================================================
@@ -873,6 +1159,8 @@
         self['inter_color'] = None
         self['lorentz'] = ''
         self['coupling'] = 'none'
+        # The Lorentz and color index used in this amplitude
+        self['coupl_key'] = (0, 0)
         # Properties relating to the vertex
         self['number'] = 0
         self['fermionfactor'] = 0
@@ -933,6 +1221,17 @@
                         "%s is not a valid coupling string" % \
                                                                 str(value)
 
+        if name == 'coupl_key':
+            if value and not isinstance(value, tuple):
+                raise self.PhysicsObjectError, \
+                      "%s is not a valid tuple" % str(value)
+            if len(value) != 2:
+                raise self.PhysicsObjectError, \
+                      "%s is not a valid tuple with 2 elements" % str(value)
+            if not isinstance(value[0], int) or not isinstance(value[1], int):
+                raise self.PhysicsObjectError, \
+                      "%s is not a valid tuple of integer" % str(value)
+
         if name == 'number':
             if not isinstance(value, int):
                 raise self.PhysicsObjectError, \
@@ -1019,8 +1318,8 @@
         """Return particle property names as a nicely sorted list."""
 
         return ['interaction_id', 'pdg_codes', 'inter_color', 'lorentz',
-                'coupling', 'number', 'color_indices', 'fermionfactor',
-                'mothers']
+                'coupling', 'coupl_key', 'number', 'color_indices',
+                'fermionfactor', 'mothers']
 
 
     # Helper functions
@@ -1039,9 +1338,10 @@
                                    wavefunctions,
                                    diagram_wavefunctions,
                                    external_wavefunctions,
-                                   'nostate',
+                                   "Nostate",
                                    wf_number)
 
+
     def needs_hermitian_conjugate(self):
         """Returns true if any of the mothers have negative
         fermionflow"""
@@ -1099,7 +1399,7 @@
             fermion_number_list.extend(in_list[1])
             fermion_number_list.extend(out_list[1])
         elif len(in_fermions) != len(out_fermions):
-            raise self.HelasWavefunctionError, \
+            raise self.PhysicsObjectError, \
                   "Error: %d incoming fermions != %d outgoing fermions" % \
                   (len(in_fermions), len(out_fermions))
 
@@ -1124,6 +1424,180 @@
 
         return (-1) ** nflips
 
+    def get_base_diagram(self, wf_dict = {}, vx_list = [], optimization = 1):
+        """Return the base_objects.Diagram which corresponds to this
+        amplitude, using a recursive method for the wavefunctions."""
+
+        vertices = base_objects.VertexList()
+
+        # Add vertices for all mothers
+        for mother in self.get('mothers'):
+            vertices.extend(mother.get_base_vertices(wf_dict, vx_list,
+                                                     optimization))
+        # Generate last vertex
+        legs = base_objects.LegList()
+        for mother in self.get('mothers'):
+            try:
+                leg = wf_dict[mother.get('number')]
+            except KeyError:
+                leg = base_objects.Leg({
+                    'id': mother.get_pdg_code(),
+                    'number': mother.get('number_external'),
+                    'state': mother.get('leg_state'),
+                    'from_group': mother.get('onshell')
+                    })
+                if optimization != 0:
+                    wf_dict[mother.get('number')] = leg
+            legs.append(leg)
+
+        vertices.append(base_objects.Vertex({
+            'id': self.get('interaction_id'),
+            'legs': legs}))
+
+        return base_objects.Diagram({'vertices': vertices})
+
+    def get_s_and_t_channels(self, ninitial):
+        """Returns two lists of vertices corresponding to the s- and
+        t-channels of this amplitude/diagram, ordered from the outermost
+        s-channel and in/down towards the highest number initial state
+        leg."""
+
+        schannels = base_objects.VertexList()
+        tchannels = base_objects.VertexList()
+
+        # Add vertices for all s-channel mothers
+        final_mothers = filter(lambda wf: wf.get('number_external') > ninitial,
+                               self.get('mothers'))
+
+        for mother in final_mothers:
+            schannels.extend(mother.get_base_vertices(optimization = 0))
+
+        # Extract initial state mothers
+        init_mothers = filter(lambda wf: wf.get('number_external') <= ninitial,
+                              self.get('mothers'))
+
+        if len(init_mothers) > 2:
+            raise self.PhysicsObjectError, \
+                  "get_s_and_t_channels can only handle up to 2 initial states"
+
+        if len(init_mothers) == 1:
+            # This is an s-channel leg. Add vertex and start stepping down
+            # towards initial state
+
+            # Create vertex
+            legs = base_objects.LegList()
+            for mother in final_mothers + init_mothers:
+                legs.append(base_objects.Leg({
+                    'id': mother.get_pdg_code(),
+                    'number': mother.get('number_external'),
+                    'state': mother.get('leg_state'),
+                    'from_group': False
+                    }))
+
+            # Renumber resulting leg according to minimum leg number
+            legs[-1].set('number', min([l.get('number') for l in legs[:-1]]))
+
+            # Add vertex to s-channels
+            schannels.append(base_objects.Vertex({
+                'id': self.get('interaction_id'),
+                'legs': legs}))
+
+            # Add s- and t-channels from further down
+            mother_s, tchannels = init_mothers[0].\
+                                  get_s_and_t_channels(ninitial, legs[-1])
+
+            schannels.extend(mother_s)
+        else:
+            # This is a t-channel leg. Start with the leg going
+            # towards external particle 1, and then do external
+            # particle 2
+            init_mothers1 = filter(lambda wf: wf.get('number_external') == 1,
+                                   init_mothers)[0]
+            init_mothers2 = filter(lambda wf: wf.get('number_external') == 2,
+                                   init_mothers)[0]
+
+            # Create vertex
+            legs = base_objects.LegList()
+            for mother in [init_mothers1] + final_mothers + [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]]))
+
+            vertex = base_objects.Vertex({
+                'id': self.get('interaction_id'),
+                'legs': legs})
+
+            # Add s- and t-channels going down towards leg 1
+            mother_s, tchannels = \
+                      init_mothers1.get_s_and_t_channels(ninitial, legs[0])
+
+            schannels.extend(mother_s)
+
+            # Add vertex to t-channels
+            tchannels.append(vertex)
+
+            # Add s- and t-channels going down towards leg 2
+            mother_s, mother_t = \
+                      init_mothers2.get_s_and_t_channels(ninitial, legs[-1])
+            schannels.extend(mother_s)
+            tchannels.extend(mother_t)
+
+        # Finally go through all vertices, sort the legs and replace
+        # leg number with propagator number -1, -2, ...
+        number_dict = {}
+        nprop = 0
+        for vertex in schannels + tchannels:
+            # Sort the legs
+            legs = vertex.get('legs')[:-1]
+            if vertex in schannels:
+                legs.sort(lambda l1, l2: l2.get('number') - \
+                          l1.get('number'))
+            else:
+                legs.sort(lambda l1, l2: l1.get('number') - \
+                          l2.get('number'))
+            for leg in legs:
+                try:
+                    leg.set('number', number_dict[leg.get('number')])
+                except KeyError:
+                    pass
+            nprop = nprop - 1
+            last_leg = vertex.get('legs')[-1]
+            number_dict[last_leg.get('number')] = nprop
+            last_leg.set('number', nprop)
+            legs.append(last_leg)
+            vertex.set('legs', base_objects.LegList(legs))
+
+        return schannels, tchannels
+
+    def get_color_indices(self):
+        """Get the color indices corresponding to
+        this amplitude and its mothers, using a recursive function."""
+
+        if not self.get('mothers'):
+            return []
+
+        color_indices = []
+
+        # Add color indices for all mothers
+        for mother in self.get('mothers'):
+            # This is where recursion happens
+            color_indices.extend(mother.get_color_indices())
+
+        # Add this amp's color index
+        if self.get('interaction_id'):
+            color_indices.append(self.get('coupl_key')[0])
+
+        return color_indices
+
     # Comparison between different amplitudes, to allow check for
     # identical processes. Note that we are then not interested in
     # interaction id, but in all other properties.
@@ -1142,8 +1616,8 @@
             return False
 
         # Check that mothers have the same numbers (only relevant info)
-        return [ mother.get('number') for mother in self['mothers'] ] == \
-               [ mother.get('number') for mother in other['mothers'] ]
+        return sorted([mother['number'] for mother in self['mothers']]) == \
+               sorted([mother['number'] for mother in other['mothers']])
 
     def __ne__(self, other):
         """Overloading the nonequality operator, to make comparison easy"""
@@ -1178,6 +1652,7 @@
         # different Lorentz or color structures associated with this
         # diagram
         self['amplitudes'] = HelasAmplitudeList()
+        self['number'] = 0
 
     def filter(self, name, value):
         """Filter for valid diagram property values."""
@@ -1241,8 +1716,13 @@
 
         self['processes'] = base_objects.ProcessList()
         self['diagrams'] = HelasDiagramList()
+        self['identical_particle_factor'] = 0
         self['color_basis'] = color_amp.ColorBasis()
         self['color_matrix'] = color_amp.ColorMatrix(color_amp.ColorBasis())
+        # base_amplitude is the Amplitude to be used in color
+        # generation, drawing etc. For decay chain processes, this is
+        # the Amplitude which corresponds to the combined process.
+        self['base_amplitude'] = None
 
     def filter(self, name, value):
         """Filter for valid diagram property values."""
@@ -1255,6 +1735,10 @@
             if not isinstance(value, HelasDiagramList):
                 raise self.PhysicsObjectError, \
                         "%s is not a valid HelasDiagramList object" % str(value)
+        if name == 'identical_particle_factor':
+            if not isinstance(value, int):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid int object" % str(value)
         if name == 'color_basis':
             if not isinstance(value, color_amp.ColorBasis):
                 raise self.PhysicsObjectError, \
@@ -1263,15 +1747,32 @@
             if not isinstance(value, color_amp.ColorMatrix):
                 raise self.PhysicsObjectError, \
                         "%s is not a valid ColorMatrix object" % str(value)
+        if name == 'base_amplitude':
+            if value != None and not \
+                   isinstance(value, diagram_generation.Amplitude):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid Amplitude object" % str(value)
         return True
 
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
 
-        return ['processes', 'diagrams', 'color_basis', 'color_matrix']
+        return ['processes', 'identical_particle_factor',
+                'diagrams', 'color_basis', 'color_matrix',
+                'base_amplitude']
+
+    # Enhanced get function
+    def get(self, name):
+        """Get the value of the property name."""
+
+        if name == 'base_amplitude' and not self[name]:
+            self['base_amplitude'] = self.get_base_amplitude()
+
+        return super(HelasMatrixElement, self).get(name)
 
     # Customized constructor
-    def __init__(self, amplitude=None, optimization=1, gen_color=True):
+    def __init__(self, amplitude=None, optimization=1,
+                 decay_ids=[], gen_color=True):
         """Constructor for the HelasMatrixElement. In particular allows
         generating a HelasMatrixElement from an Amplitude, with
         automatic generation of the necessary wavefunctions
@@ -1281,10 +1782,11 @@
             if isinstance(amplitude, diagram_generation.Amplitude):
                 super(HelasMatrixElement, self).__init__()
                 self.get('processes').append(amplitude.get('process'))
-                self.generate_helas_diagrams(amplitude, optimization)
+                self.generate_helas_diagrams(amplitude, optimization, decay_ids)
                 self.calculate_fermionfactors()
+                self.calculate_identical_particle_factors()
                 if gen_color:
-                    self.get('color_basis').build(amplitude)
+                    self.get('color_basis').build(self.get('base_amplitude'))
                     self.set('color_matrix',
                              color_amp.ColorMatrix(self.get('color_basis')))
             else:
@@ -1304,10 +1806,20 @@
         if not isinstance(other, HelasMatrixElement):
             return False
 
+        # If no processes, this is an empty matrix element
+        if not self['processes'] and not other['processes']:
+            return True
+
         # Should only check if diagrams and process id are identical
+        # Except in case of decay processes: then also initial state
+        # must be the same
         if self['processes'] and not other['processes'] or \
                self['processes'] and \
                self['processes'][0]['id'] != other['processes'][0]['id'] or \
+               self['processes'][0]['is_decay_chain'] or \
+               other['processes'][0]['is_decay_chain'] or \
+               self['identical_particle_factor'] != \
+                           other['identical_particle_factor'] or \
                self['diagrams'] != other['diagrams']:
             return False
 
@@ -1317,13 +1829,17 @@
         """Overloading the nonequality operator, to make comparison easy"""
         return not self.__eq__(other)
 
-    def generate_helas_diagrams(self, amplitude, optimization=1):
+    def generate_helas_diagrams(self, amplitude, optimization=1,
+                                decay_ids=[]):
         """Starting from a list of Diagrams from the diagram
         generation, generate the corresponding HelasDiagrams, i.e.,
         the wave functions and amplitudes. Choose between default
         optimization (= 1, maximum recycling of wavefunctions) or no
         optimization (= 0, no recycling of wavefunctions, useful for
         GPU calculations with very restricted memory).
+
+        Note that we need special treatment for decay chains, since
+        the end product then is a wavefunction, not an amplitude.
         """
 
         if not isinstance(amplitude, diagram_generation.Amplitude) or \
@@ -1333,6 +1849,7 @@
 
         diagram_list = amplitude.get('diagrams')
         process = amplitude.get('process')
+
         model = process.get('model')
         if not diagram_list:
             return
@@ -1347,26 +1864,36 @@
 
         # Generate wavefunctions for the external particles
         external_wavefunctions = dict([(leg.get('number'),
-                                        HelasWavefunction(leg, 0, model)) \
+                                        HelasWavefunction(leg, 0, model,
+                                                          decay_ids)) \
                                        for leg in process.get('legs')])
+
         # Initially, have one wavefunction for each external leg.
         wf_number = len(process.get('legs'))
 
-        # For initial state bosons, need to flip PDG code (if has antipart)
+        # For initial state bosons, need to flip part-antipart
         # since all bosons should be treated as outgoing
         for key in external_wavefunctions.keys():
             wf = external_wavefunctions[key]
             if wf.is_boson() and wf.get('state') == 'initial' and \
                not wf.get('self_antipart'):
-                wf.set('pdg_code', -wf.get('pdg_code'))
                 wf.set('is_part', not wf.get('is_part'))
 
+        # For initial state particles, need to flip PDG code (if has
+        # antipart)
+        for key in external_wavefunctions.keys():
+            wf = external_wavefunctions[key]
+            if wf.get('leg_state') == False and \
+               not wf.get('self_antipart'):
+                wf.flip_part_antipart()
+
         # Now go through the diagrams, looking for undefined wavefunctions
 
         helas_diagrams = HelasDiagramList()
 
-        # Keep track of amplitude number
+        # Keep track of amplitude number and diagram number
         amplitude_number = 0
+        diagram_number = 0
 
         for diagram in diagram_list:
 
@@ -1387,7 +1914,7 @@
             lastvx = vertices.pop()
 
             # Check if last vertex is identity vertex
-            if lastvx.get('id') == 0:
+            if not process.get('is_decay_chain') and lastvx.get('id') == 0:
                 # Need to "glue together" last and next-to-last
                 # vertext, by replacing the (incoming) last leg of the
                 # next-to-last vertex with the (outgoing) leg in the
@@ -1432,6 +1959,7 @@
                     # Need one amplitude for each Lorentz/color structure,
                     # i.e. for each coupling
                     for coupl_key in inter.get('couplings').keys():
+
                         wf = HelasWavefunction(last_leg, vertex.get('id'), model)
                         wf.set('coupling', inter.get('couplings')[coupl_key])
                         # Special feature: For HVS vertices with the two
@@ -1462,7 +1990,10 @@
                         new_number_wf_dict = copy.copy(number_wf_dict)
 
                         # Store wavefunction
-                        if not wf in diagram_wavefunctions:
+                        try:
+                            wf = diagram_wavefunctions[\
+                                    diagram_wavefunctions.index(wf)]
+                        except ValueError:
                             # Update wf number
                             wf_number = wf_number + 1
                             wf.set('number', wf_number)
@@ -1477,7 +2008,7 @@
                             except ValueError:
                                 diagram_wavefunctions.append(wf)
 
-                            new_number_wf_dict[last_leg.get('number')] = wf
+                        new_number_wf_dict[last_leg.get('number')] = wf
 
                         # Store the new copy of number_wf_dict
                         new_number_to_wavefunctions.append(\
@@ -1493,6 +2024,8 @@
             # Generate all amplitudes corresponding to the different
             # copies of this diagram
             helas_diagram = HelasDiagram()
+            diagram_number = diagram_number + 1
+            helas_diagram.set('number', diagram_number)
             for number_wf_dict, color_list in zip(number_to_wavefunctions,
                                                   color_lists):
                 # Find mothers for the amplitude
@@ -1506,29 +2039,42 @@
                 wf_number = mothers.check_and_fix_fermion_flow(wavefunctions,
                                               diagram_wavefunctions,
                                               external_wavefunctions,
-                                              'nostate',
+                                              "Nostate",
                                               wf_number)
+
                 # Sort the wavefunctions according to number
                 diagram_wavefunctions.sort(lambda wf1, wf2: \
                               wf1.get('number') - wf2.get('number'))
 
                 # Now generate HelasAmplitudes from the last vertex.
-                inter = model.get('interaction_dict')[lastvx.get('id')]
-                for i, coupl_key in enumerate(inter.get('couplings').keys()):
+                if lastvx.get('id'):
+                    inter = model.get('interaction_dict')[lastvx.get('id')]
+                    keys = 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)
-                    amp.set('coupling', inter.get('couplings')[coupl_key])
-                    if inter.get('color'):
-                        amp.set('inter_color', inter.get('color')[coupl_key[0]])
-                    amp.set('lorentz', inter.get('lorentz')[coupl_key[1]])
+                    if inter:
+                        amp.set('coupling', inter.get('couplings')[coupl_key])
+                        amp.set('lorentz', inter.get('lorentz')[\
+                                coupl_key[1]])
+                        if inter.get('color'):
+                            amp.set('inter_color', inter.get('color')[\
+                                coupl_key[0]])
+                        amp.set('coupl_key', coupl_key)
                     amp.set('mothers', mothers)
                     amplitude_number = amplitude_number + 1
                     amp.set('number', amplitude_number)
                     # Add the list with color indices to the amplitude
                     new_color_list = copy.copy(color_list)
-                    new_color_list.append(coupl_key[0])
+                    if inter:
+                        new_color_list.append(coupl_key[0])
                     amp.set('color_indices', new_color_list)
+
                     # Generate HelasDiagram
-
                     helas_diagram.get('amplitudes').append(amp)
                     if diagram_wavefunctions and not \
                                        helas_diagram.get('wavefunctions'):
@@ -1546,15 +2092,714 @@
 
 
         self.set('diagrams', helas_diagrams)
+        # Sort all mothers according to the order wanted in Helas calls
+
+        for wf in self.get_all_wavefunctions():
+            wf.set('mothers', HelasMatrixElement.sorted_mothers(wf))
+        for amp in self.get_all_amplitudes():
+            amp.set('mothers', HelasMatrixElement.sorted_mothers(amp))
+            amp.set('color_indices', amp.get_color_indices())
+
+    def insert_decay_chains(self, decay_dict):
+        """Iteratively insert decay chains decays into this matrix
+        element.        
+        * decay_dict: a dictionary from external leg number
+          to decay matrix element.
+        """
+
+        # We need to keep track of how the
+        # wavefunction numbers change
+        replace_dict = {}
+        for number in decay_dict.keys():
+            # Find all wavefunctions corresponding to this external
+            # leg number
+            replace_dict[number] = [wf for wf in \
+                          filter(lambda wf: not wf.get('mothers') and \
+                                 wf.get('number_external') == number,
+                                 self.get_all_wavefunctions())]
+
+        # Keep track of wavefunction and amplitude numbers, to ensure
+        # unique numbers for all new wfs and amps during manipulations
+        numbers = [self.get_all_wavefunctions()[-1].get('number'),
+                   self.get_all_amplitudes()[-1].get('number')]
+
+        # Check if there are any Majorana particles present in any diagrams
+        got_majoranas = False
+        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():
+                got_majoranas = True
+
+        # Now insert decays for all legs that have decays
+        for number in decay_dict.keys():
+
+                self.insert_decay(replace_dict[number],
+                                  decay_dict[number],
+                                  numbers,
+                                  got_majoranas)
+
+        # Final cleaning out duplicate wavefunctions - needed only if
+        # we have multiple fermion flows, i.e., either multiple replaced
+        # wavefunctions or  majorana fermions and multiple diagrams
+        flows = reduce(lambda i1, i2: i1 * i2,
+                       [len(replace_dict[i]) for i in decay_dict.keys()], 1)
+        diagrams = reduce(lambda i1, i2: i1 * i2,
+                               [len(decay_dict[i].get('diagrams')) for i in \
+                                decay_dict.keys()], 1)
+
+        if flows > 1 or (diagrams > 1 and got_majoranas):
+
+            # Clean out any doublet wavefunctions
+
+            earlier_wfs = []
+
+            earlier_wf_arrays = []
+
+            mothers = self.get_all_wavefunctions() + self.get_all_amplitudes()
+            mother_arrays = [w['mothers'].to_array() \
+                             for w in mothers]
+
+            for diagram in self.get('diagrams'):
+
+                if diagram.get('number') > 1:
+                    earlier_wfs.extend(self.get('diagrams')[\
+                        diagram.get('number') - 2].get('wavefunctions'))
+
+                i = 0
+                diag_wfs = diagram.get('wavefunctions')
+
+
+                # Remove wavefunctions and replace mothers
+                while diag_wfs[i:]:
+                    try:
+                        new_wf = earlier_wfs[\
+                            earlier_wfs.index(diag_wfs[i])]
+                        wf = diag_wfs.pop(i)
+
+                        self.update_later_mothers(wf, new_wf, mothers,
+                                                  mother_arrays)
+                    except ValueError:
+                        i = i + 1
+
+        # When we are done with all decays, set wavefunction and
+        # amplitude numbers
+        for i, wf in enumerate(self.get_all_wavefunctions()):
+            wf.set('number', i + 1)
+        for i, amp in enumerate(self.get_all_amplitudes()):
+            amp.set('number', i + 1)
+            # Update fermion factors for all amplitudes
+            amp.calculate_fermionfactor()
+            # Update color indices
+            amp.set('color_indices', amp.get_color_indices())
+
+        # Calculate identical particle factors for
+        # this matrix element
+        self.identical_decay_chain_factor(decay_dict.values())
+
+    def insert_decay(self, old_wfs, decay, numbers, got_majoranas):
+        """Insert a decay chain matrix element into the matrix element.
+        * old_wfs: the wavefunctions to be replaced.
+          They all correspond to the same external particle, but might
+          have different fermion flow directions
+        * decay: the matrix element for the decay chain
+        * numbers: the present wavefunction and amplitude number,
+          to allow for unique numbering
+          
+        Note that:
+        1) All amplitudes and all wavefunctions using the decaying wf
+           must be copied as many times as there are amplitudes in the
+           decay matrix element
+        2) In the presence of Majorana particles, we must make sure
+           to flip fermion flow for the decay process if needed.
+
+        The algorithm is the following:
+        1) Multiply the diagrams with the number of diagrams Ndiag in
+           the decay element
+        2) For each diagram in the decay element, work on the diagrams
+           which corresponds to it
+        3) Flip fermion flow for the decay wavefunctions if needed
+        4) Insert all auxiliary wavefunctions into the diagram (i.e., all 
+           except the final wavefunctions, which directly replace the
+           original final state wavefunctions)
+        4) Replace the wavefunctions recursively, so that we always replace
+           each old wavefunctions with Namp new ones, where Namp is
+           the number of amplitudes in this decay element
+           diagram. Do recursion for wavefunctions which have this
+           wavefunction as mother. Simultaneously replace any
+           amplitudes which have this wavefunction as mother.
+        """
+
+        len_decay = len(decay.get('diagrams'))
+
+        number_external = old_wfs[0].get('number_external')
+
+        # Insert the decay process in the process
+        for process in self.get('processes'):
+            process.get('decay_chains').append(\
+                   decay.get('processes')[0])
+
+        # We need one copy of the decay element diagrams for each
+        # old_wf to be replaced, since we need different wavefunction
+        # numbers for them
+        decay_elements = [copy.deepcopy(d) for d in \
+                          [ decay.get('diagrams') ] * len(old_wfs)]
+
+        # Need to replace Particle in all wavefunctions to avoid
+        # deepcopy
+        idecay = 0
+        for decay_element in decay_elements:
+            for idiag, diagram in enumerate(decay.get('diagrams')):
+                wfs = diagram.get('wavefunctions')
+                decay_diag = decay_element[idiag]
+                for i, wf in enumerate(decay_diag.get('wavefunctions')):
+                    wf.set('particle', wfs[i].get('particle'))
+                    wf.set('antiparticle', wfs[i].get('antiparticle'))
+
+        for decay_element in decay_elements:
+
+            # Remove the unwanted initial state wavefunctions from decay
+            for decay_diag in decay_element:
+                for wf in filter(lambda wf: wf.get('number_external') == 1,
+                                 decay_diag.get('wavefunctions')):
+                    decay_diag.get('wavefunctions').remove(wf)
+
+            decay_wfs = sum([d.get('wavefunctions') for d in decay_element], [])
+
+            # External wavefunction offset for new wfs
+            incr_new = number_external - \
+                       decay_wfs[0].get('number_external')
+
+            for wf in decay_wfs:
+                # Set number_external for new wavefunctions
+                wf.set('number_external', wf.get('number_external') + incr_new)
+                # Set unique number for new wavefunctions
+                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()
+        incr_old = nex - 2 # due to the incoming particle
+        wavefunctions = self.get_all_wavefunctions()
+        for wf in wavefunctions:
+            # Only modify number_external for wavefunctions above old_wf
+            if wf.get('number_external') > number_external:
+                wf.set('number_external',
+                       wf.get('number_external') + incr_old)
+
+        # Multiply the diagrams by Ndiag
+
+        diagrams = HelasDiagramList()
+        for diagram in self.get('diagrams'):
+            new_diagrams = [copy.copy(diag) for diag in \
+                            [ diagram ] * (len_decay - 1)]
+            # Update diagram number
+            diagram.set('number', (diagram.get('number') - 1) * \
+                        len_decay + 1)
+
+            for i, diag in enumerate(new_diagrams):
+                # Set diagram number
+                diag.set('number', diagram.get('number') + i + 1)
+                # Copy over all wavefunctions
+                diag.set('wavefunctions',
+                         copy.copy(diagram.get('wavefunctions')))
+                # Copy over the amplitudes
+                amplitudes = HelasAmplitudeList(\
+                                [copy.copy(amp) for amp in \
+                                 diag.get('amplitudes')])
+                # Renumber amplitudes
+                for amp in amplitudes:
+                    numbers[1] = numbers[1] + 1
+                    amp.set('number', numbers[1])
+                diag.set('amplitudes', amplitudes)
+            # Add old and new diagram to diagrams
+            diagrams.append(diagram)
+            diagrams.extend(new_diagrams)
+
+        self.set('diagrams', diagrams)
+
+        # Now we work by decay process diagram, parameterized by numdecay
+        for numdecay in range(len_decay):
+
+            # Pick out the diagrams which correspond to this decay diagram
+            diagrams = [self.get('diagrams')[i] for i in \
+                        range(numdecay, len(self.get('diagrams')), len_decay)]
+
+            # Perform replacement for each of the wavefunctions in old_wfs
+            for decay_element, old_wf in zip(decay_elements, old_wfs):
+
+                decay_diag = decay_element[numdecay]
+
+                # Find the diagrams which have old_wf
+                my_diagrams = filter(lambda diag: (old_wf.get('number') in \
+                                            [wf.get('number') for wf in \
+                                            diag.get('wavefunctions')]),
+                                     diagrams)
+
+                # Ignore possibility for unoptimizated generation for now
+                if len(my_diagrams) > 1:
+                    raise self.PhysicsObjectError, \
+                          "Decay chains not yet prepared for GPU"
+
+                for diagram in my_diagrams:
+
+                    if got_majoranas:
+                        # If there are Majorana particles in any of
+                        # the matrix elements, we need to check for
+                        # fermion flow
+
+                        # Earlier wavefunctions, will be used for fermion flow
+                        index = [d.get('number') for d in diagrams].\
+                                index(diagram.get('number'))
+                        earlier_wavefunctions = \
+                                      sum([d.get('wavefunctions') for d in \
+                                           diagrams[:index]], [])
+
+                        # Don't want to affect original decay
+                        # wavefunctions, so need to deepcopy
+                        decay_diag_wfs = copy.deepcopy(\
+                                                decay_diag.get('wavefunctions'))
+                        # Need to replace Particle in all
+                        # wavefunctions to avoid deepcopy
+                        idecay = 0
+                        for i, wf in enumerate(decay_diag.get('wavefunctions')):
+                            decay_diag_wfs[i].set('particle', \
+                                                  wf.get('particle'))
+                            decay_diag_wfs[i].set('antiparticle', \
+                                                  wf.get('antiparticle'))
+
+                        # Complete decay_diag_wfs with the mother wavefunctions
+                        # to allow for independent fermion flow flips
+                        decay_diag_wfs = decay_diag_wfs.insert_own_mothers()
+
+                        # These are the wavefunctions which directly replace old_wf
+                        final_decay_wfs = [amp.get('mothers')[1] for amp in \
+                                              decay_diag.get('amplitudes')]
+
+                        # Since we made deepcopy, need to syncronize
+                        for i, wf in enumerate(final_decay_wfs):
+                            final_decay_wfs[i] = \
+                                               decay_diag_wfs[decay_diag_wfs.index(wf)]
+                        # Remove final wavefunctions from decay_diag_wfs,
+                        # since these will be replaced separately by
+                        # replace_wavefunctions
+                        for wf in final_decay_wfs:
+                            decay_diag_wfs.remove(wf)
+
+                        # Check fermion flow direction
+                        if old_wf.is_fermion() and \
+                               old_wf.get_with_flow('state') != \
+                                     final_decay_wfs[0].get_with_flow('state'):
+
+                            # Not same flow state - need to flip flow of wf
+
+                            for i, wf in enumerate(final_decay_wfs):
+
+                                # We use the function
+                                # check_majorana_and_flip_flow, as in the
+                                # helas diagram generation.  Since we have
+                                # different flow, there is already a Majorana
+                                # particle along the fermion line.
+
+                                final_decay_wfs[i], numbers[0] = \
+                                                wf.check_majorana_and_flip_flow(\
+                                                         True,
+                                                         earlier_wavefunctions,
+                                                         decay_diag_wfs,
+                                                         {},
+                                                         numbers[0])
+
+                        # Remove wavefunctions which are already present in
+                        # earlier_wavefunctions
+                        i = 0
+                        earlier_wavefunctions = \
+                            sum([d.get('wavefunctions') for d in \
+                                 self.get('diagrams')[:diagram.get('number') - 1]], \
+                                [])
+                        earlier_wf_numbers = [wf.get('number') for wf in \
+                                              earlier_wavefunctions]
+
+                        i = 0
+                        mother_arrays = [w.get('mothers').to_array() for \
+                                         w in final_decay_wfs]
+                        while decay_diag_wfs[i:]:
+                            wf = decay_diag_wfs[i]
+                            try:
+                                new_wf = earlier_wavefunctions[\
+                                    earlier_wf_numbers.index(wf.get('number'))]
+                                # If the wavefunctions are not identical,
+                                # then we should keep this wavefunction,
+                                # and update its number so it is unique
+                                if wf != new_wf:
+                                    numbers[0] = numbers[0] + 1
+                                    wf.set('number', numbers[0])
+                                    continue
+                                decay_diag_wfs.pop(i)
+                                pres_mother_arrays = [w.get('mothers').to_array() for \
+                                                      w in decay_diag_wfs[i:]] + \
+                                                      mother_arrays
+                                self.update_later_mothers(wf, new_wf,
+                                                          decay_diag_wfs[i:] + \
+                                                          final_decay_wfs,
+                                                          pres_mother_arrays)
+                            except ValueError:
+                                i = i + 1
+
+                        # Since we made deepcopy, go through mothers and make
+                        # sure we are using the ones in earlier_wavefunctions
+                        for decay_wf in decay_diag_wfs + final_decay_wfs:
+                            mothers = decay_wf.get('mothers')
+                            for i, wf in enumerate(mothers):
+                                try:
+                                    mothers[i] = earlier_wavefunctions[\
+                                        earlier_wf_numbers.index(wf.get('number'))]
+                                except ValueError:
+                                    pass
+                    else:
+                        # If there are no Majorana particles, the
+                        # treatment is much simpler
+                        decay_diag_wfs = \
+                                       copy.copy(decay_diag.get('wavefunctions'))
+
+                        # These are the wavefunctions which directly
+                        # replace old_wf
+                        final_decay_wfs = [amp.get('mothers')[1] for amp in \
+                                              decay_diag.get('amplitudes')]
+
+                        # Remove final wavefunctions from decay_diag_wfs,
+                        # since these will be replaced separately by
+                        # replace_wavefunctions
+                        for wf in final_decay_wfs:
+                            decay_diag_wfs.remove(wf)
+
+
+                    diagram_wfs = diagram.get('wavefunctions')
+
+                    old_wf_index = [wf.get('number') for wf in \
+                                    diagram_wfs].index(old_wf.get('number'))
+
+                    diagram_wfs = diagram_wfs[0:old_wf_index] + \
+                                  decay_diag_wfs + diagram_wfs[old_wf_index:]
+
+                    diagram.set('wavefunctions', HelasWavefunctionList(diagram_wfs))
+
+                    # Set the decay flag for final_decay_wfs, to
+                    # indicate that these correspond to decayed
+                    # particles
+                    for wf in final_decay_wfs:
+                        wf.set('onshell', True)
+
+                    if len_decay == 1:
+                        # Can use simplified treatment, by just modifying old_wf
+                        self.replace_single_wavefunction(old_wf,
+                                                         final_decay_wfs[0])
+                    else:
+                        # Multiply wavefunctions and amplitudes and
+                        # insert mothers using a recursive function
+                        self.replace_wavefunctions(old_wf,
+                                                   final_decay_wfs,
+                                                   diagrams,
+                                                   numbers)
+
+            # Now that we are done with this set of diagrams, we need
+            # to clean out duplicate wavefunctions (i.e., remove
+            # identical wavefunctions which are already present in
+            # earlier diagrams)
+            for diagram in diagrams:
+
+                # We can have duplicate wfs only from previous copies of
+                # this diagram
+                earlier_wfs = sum([d.get('wavefunctions') for d in \
+                                   self.get('diagrams')[\
+                                     diagram.get('number') - numdecay - 1:\
+                                     diagram.get('number') - 1]], [])
+
+                later_wfs = sum([d.get('wavefunctions') for d in \
+                                   self.get('diagrams')[\
+                                     diagram.get('number'):]], [])
+
+                later_amps = sum([d.get('amplitudes') for d in \
+                                   self.get('diagrams')[\
+                                     diagram.get('number') - 1:]], [])
+
+                i = 0
+                diag_wfs = diagram.get('wavefunctions')
+
+                # Remove wavefunctions and replace mothers, to make
+                # sure we only have one copy of each wavefunction
+                # number
+
+                mother_arrays = [w.get('mothers').to_array() for \
+                                 w in later_wfs + later_amps]
+
+                while diag_wfs[i:]:
+                    try:
+                        index = [w.get('number') for w in earlier_wfs].\
+                                index(diag_wfs[i].get('number'))
+                        wf = diag_wfs.pop(i)
+                        pres_mother_arrays = [w.get('mothers').to_array() for \
+                                              w in diag_wfs[i:]] + \
+                                              mother_arrays
+                        self.update_later_mothers(wf, earlier_wfs[index],
+                                              diag_wfs[i:] + later_wfs + later_amps,
+                                              pres_mother_arrays)
+                    except ValueError:
+                        i = i + 1
+
+    def update_later_mothers(self, wf, new_wf, later_wfs, later_wf_arrays):
+        """Update mothers for all later wavefunctions"""
+
+        daughters = filter(lambda tup: wf.get('number') in tup[1],
+                              enumerate(later_wf_arrays))
+
+        for (index, mothers) in daughters:
+            try:
+                # Replace mother
+                later_wfs[index].get('mothers')[\
+                    mothers.index(wf.get('number'))] = new_wf
+            except ValueError:
+                pass
+
+    def replace_wavefunctions(self, old_wf, new_wfs,
+                              diagrams, numbers):
+        """Recursive function to replace old_wf with new_wfs, and
+        multiply all wavefunctions or amplitudes that use old_wf
+
+        * old_wf: The wavefunction to be replaced
+        * new_wfs: The replacing wavefunction
+        * diagrams - the diagrams that are relevant for these new
+          wavefunctions.
+        * numbers: the present wavefunction and amplitude number,
+          to allow for unique numbering
+        """
+
+        # Pick out the diagrams which have the old_wf
+        my_diagrams = filter(lambda diag: old_wf.get('number') in \
+                         [wf.get('number') for wf in diag.get('wavefunctions')],
+                         diagrams)
+
+        # Replace old_wf with new_wfs in the diagrams
+        for diagram in my_diagrams:
+
+            diagram_wfs = diagram.get('wavefunctions')
+
+            old_wf_index = [wf.get('number') for wf in \
+                            diagram.get('wavefunctions')].index(old_wf.get('number'))
+            diagram_wfs = diagram_wfs[:old_wf_index] + \
+                          new_wfs + diagram_wfs[old_wf_index + 1:]
+
+            diagram.set('wavefunctions', HelasWavefunctionList(diagram_wfs))
+
+        # Now need to take care of amplitudes and wavefunctions which
+        # are daughters of old_wf (only among the relevant diagrams)
+
+        # Pick out diagrams with amplitudes which are daughters of old_wf
+        amp_diagrams = filter(lambda diag: old_wf.get('number') in \
+                          sum([[wf.get('number') for wf in amp.get('mothers')] \
+                               for amp in diag.get('amplitudes')], []),
+                              diagrams)
+
+        for diagram in amp_diagrams:
+
+            # Amplitudes in this diagram that are daughters of old_wf
+            daughter_amps = filter(lambda amp: old_wf.get('number') in \
+                                [wf.get('number') for wf in amp.get('mothers')],
+                                diagram.get('amplitudes'))
+
+            new_amplitudes = copy.copy(diagram.get('amplitudes'))
+
+            # Loop over daughter_amps, to multiply each amp by the
+            # number of replacement wavefunctions and substitute mothers
+            for old_amp in daughter_amps:
+                # Create copies of this amp
+                new_amps = [copy.copy(amp) for amp in \
+                            [ old_amp ] * len(new_wfs)]
+                # Replace the old mother with the new ones
+                for i, (new_amp, new_wf) in enumerate(zip(new_amps, new_wfs)):
+                    mothers = copy.copy(new_amp.get('mothers'))
+                    old_wf_index = [wf.get('number') for wf in mothers].index(\
+                         old_wf.get('number'))
+                    # Update mother
+                    mothers[old_wf_index] = new_wf
+                    new_amp.set('mothers', mothers)
+                    # Update amp numbers for replaced amp
+                    numbers[1] = numbers[1] + 1
+                    new_amp.set('number', numbers[1])
+
+                # Insert the new amplitudes in diagram amplitudes
+                index = [a.get('number') for a in new_amplitudes].\
+                                   index(old_amp.get('number'))
+                new_amplitudes = new_amplitudes[:index] + \
+                                 new_amps + new_amplitudes[index + 1:]
+
+            # Replace diagram amplitudes with the new ones
+            diagram.set('amplitudes', HelasAmplitudeList(new_amplitudes))
+
+        # Find wavefunctions that are daughters of old_wf
+        daughter_wfs = filter(lambda wf: old_wf.get('number') in \
+                              [wf1.get('number') for wf1 in wf.get('mothers')],
+                              sum([diag.get('wavefunctions') for diag in \
+                                   diagrams], []))
+
+        # Loop over daughter_wfs, multiply them and replace mothers
+        for daughter_wf in daughter_wfs:
+
+            # Pick out the diagrams where daughter_wf occurs
+            wf_diagrams = filter(lambda diag: daughter_wf.get('number') in \
+                                 [wf.get('number') for wf in \
+                                  diag.get('wavefunctions')],
+                                 diagrams)
+
+            if len(wf_diagrams) > 1:
+                raise self.PhysicsObjectError, \
+                      "Decay chains not yet prepared for GPU"
+
+            for diagram in wf_diagrams:
+
+                # Now create new wfs with updated mothers
+                replace_daughters = [ copy.copy(wf) for wf in \
+                                      [daughter_wf] * len(new_wfs) ]
+
+                index = [wf.get('number') for wf in \
+                         daughter_wf.get('mothers')].index(old_wf.get('number'))
+
+                # Replace the old mother with the new ones, update wf numbers
+                for i, (new_daughter, new_wf) in \
+                        enumerate(zip(replace_daughters, new_wfs)):
+                    mothers = copy.copy(new_daughter.get('mothers'))
+                    mothers[index] = new_wf
+                    new_daughter.set('mothers', mothers)
+                    numbers[0] = numbers[0] + 1
+                    new_daughter.set('number', numbers[0])
+
+                # This is where recursion happens.  We need to replace
+                # the daughter wavefunction, and fix amplitudes and
+                # wavefunctions which have it as mothers.
+
+                self.replace_wavefunctions(daughter_wf,
+                                           replace_daughters,
+                                           diagrams,
+                                           numbers)
+
+    def replace_single_wavefunction(self, old_wf, new_wf):
+        """Insert decay chain by simply modifying wavefunction. This
+        is possible only if there is only one diagram in the decay."""
+
+        for key in old_wf.keys():
+            old_wf.set(key, new_wf.get(key))
+
+    def identical_decay_chain_factor(self, decay_chains):
+        """Calculate the denominator factor from identical decay chains"""
+
+        final_legs = [leg.get('id') for leg in \
+                      filter(lambda leg: leg.get('state') == True, \
+                              self.get('processes')[0].get('legs'))]
+
+        # Leg ids for legs being replaced by decay chains
+        decay_ids = [decay.get('legs')[0].get('id') for decay in \
+                     self.get('processes')[0].get('decay_chains')]
+
+        # Find all leg ids which are not being replaced by decay chains
+        non_decay_legs = filter(lambda id: id not in decay_ids,
+                                final_legs)
+
+        # Identical particle factor for legs not being decayed
+        identical_indices = {}
+        for id in non_decay_legs:
+            if id in identical_indices:
+                identical_indices[id] = \
+                                    identical_indices[id] + 1
+            else:
+                identical_indices[id] = 1
+        non_chain_factor = reduce(lambda x, y: x * y,
+                                  [ math.factorial(val) for val in \
+                                    identical_indices.values() ], 1)
+
+        # Identical particle factor for decay chains
+        # Go through chains to find identical ones
+        chains = copy.copy(decay_chains)
+        iden_chains_factor = 1
+        while chains:
+            ident_copies = 1
+            first_chain = chains.pop(0)
+            i = 0
+            while i < len(chains):
+                chain = chains[i]
+                if HelasMatrixElement.check_equal_decay_processes(\
+                                                 first_chain, chain):
+                    ident_copies = ident_copies + 1
+                    chains.pop(i)
+                else:
+                    i = i + 1
+            iden_chains_factor = iden_chains_factor * \
+                                 math.factorial(ident_copies)
+
+        self['identical_particle_factor'] = non_chain_factor * \
+                                    iden_chains_factor * \
+                                    reduce(lambda x1, x2: x1 * x2,
+                                    [me.get('identical_particle_factor') \
+                                     for me in decay_chains], 1)
 
     def calculate_fermionfactors(self):
-        """Generate the fermion factors for all diagrams in the amplitude
+        """Generate the fermion factors for all diagrams in the matrix element
         """
 
         for diagram in self.get('diagrams'):
             for amplitude in diagram.get('amplitudes'):
                 amplitude.get('fermionfactor')
 
+    def calculate_identical_particle_factors(self):
+        """Calculate the denominator factor for identical final state particles
+        """
+
+        final_legs = filter(lambda leg: leg.get('state') == True, \
+                              self.get('processes')[0].get('legs'))
+
+        identical_indices = {}
+        for leg in final_legs:
+            if leg.get('id') in identical_indices:
+                identical_indices[leg.get('id')] = \
+                                    identical_indices[leg.get('id')] + 1
+            else:
+                identical_indices[leg.get('id')] = 1
+        self["identical_particle_factor"] = reduce(lambda x, y: x * y,
+                                          [ math.factorial(val) for val in \
+                                            identical_indices.values() ])
+
+    def get_base_amplitude(self):
+        """Generate a diagram_generation.Amplitude from a
+        HelasMatrixElement. This is used to generate both color
+        amplitudes and diagram drawing."""
+
+        # Need to take care of diagram numbering for decay chains
+        # before this can be used for those!
+
+        optimization = 1
+        if len(filter(lambda wf: wf.get('number') == 1,
+                      self.get_all_wavefunctions())) > 1:
+            optimization = 0
+
+        model = self.get('processes')[0].get('model')
+
+        wf_dict = {}
+        vx_list = []
+        diagrams = base_objects.DiagramList()
+        for diag in self.get('diagrams'):
+            diagrams.append(diag.get('amplitudes')[0].get_base_diagram(\
+                wf_dict, vx_list, optimization))
+
+        return diagram_generation.Amplitude({\
+            'process': self.get('processes')[0],
+            'diagrams': diagrams})
+
     # Helper methods
 
     def getmothers(self, legs, number_to_wavefunctions,
@@ -1585,6 +2830,36 @@
         return sum([ len(d.get('wavefunctions')) for d in \
                        self.get('diagrams')])
 
+    def get_all_wavefunctions(self):
+        """Gives a list of all wavefunctions for this ME"""
+
+        return sum([d.get('wavefunctions') for d in \
+                       self.get('diagrams')], [])
+
+    def get_all_amplitudes(self):
+        """Gives a list of all amplitudes for this ME"""
+
+        return sum([d.get('amplitudes') for d in \
+                       self.get('diagrams')], [])
+
+    def get_external_wavefunctions(self):
+        """Gives the external wavefunctions for this ME"""
+
+        external_wfs = filter(lambda wf: not wf.get('mothers'),
+                              self.get('diagrams')[0].get('wavefunctions'))
+
+        external_wfs.sort(lambda w1, w2: w1.get('number_external') - \
+             w1.get('number_external'))
+
+        i = 1
+        while i < len(external_wfs):
+            if external_wfs[i].get('number_external') == \
+               external_wfs[i - 1].get('number_external'):
+                external_wfs.pop(i)
+            else:
+                i = i + 1
+        return external_wfs
+
     def get_number_of_amplitudes(self):
         """Gives the total number of amplitudes for this ME"""
 
@@ -1595,9 +2870,15 @@
         """Gives (number or external particles, number of
         incoming particles)"""
 
-        return (len(self.get('processes')[0].get('legs')),
-                len(filter(lambda leg: leg.get('state') == 'initial',
-                           self.get('processes')[0].get('legs'))))
+        external_wfs = filter(lambda wf: wf.get('leg_state') != \
+                              'intermediate',
+                              self.get_all_wavefunctions())
+
+        return (len(set([wf.get('number_external') for wf in \
+                         external_wfs])),
+                len(set([wf.get('number_external') for wf in \
+                         filter(lambda wf: wf.get('leg_state') == False,
+                                external_wfs)])))
 
     def get_helicity_combinations(self):
         """Gives the number of helicity combinations for external
@@ -1609,9 +2890,9 @@
         model = self.get('processes')[0].get('model')
 
         return reduce(lambda x, y: x * y,
-                      [ len(model.get('particle_dict')[leg.get('id')].\
+                      [ len(model.get('particle_dict')[wf.get('pdg_code')].\
                             get_helicity_states())\
-                        for leg in self.get('processes')[0].get('legs') ])
+                        for wf in self.get_external_wavefunctions() ], 1)
 
     def get_helicity_matrix(self):
         """Gives the helicity matrix for external wavefunctions"""
@@ -1623,8 +2904,8 @@
         model = process.get('model')
 
         return apply(itertools.product, [ model.get('particle_dict')[\
-                                         leg.get('id')].get_helicity_states()\
-                                         for leg in process.get('legs') ])
+                                  wf.get('pdg_code')].get_helicity_states()\
+                                  for wf in self.get_external_wavefunctions()])
 
     def get_denominator_factor(self):
         """Calculate the denominator factor due to:
@@ -1633,7 +2914,7 @@
 
         model = self.get('processes')[0].get('model')
 
-        initial_legs = filter(lambda leg: leg.get('state') == 'initial', \
+        initial_legs = filter(lambda leg: leg.get('state') == False, \
                               self.get('processes')[0].get('legs'))
 
         spin_factor = reduce(lambda x, y: x * y,
@@ -1646,21 +2927,244 @@
                                     get('color')\
                                 for leg in initial_legs ])
 
-        final_legs = filter(lambda leg: leg.get('state') == 'final', \
-                              self.get('processes')[0].get('legs'))
-
-        identical_indices = {}
-        for leg in final_legs:
-            try:
-                identical_indices[leg.get('id')] = \
-                                    identical_indices[leg.get('id')] + 1
-            except KeyError:
-                identical_indices[leg.get('id')] = 1
-        identical_factor = reduce(lambda x, y: x * y,
-                                  [ math.factorial(val) for val in \
-                                    identical_indices.values() ])
-
-        return spin_factor * color_factor * identical_factor
+        return spin_factor * color_factor * self['identical_particle_factor']
+
+    def get_color_amplitudes(self):
+        """Return a list of (coefficient, amplitude number) lists,
+        corresponding to the JAMPs for this matrix element. The
+        coefficients are given in the format (fermion factor, color
+        coeff (frac), imaginary, Nc power)."""
+
+        if not self.get('color_basis'):
+            # No color, simply add all amplitudes with correct factor
+            # for first color amplitude
+            col_amp = []
+            for diagram in self.get('diagrams'):
+                for amplitude in diagram.get('amplitudes'):
+                    col_amp.append(((amplitude.get('fermionfactor'),
+                                    1, False, 0),
+                                    amplitude.get('number')))
+            return [col_amp]
+
+        # There is a color basis - create a list of coefficients and
+        # amplitude numbers
+
+        col_amp_list = []
+        for i, col_basis_elem in \
+                enumerate(sorted(self.get('color_basis').keys())):
+
+            col_amp = []
+            for diag_tuple in self.get('color_basis')[col_basis_elem]:
+                res_amp = filter(lambda amp: \
+                          tuple(amp.get('color_indices')) == diag_tuple[1],
+                          self.get('diagrams')[diag_tuple[0]].get('amplitudes'))
+                if not res_amp:
+                    raise FortranWriter.FortranWriterError, \
+                          """No amplitude found for color structure
+                            %s and color index chain (%s) (diagram %i)""" % \
+                            (col_basis_elem,
+                             str(diag_tuple[1]),
+                             diag_tuple[0])
+
+                if len(res_amp) > 1:
+                    raise FortranWriter.FortranWriterError, \
+                        """More than one amplitude found for color structure
+                        %s and color index chain (%s) (diagram %i)""" % \
+                        (col_basis_elem,
+                         str(diag_tuple[1]),
+                         diag_tuple[0])
+
+                col_amp.append(((res_amp[0].get('fermionfactor'),
+                                 diag_tuple[2],
+                                 diag_tuple[3],
+                                 diag_tuple[4]),
+                                res_amp[0].get('number')))
+
+            col_amp_list.append(col_amp)
+
+        return col_amp_list
+
+
+
+    @staticmethod
+    def check_equal_decay_processes(decay1, decay2):
+        """Check if two single-sided decay processes
+        (HelasMatrixElements) are equal.
+
+        Note that this has to be called before any combination of
+        processes has occured.
+        
+        Since a decay processes for a decay chain is always generated
+        such that all final state legs are completely contracted
+        before the initial state leg is included, all the diagrams
+        will have identical wave function, independently of the order
+        of final state particles.
+        
+        Note that we assume that the process definitions have all
+        external particles, corresponding to the external
+        wavefunctions.
+        """
+
+        if len(decay1.get('processes')) != 1 or \
+           len(decay2.get('processes')) != 1:
+            raise HelasMatrixElement.PhysicsObjectError, \
+                  "Can compare only single process HelasMatrixElements"
+
+        if len(filter(lambda leg: leg.get('state') == False, \
+                      decay1.get('processes')[0].get('legs'))) != 1 or \
+           len(filter(lambda leg: leg.get('state') == False, \
+                      decay2.get('processes')[0].get('legs'))) != 1:
+            raise HelasMatrixElement.PhysicsObjectError, \
+                  "Call to check_decay_processes_equal requires " + \
+                  "both processes to be unique"
+
+        # Compare bulk process properties (number of external legs,
+        # identity factors, number of diagrams, number of wavefunctions
+        # initial leg, final state legs
+        if len(decay1.get('processes')[0].get("legs")) != \
+           len(decay1.get('processes')[0].get("legs")) or \
+           len(decay1.get('diagrams')) != len(decay2.get('diagrams')) or \
+           decay1.get('identical_particle_factor') != \
+           decay2.get('identical_particle_factor') or \
+           sum(len(d.get('wavefunctions')) for d in \
+               decay1.get('diagrams')) != \
+           sum(len(d.get('wavefunctions')) for d in \
+               decay2.get('diagrams')) or \
+           decay1.get('processes')[0].get('legs')[0].get('id') != \
+           decay2.get('processes')[0].get('legs')[0].get('id') or \
+           sorted([leg.get('id') for leg in \
+                   decay1.get('processes')[0].get('legs')[1:]]) != \
+           sorted([leg.get('id') for leg in \
+                   decay2.get('processes')[0].get('legs')[1:]]):
+            return False
+
+        # Run a quick check to see if the processes are already
+        # identical (i.e., the wavefunctions are in the same order)
+        if [leg.get('id') for leg in \
+            decay1.get('processes')[0].get('legs')] == \
+           [leg.get('id') for leg in \
+            decay2.get('processes')[0].get('legs')] and \
+            decay1 == decay2:
+            return True
+
+        # Now check if all diagrams are identical. This is done by a
+        # recursive function starting from the last wavefunction
+        # (corresponding to the initial state), since it is the
+        # same steps for each level in mother wavefunctions
+
+        amplitudes2 = copy.copy(reduce(lambda a1, d2: a1 + \
+                                       d2.get('amplitudes'),
+                                       decay2.get('diagrams'), []))
+
+        for amplitude1 in reduce(lambda a1, d2: a1 + d2.get('amplitudes'),
+                                  decay1.get('diagrams'), []):
+            foundamplitude = False
+            for amplitude2 in amplitudes2:
+                if HelasMatrixElement.check_equal_wavefunctions(\
+                   amplitude1.get('mothers')[-1],
+                   amplitude2.get('mothers')[-1]):
+                    foundamplitude = True
+                    # Remove amplitude2, since it has already been matched
+                    amplitudes2.remove(amplitude2)
+                    break
+            if not foundamplitude:
+                return False
+
+        return True
+
+    @staticmethod
+    def check_equal_wavefunctions(wf1, wf2):
+        """Recursive function to check if two wavefunctions are equal.
+        First check that mothers have identical pdg codes, then repeat for
+        all mothers with identical pdg codes."""
+
+        # End recursion with False if the wavefunctions do not have
+        # the same mother pdgs
+        if sorted([wf.get('pdg_code') for wf in wf1.get('mothers')]) != \
+           sorted([wf.get('pdg_code') for wf in wf2.get('mothers')]):
+            return False
+
+        # End recursion with True if these are external wavefunctions
+        # (note that we have already checked that the pdgs are
+        # identical)
+        if not wf1.get('mothers') and not wf2.get('mothers'):
+            return True
+
+        mothers2 = copy.copy(wf2.get('mothers'))
+
+        for mother1 in wf1.get('mothers'):
+            # Compare mother1 with all mothers in wf2 that have not
+            # yet been used and have identical pdg codes
+            equalmothers = filter(lambda wf: wf.get('pdg_code') == \
+                                  mother1.get('pdg_code'),
+                                  mothers2)
+            foundmother = False
+            for mother2 in equalmothers:
+                if HelasMatrixElement.check_equal_wavefunctions(\
+                    mother1, mother2):
+                    foundmother = True
+                    # Remove mother2, since it has already been matched
+                    mothers2.remove(mother2)
+                    break
+            if not foundmother:
+                return False
+
+        return True
+
+    # This gives the order in which the different spin states will be
+    # written in all Helas calls. Note that this is 
+    sort_spin_dict = {1: 1, -2: 4, 2: 3, 3: 2, 5: 0}
+
+    @staticmethod
+    def sorted_mothers(arg):
+        """Gives a list of mother wavefunctions sorted according to
+        1. the spin order needed in the Fortran Helas calls and
+        2. the order of the particles in the interaction (cyclic)"""
+
+        if not isinstance(arg, HelasWavefunction) and \
+               not isinstance(arg, HelasAmplitude):
+            raise base_objects.PhysicsObject.PhysicsObjectError, \
+                  "%s is not a valid HelasWavefunction or HelasAmplitude" % \
+                  repr(arg)
+
+        if not arg.get('interaction_id'):
+            return arg.get('mothers')
+
+        sorted_mothers1 = copy.copy(arg.get('mothers'))
+
+        # Next sort according to interaction pdg codes
+
+        mother_codes = [ wf.get_pdg_code() for wf \
+                         in sorted_mothers1 ]
+        pdg_codes = copy.copy(arg.get('pdg_codes'))
+        if isinstance(arg, HelasWavefunction):
+            my_code = arg.get_anti_pdg_code()
+            # We need to create the cyclic pdg_codes
+            missing_index = pdg_codes.index(my_code)
+            pdg_codes_cycl = pdg_codes[missing_index + 1:] + \
+                             pdg_codes[:missing_index]
+        else:
+            pdg_codes_cycl = pdg_codes
+
+        sorted_mothers2 = HelasWavefunctionList()
+        for code in pdg_codes_cycl:
+            index = mother_codes.index(code)
+            mother_codes.pop(index)
+            sorted_mothers2.append(sorted_mothers1.pop(index))
+
+        if sorted_mothers1:
+            raise base_objects.PhysicsObject.PhysicsObjectError, \
+                  "Mismatch of pdg codes, %s != %s" % \
+                  (repr(mother_codes), repr(pdg_codes_cycl))
+
+        # Next sort according to spin_state_number
+        return HelasWavefunctionList(\
+                  sorted(sorted_mothers2, lambda wf1, wf2: \
+                         HelasMatrixElement.sort_spin_dict[\
+                                          wf2.get_spin_state_number()]\
+                         - HelasMatrixElement.sort_spin_dict[\
+                                          wf1.get_spin_state_number()]))
+
 
 #===============================================================================
 # HelasMatrixElementList
@@ -1675,6 +3179,284 @@
         return isinstance(obj, HelasMatrixElement)
 
 #===============================================================================
+# HelasDecayChainProcess
+#===============================================================================
+class HelasDecayChainProcess(base_objects.PhysicsObject):
+    """HelasDecayChainProcess: If initiated with a DecayChainAmplitude
+    object, generates the HelasMatrixElements for the core process(es)
+    and decay chains. Then call combine_decay_chain_processes in order
+    to generate the matrix elements for all combined processes."""
+
+    def default_setup(self):
+        """Default values for all properties"""
+
+        self['core_processes'] = HelasMatrixElementList()
+        self['decay_chains'] = HelasDecayChainProcessList()
+
+    def filter(self, name, value):
+        """Filter for valid process property values."""
+
+        if name == 'core_processes':
+            if not isinstance(value, HelasMatrixElementList):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid HelasMatrixElementList object" % \
+                        str(value)
+
+        if name == 'decay_chains':
+            if not isinstance(value, HelasDecayChainProcessList):
+                raise self.PhysicsObjectError, \
+                     "%s is not a valid HelasDecayChainProcessList object" % \
+                     str(value)
+
+        return True
+
+    def get_sorted_keys(self):
+        """Return process property names as a nicely sorted list."""
+
+        return ['core_processes', 'decay_chains']
+
+    def __init__(self, argument=None):
+        """Allow initialization with DecayChainAmplitude"""
+
+        if isinstance(argument, diagram_generation.DecayChainAmplitude):
+            super(HelasDecayChainProcess, self).__init__()
+            self.generate_matrix_elements(argument)
+        elif argument:
+            # call the mother routine
+            super(HelasDecayChainProcess, self).__init__(argument)
+        else:
+            # call the mother routine
+            super(HelasDecayChainProcess, self).__init__()
+
+    def generate_matrix_elements(self, dc_amplitude):
+        """Generate the HelasMatrixElements for the core processes and
+        decay processes (separately)"""
+
+        if not isinstance(dc_amplitude, diagram_generation.DecayChainAmplitude):
+            raise base_objects.PhysicsObjectError, \
+                  "%s is not a valid DecayChainAmplitude" % dc_amplitude
+
+        matrix_elements = self['core_processes']
+
+        # Extract the pdg codes of all particles decayed by decay chains
+        # since these should not be combined in a MultiProcess
+        decay_ids = dc_amplitude.get_decay_ids()
+
+        while dc_amplitude.get('amplitudes'):
+            # Pop the amplitude to save memory space
+            amplitude = dc_amplitude.get('amplitudes').pop(0)
+
+            logger.info("Generating Helas calls for %s" % \
+                        amplitude.get('process').nice_string().\
+                                            replace('Process', 'process'))
+            matrix_element = HelasMatrixElement(amplitude,
+                                                decay_ids=decay_ids,
+                                                gen_color=False)
+
+            try:
+                # If an identical matrix element is already in the list,
+                # then simply add this process to the list of
+                # processes for that matrix element
+                other_processes = matrix_elements[\
+                    matrix_elements.index(matrix_element)].get('processes')
+                logger.info("Combining process with %s" % \
+                      other_processes[0].nice_string().replace('Process: ', ''))
+                other_processes.append(amplitude.get('process'))
+            except ValueError:
+                # Otherwise, if the matrix element has any diagrams,
+                # add this matrix element.
+                if matrix_element.get('processes') and \
+                       matrix_element.get('diagrams'):
+                    matrix_elements.append(matrix_element)
+
+        while dc_amplitude.get('decay_chains'):
+            # Pop the amplitude to save memory space
+            decay_chain = dc_amplitude.get('decay_chains').pop(0)
+            self['decay_chains'].append(HelasDecayChainProcess(\
+                decay_chain))
+
+    def combine_decay_chain_processes(self):
+        """Recursive function to generate complete
+        HelasMatrixElements, combining the core process with the decay
+        chains.
+
+        * If there are several identical final state particles and only
+        one decay chain defined, apply this decay chain to all
+        copies.
+        * If there are several decay chains defined for the same
+        particle, apply them in order of the FS particles and the
+        defined decay chains."""
+
+        # End recursion when there are no more decay chains
+        if not self['decay_chains']:
+            # Just return the list of matrix elements
+            return self['core_processes']
+
+        # decay_elements is a list of HelasMatrixElementLists with
+        # all decay processes
+        decay_elements = []
+
+        for decay_chain in self['decay_chains']:
+            # This is where recursion happens
+            decay_elements.append(decay_chain.combine_decay_chain_processes())
+
+        # Store the result in matrix_elements
+        matrix_elements = HelasMatrixElementList()
+
+        # List of list of ids for the initial state legs in all decay
+        # processes
+        decay_is_ids = [[element.get('processes')[0].get_initial_ids()[0] \
+                         for element in elements]
+                         for elements in decay_elements]
+
+        while self['core_processes']:
+            # Pop the process to save memory space
+            core_process = self['core_processes'].pop(0)
+            # Get all final state legs
+            fs_legs = core_process.get('processes')[0].get_final_legs()
+            fs_ids = [leg.get('id') for leg in fs_legs]
+            decay_lists = []
+            # Loop over unique final state particle ids
+            for fs_id in set(fs_ids):
+                # Check if the particle id for this leg has a decay
+                # chain defined
+                if not any([any([id == fs_id for id \
+                            in is_ids]) for is_ids in decay_is_ids]):
+                    continue
+                # decay_list has the leg numbers and decays for this
+                # fs particle id
+                decay_list = []
+                # Now check if the number of decay chains with
+                # this particle id is the same as the number of
+                # identical particles in the core process - if so,
+                # use one chain for each of the identical
+                # particles. Otherwise, use all combinations of
+                # decay chains for the particles.
+
+                # Indices for the decay chain lists which contain at
+                # least one decay for this final state
+                chain_indices = filter(lambda index: fs_id in \
+                                                   decay_is_ids[index],
+                                       range(len(decay_is_ids)))
+
+                my_fs_legs = filter(lambda leg: leg.get('id') == fs_id,
+                                    fs_legs)
+                leg_numbers = [leg.get('number') for leg in my_fs_legs]
+
+                if len(leg_numbers) > 1 and \
+                       len(leg_numbers) == len(chain_indices):
+
+                    # The decay of the different fs parts is given
+                    # by the different decay chains, respectively
+                    # Chains is a list of matrix element lists
+                    chains = []
+                    for index in chain_indices:
+                        decay_chains = decay_elements[index]
+                        chains.append(filter(lambda me: \
+                                             me.get('processes')[0].\
+                                             get_initial_ids()[0] == fs_id,
+                                             decay_chains))
+
+                    # Combine decays for this final state type
+                    for element in itertools.product(*chains):
+                        decay_list.append([[n, d] for [n, d] in \
+                                           zip(leg_numbers, element)])
+                else:
+                    # We let the particles decay according to the
+                    # first decay list only
+                    proc_index = chain_indices[0]
+
+                    # Generate all combinations of decay chains with
+                    # the given decays, without double counting
+                    decay_indices = filter(lambda index: fs_id == \
+                                           decay_is_ids[proc_index][index],
+                                       range(len(decay_is_ids[proc_index])))
+
+
+                    red_decay_ids = []
+                    decay_ids = [decay_indices] * len(leg_numbers)
+                    # Combine all decays for this final state type,
+                    # without double counting
+                    for prod in itertools.product(*decay_ids):
+
+                        # Remove double counting between final states
+                        if tuple(sorted(prod)) in red_decay_ids:
+                            continue
+
+                        # Specify decay processes in the matrix element process
+                        red_decay_ids.append(tuple(sorted(prod)));
+
+                        # Pick out the decays for this iteration
+                        decays = [decay_elements[proc_index][chain_index] \
+                                  for chain_index in prod]
+
+                        decay_list.append([[n, d] for [n, d] in \
+                                           zip(leg_numbers, decays)])
+
+                decay_lists.append(decay_list)
+
+            # Finally combine all decays for this process,
+            # and combine them, decay by decay
+            for decays in itertools.product(*decay_lists):
+
+                # Generate a dictionary from leg number to decay process
+                decay_dict = dict(sum(decays, []))
+
+                # Make sure to not modify the original matrix element
+                matrix_element = copy.deepcopy(core_process)
+                # Avoid Python copying the complete model every time
+                for i, process in enumerate(matrix_element.get('processes')):
+                    process.set('model',
+                            core_process.get('processes')[i].get('model'))
+                # Need to replace Particle in all wavefunctions to avoid
+                # deepcopy
+                idecay = 0
+                org_wfs = core_process.get_all_wavefunctions()
+                for i, wf in enumerate(matrix_element.get_all_wavefunctions()):
+                    wf.set('particle', org_wfs[i].get('particle'))
+                    wf.set('antiparticle', org_wfs[i].get('antiparticle'))
+
+                # Insert the decay chains
+                logger.info("Combine %s with decays %s" % \
+                            (core_process.get('processes')[0].nice_string().\
+                             replace('Process: ', ''), \
+                             ", ".join([d.get('processes')[0].nice_string().\
+                                        replace('Process: ', '') \
+                                        for d in decay_dict.values()])))
+
+                matrix_element.insert_decay_chains(decay_dict)
+
+                try:
+                    # If an identical matrix element is already in the list,
+                    # then simply add this process to the list of
+                    # processes for that matrix element
+                    other_processes = matrix_elements[\
+                    matrix_elements.index(matrix_element)].get('processes')
+                    logger.info("Combining process with %s" % \
+                      other_processes[0].nice_string().replace('Process: ', ''))
+                    other_processes.extend(matrix_element.get('processes'))
+                except ValueError:
+                    # Otherwise, if the matrix element has any diagrams,
+                    # add this matrix element.
+                    if matrix_element.get('processes') and \
+                           matrix_element.get('diagrams'):
+                        matrix_elements.append(matrix_element)
+
+        return matrix_elements
+
+#===============================================================================
+# HelasDecayChainProcessList
+#===============================================================================
+class HelasDecayChainProcessList(base_objects.PhysicsObjectList):
+    """List of HelasDecayChainProcess objects
+    """
+
+    def is_valid_element(self, obj):
+        """Test if object obj is a valid HelasDecayChainProcess for the list."""
+
+        return isinstance(obj, HelasDecayChainProcess)
+
+#===============================================================================
 # HelasMultiProcess
 #===============================================================================
 class HelasMultiProcess(base_objects.PhysicsObject):
@@ -1711,6 +3493,10 @@
         elif isinstance(argument, diagram_generation.MultiProcess):
             super(HelasMultiProcess, self).__init__()
             self.generate_matrix_elements(argument.get('amplitudes'))
+        elif isinstance(argument, diagram_generation.Amplitude):
+            super(HelasMultiProcess, self).__init__()
+            self.generate_matrix_elements(\
+                diagram_generation.AmplitudeList([argument]))
         elif argument:
             # call the mother routine
             super(HelasMultiProcess, self).__init__(argument)
@@ -1724,7 +3510,7 @@
         defined by HelasMatrixElement.__eq__"""
 
         if not isinstance(amplitudes, diagram_generation.AmplitudeList):
-            raise self.HelasMultiProcessError, \
+            raise self.PhysicsObjectError, \
                   "%s is not valid AmplitudeList" % repr(amplitudes)
 
         # Keep track of already generated color objects, to reuse as
@@ -1735,56 +3521,73 @@
 
         matrix_elements = self.get('matrix_elements')
 
-        for amplitude in amplitudes:
-            logger.info("Generating Helas calls for %s" % \
-                         amplitude.get('process').nice_string().replace('Process', 'process'))
-            matrix_element = HelasMatrixElement(amplitude, gen_color=False)
-            try:
-                # If an identical matrix element is already in the list,
-                # then simply add this process to the list of
-                # processes for that matrix element
-                other_processes = matrix_elements[\
+        while amplitudes:
+            # Pop the amplitude to save memory space
+            amplitude = amplitudes.pop(0)
+            if isinstance(amplitude, diagram_generation.DecayChainAmplitude):
+                matrix_element_list = HelasDecayChainProcess(amplitude).\
+                                      combine_decay_chain_processes()
+            else:
+                logger.info("Generating Helas calls for %s" % \
+                         amplitude.get('process').nice_string().\
+                                           replace('Process', 'process'))
+                matrix_element_list = [HelasMatrixElement(amplitude,
+                                                          gen_color=False)]
+            for matrix_element in matrix_element_list:
+                if not isinstance(matrix_element, HelasMatrixElement):
+                    raise self.PhysicsObjectError, \
+                          "Not a HelasMatrixElement: ", matrix_element
+                try:
+                    # If an identical matrix element is already in the list,
+                    # then simply add this process to the list of
+                    # processes for that matrix element
+                    other_processes = matrix_elements[\
                     matrix_elements.index(matrix_element)].get('processes')
-                logger.info("Combining process with %s" % \
-                             other_processes[0].nice_string().replace('Process: ', ''))
-                other_processes.append(amplitude.get('process'))
-
-            except ValueError:
-                # Otherwise, if the matrix element has any diagrams,
-                # add this matrix element.
-                if matrix_element.get('processes') and \
-                       matrix_element.get('diagrams'):
-                    matrix_elements.append(matrix_element)
-
-                # Always create an empty color basis, and the list of raw
-                # colorize objects (before simplification) associated with amplitude
-                col_basis = color_amp.ColorBasis()
-                colorize_obj = col_basis.create_color_dict_list(amplitude)
-
-                try:
-                    # If the color configuration of the ME has already been 
-                    # considered before, recycle the information
-                    col_index = list_colorize.index(colorize_obj)
-                    logger.info(\
-                        "Reusing existing color information for %s" % \
-                        amplitude.get('process').nice_string().replace('Process',
-                                                                   'process'))
+                    logger.info("Combining process with %s" % \
+                      other_processes[0].nice_string().replace('Process: ', ''))
+                    other_processes.extend(matrix_element.get('processes'))
                 except ValueError:
-                    # If not, create color basis and color matrix accordingly
-                    list_colorize.append(colorize_obj)
-                    col_basis.build()
-                    list_color_basis.append(col_basis)
-                    col_matrix = color_amp.ColorMatrix(col_basis)
-                    list_color_matrices.append(col_matrix)
-                    col_index = -1
-                    logger.info(\
-                        "Processing color information for %s" % \
-                        amplitude.get('process').nice_string().replace('Process',
-                                                                   'process'))
+                    # Otherwise, if the matrix element has any diagrams,
+                    # add this matrix element.
+                    if matrix_element.get('processes') and \
+                           matrix_element.get('diagrams'):
+                        matrix_elements.append(matrix_element)
+
+                        # Always create an empty color basis, and the
+                        # list of raw colorize objects (before
+                        # simplification) associated with amplitude
+                        col_basis = color_amp.ColorBasis()
+                        new_amp = matrix_element.get_base_amplitude()
+                        matrix_element.set('base_amplitude', new_amp)
+                        colorize_obj = col_basis.create_color_dict_list(new_amp)
+
+                        try:
+                            # If the color configuration of the ME has
+                            # already been considered before, recycle
+                            # the information
+                            col_index = list_colorize.index(colorize_obj)
+                            logger.info(\
+                              "Reusing existing color information for %s" % \
+                              matrix_element.get('processes')[0].nice_string().\
+                                                 replace('Process', 'process'))
+                        except ValueError:
+                            # If not, create color basis and color
+                            # matrix accordingly
+                            list_colorize.append(colorize_obj)
+                            col_basis.build()
+                            list_color_basis.append(col_basis)
+                            col_matrix = color_amp.ColorMatrix(col_basis)
+                            list_color_matrices.append(col_matrix)
+                            col_index = -1
+                            logger.info(\
+                              "Processing color information for %s" % \
+                              matrix_element.get('processes')[0].nice_string().\
+                                             replace('Process', 'process'))
 
                 matrix_element.set('color_basis', list_color_basis[col_index])
-                matrix_element.set('color_matrix', list_color_matrices[col_index])
-
+                matrix_element.set('color_matrix',
+                                   list_color_matrices[col_index])
+            
 #===============================================================================
 # HelasModel
 #===============================================================================
@@ -1850,10 +3653,11 @@
                   repr(matrix_element)
 
         res = []
-        for n, diagram in enumerate(matrix_element.get('diagrams')):
+        for diagram in matrix_element.get('diagrams'):
             res.extend([ self.get_wavefunction_call(wf) for \
                          wf in diagram.get('wavefunctions') ])
-            res.append("# Amplitude(s) for diagram number %d" % (n + 1))
+            res.append("# Amplitude(s) for diagram number %d" % \
+                       diagram.get('number'))
             for amplitude in diagram.get('amplitudes'):
                 res.append(self.get_amplitude_call(amplitude))
 
@@ -1926,3 +3730,4 @@
             self.set('name', argument.get('name'))
         else:
             super(HelasModel, self).__init__(argument)
+

=== modified file 'madgraph/interface/cmd_interface.py'
--- madgraph/interface/cmd_interface.py	2010-02-12 01:58:30 +0000
+++ madgraph/interface/cmd_interface.py	2010-05-07 07:20:46 +0000
@@ -18,6 +18,7 @@
 """
 
 import cmd
+import copy
 import os
 import subprocess
 import sys
@@ -59,10 +60,17 @@
                       'interactions',
                       'processes',
                       'multiparticles']
+    __add_opts = ['process']
     __save_opts = ['model',
                    'processes']
     __import_formats = ['v4']
-    __export_formats = ['v4standalone', 'v4sa_dirs']
+    __export_formats = ['v4standalone', 'v4sa_dirs', 'v4madevent']
+
+
+    class MadGraphCmdError(Exception):
+        """Exception raised if an error occurs in the execution
+        of command."""
+        pass
 
     def split_arg(self, line):
         """Split a line of arguments"""
@@ -145,6 +153,30 @@
         "*                                                          *\n" + \
         "************************************************************"
 
+    def precmd(self, line):
+        """ force the printing of the line if this is executed with an stdin """
+        if not self.use_rawinput:
+            print line
+
+        return line
+    
+
+    def emptyline(self):
+        """If empty line, do nothing. Default is repeat previous command."""
+
+        pass
+    
+    def default(self, line):
+        """Default action if line is not recognized"""
+
+        if line[0] == "#":
+            # This is a comment - do nothing
+            return
+        else:
+            # Faulty command
+            print "Command \"%s\" not recognized, please try again" % \
+                  line.split()[0]
+
     # Import files
     def do_import(self, line):
         """Import files with external formats"""
@@ -231,18 +263,18 @@
             if self.__curr_model:
                 #save_model.save_model(args[1], self.__curr_model)
                 if save_load_object.save_to_file(args[1], self.__curr_model):
-                    print 'Saved model to file ',args[1]
+                    print 'Saved model to file ', args[1]
             else:
                 print 'No model to save!'
         elif args[0] == 'processes':
             if self.__curr_amps:
                 if save_load_object.save_to_file(args[1], self.__curr_amps):
-                    print 'Saved processes to file ',args[1]
+                    print 'Saved processes to file ', args[1]
             else:
                 print 'No processes to save!'
         else:
             self.help_save()
-                
+
     def do_load(self, line):
         """Load information from file"""
 
@@ -260,7 +292,7 @@
                 print "Loaded model from file in %0.3f s" % \
                       (cpu_time2 - cpu_time1)
             else:
-                print 'Error: Could not load model from file ',args[1]
+                print 'Error: Could not load model from file ', args[1]
         elif args[0] == 'processes':
             self.__curr_amps = save_load_object.load_from_file(args[1])
             if isinstance(self.__curr_amps, diagram_generation.AmplitudeList):
@@ -272,10 +304,10 @@
                                         get('process').get('model')
                     print "Model set from process."
             else:
-                print 'Error: Could not load processes from file ',args[1]
+                print 'Error: Could not load processes from file ', args[1]
         else:
             self.help_save()
-                
+
     def complete_save(self, text, line, begidx, endidx):
         "Complete the save command"
 
@@ -320,10 +352,10 @@
 
         #option
         if len(self.split_arg(line[0:begidx])) >= 2:
-            option=['external=', 'horizontal=', 'add_gap=','max_size=', \
+            option = ['external=', 'horizontal=', 'add_gap=', 'max_size=', \
                                 'contract_non_propagating=']
             return self.list_completion(text, option)
-        
+
     # Display
     def do_display(self, line):
         """Display current internal status"""
@@ -362,9 +394,7 @@
 
         if args[0] == 'processes':
             for amp in self.__curr_amps:
-                print amp.get('process').nice_string()
-                print amp.get('diagrams').nice_string()
-
+                print amp.nice_string()
         if args[0] == 'multiparticles':
             print 'Multiparticle labels:'
             for key in self.__multiparticles:
@@ -404,17 +434,129 @@
             " please create one first!"
             return False
 
+        # Reset Helas matrix elements
+        self.__curr_matrix_elements = helas_objects.HelasMultiProcess()
+
+        try:
+            if line.find(',') == -1:
+                myprocdef = self.extract_process(line)
+            else:
+                myprocdef, line = self.extract_decay_chain_process(line)
+        except self.MadGraphCmdError as error:
+            print "Empty or wrong format process, please try again. Error:\n" \
+                  + str(error)
+            myprocdef = None
+            
+        if myprocdef:
+
+            cpu_time1 = time.time()
+
+            myproc = diagram_generation.MultiProcess(myprocdef)
+
+            self.__curr_amps = myproc.get('amplitudes')
+
+            cpu_time2 = time.time()
+
+            nprocs = len(self.__curr_amps)
+            ndiags = sum([amp.get_number_of_diagrams() for \
+                              amp in self.__curr_amps])
+            print "%i processes with %i diagrams generated in %0.3f s" % \
+                  (nprocs, ndiags, (cpu_time2 - cpu_time1))
+
+        else:
+            print "Empty or wrong format process, please try again."
+
+
+    # Helper functions
+    def extract_decay_chain_process(self, line, level_down = False):
+        """Recursively extract a decay chain process definition from a
+        string. Returns a ProcessDefinition."""
+
+        # Start with process number (identified by "@")
+        proc_number_pattern = re.compile("^(.+)@\s*(\d+)\s*(.*)$")
+        proc_number_re = proc_number_pattern.match(line)
+        proc_number = 0
+        if proc_number_re:
+            proc_number = int(proc_number_re.group(2))
+            line = proc_number_re.group(1) + \
+                   proc_number_re.group(3)
+            print line
+            
+        index_comma = line.find(",")
+        index_par = line.find(")")
+        min_index = index_comma
+        if index_par > -1 and (index_par < min_index or min_index == -1):
+            min_index = index_par
+        
+        if min_index > -1:
+            core_process = self.extract_process(line[:min_index], proc_number)
+        else:
+            core_process = self.extract_process(line, proc_number)
+
+        #level_down = False
+
+        while index_comma > -1:
+            line = line[index_comma + 1:]
+            index_par = line.find(')')
+            if line.lstrip()[0] == '(':
+                # Go down one level in process hierarchy
+                #level_down = True
+                line = line.lstrip()[1:]
+                # This is where recursion happens
+                decay_process, line = \
+                            self.extract_decay_chain_process(line,
+                                                             level_down = True)
+                index_comma = line.find(",")
+                index_par = line.find(')')
+            else:
+                index_comma = line.find(",")
+                min_index = index_comma
+                if index_par > -1 and \
+                       (index_par < min_index or min_index == -1):
+                    min_index = index_par
+                if min_index > -1:
+                    decay_process = self.extract_process(line[:min_index])
+                else:
+                    decay_process = self.extract_process(line)
+
+            core_process.get('decay_chains').append(decay_process)
+
+            if level_down:
+                if index_par == -1:
+                    raise self.MadGraphCmdError,\
+                      "Missing ending parenthesis for decay process"
+
+                if index_par < index_comma:
+                    line = line[index_par + 1:]
+                    level_down = False
+                    break
+
+        if level_down:
+            index_par = line.find(')')
+            if index_par == -1:
+                raise self.MadGraphCmdError,\
+                      "Missing ending parenthesis for decay process"
+            line = line[index_par + 1:]
+            
+        # Return the core process (ends recursion when there are no
+        # more decays)
+        return core_process, line
+        
+    def extract_process(self, line, proc_number = 0):
+        """Extract a process definition from a string. Returns
+        a ProcessDefinition."""
+
         # Use regular expressions to extract s-channel propagators,
         # forbidden s-channel propagators/particles, coupling orders
         # and process number, starting from the back
 
         # Start with process number (identified by "@")
-        proc_number_pattern = re.compile("^(.+)\s*@\s*(\d+)\s*$")
+        proc_number_pattern = re.compile("^(.+)@\s*(\d+)\s*(.*)$")
         proc_number_re = proc_number_pattern.match(line)
-        proc_number = 0
         if proc_number_re:
             proc_number = int(proc_number_re.group(2))
-            line = proc_number_re.group(1)
+            line = proc_number_re.group(1) + \
+                   proc_number_re.group(3)
 
         # Start with coupling orders (identified by "=")
         order_pattern = re.compile("^(.+)\s+(\w+)\s*=\s*(\d+)\s*$")
@@ -429,11 +571,19 @@
         line = line.lower()
 
         # Now check for forbidden particles, specified using "/"
-        forbidden_particles_re = re.match("^(.+)\s*/\s*(.+)\s*$", line)
+        slash = line.find("/")
+        dollar = line.find("$")
         forbidden_particles = ""
-        if forbidden_particles_re:
-            forbidden_particles = forbidden_particles_re.group(2)
-            line = forbidden_particles_re.group(1)
+        if slash > 0:
+            if dollar > slash:
+                forbidden_particles_re = re.match("^(.+)\s*/\s*(.+\s*)(\$.*)$", line)
+            else:
+                forbidden_particles_re = re.match("^(.+)\s*/\s*(.+\s*)$", line)
+            if forbidden_particles_re:
+                forbidden_particles = forbidden_particles_re.group(2)
+                line = forbidden_particles_re.group(1)
+                if len(forbidden_particles_re.groups()) > 2:
+                    line = line + forbidden_particles_re.group(3)
 
         # Now check for forbidden schannels, specified using "$"
         forbidden_schannels_re = re.match("^(.+)\s*\$\s*(.+)\s*$", line)
@@ -452,21 +602,18 @@
 
         args = self.split_arg(line)
 
-        # Reset Helas matrix elements
-        self.__curr_matrix_elements = helas_objects.HelasMultiProcess()
-
         myleglist = base_objects.MultiLegList()
-        state = 'initial'
+        state = False
         number = 1
 
         # Extract process
         for part_name in args:
 
             if part_name == '>':
-                if not len(myleglist):
-                    print "Empty or wrong format process, please try again."
-                    return False
-                state = 'final'
+                if not myleglist:
+                    raise self.MadGraphCmdError, \
+                          "No final state particles"
+                state = True
                 continue
 
             mylegids = []
@@ -481,9 +628,10 @@
                 myleglist.append(base_objects.MultiLeg({'ids':mylegids,
                                                         'state':state}))
             else:
-                print "No particle %s in model: skipped" % part_name
+                raise self.MadGraphCmdError,\
+                      "No particle %s in model" % part_name
 
-        if filter(lambda leg: leg.get('state') == 'final', myleglist):
+        if filter(lambda leg: leg.get('state') == True, myleglist):
             # We have a valid process
 
             # Now extract restrictions
@@ -491,6 +639,10 @@
             forbidden_schannel_ids = []
             required_schannel_ids = []
 
+            decay_process = len(filter(lambda leg: \
+                                       leg.get('state') == False,
+                                       myleglist)) == 1
+
             if forbidden_particles:
                 args = self.split_arg(forbidden_particles)
                 for part_name in args:
@@ -521,37 +673,96 @@
                         if mypart:
                             required_schannel_ids.append(mypart.get_pdg_code())
 
-
-
-            myprocdef = base_objects.ProcessDefinitionList([\
+            return \
                 base_objects.ProcessDefinition({'legs': myleglist,
-                                                'model': self.__curr_model,
-                                                'id': proc_number,
-                                                'orders': orders,
-                                                'forbidden_particles': forbidden_particle_ids,
-                                                'forbidden_s_channels': forbidden_schannel_ids,
-                                                'required_s_channels': required_schannel_ids \
-                                                })])
-            myproc = diagram_generation.MultiProcess({'process_definitions':\
-                                                      myprocdef})
-
-            cpu_time1 = time.time()
-
-            self.__curr_amps = myproc.get('amplitudes')
-
-            cpu_time2 = time.time()
-
-            nprocs = len(filter(lambda amp: amp.get("diagrams"),
-                                self.__curr_amps))
-            ndiags = sum([len(amp.get('diagrams')) for \
-                              amp in self.__curr_amps])
-            print "%i processes with %i diagrams generated in %0.3f s" % \
-                  (nprocs, ndiags, (cpu_time2 - cpu_time1))
-
+                                'model': self.__curr_model,
+                                'id': proc_number,
+                                'orders': orders,
+                                'forbidden_particles': forbidden_particle_ids,
+                                'forbidden_s_channels': forbidden_schannel_ids,
+                                'required_s_channels': required_schannel_ids
+                                 })
+        #                       'is_decay_chain': decay_process\
         else:
-            print "Empty or wrong format process, please try again."
+            return None
 
+    # Add a process to the existing multiprocess definition
     # Generate a new amplitude
+    def do_add(self, line):
+        """Generate an amplitude for a given process and add to
+        existing amplitudes"""
+
+        if len(line) < 1:
+            self.help_add()
+            return False
+
+        if len(self.__curr_model['particles']) == 0:
+            print "No particle list currently active, please create one first!"
+            return False
+
+        if len(self.__curr_model['interactions']) == 0:
+            print "No interaction list currently active," + \
+            " please create one first!"
+            return False
+
+        args = self.split_arg(line)
+        if len(args) < 2:
+            self.help_import()
+            return False
+
+        if args[0] == 'process':
+            # Rejoin line
+            line = ' '.join(args[1:])
+
+            # Reset Helas matrix elements
+            self.__curr_matrix_elements = helas_objects.HelasMultiProcess()
+
+            try:
+                if line.find(',') == -1:
+                    myprocdef = self.extract_process(line)
+                else:
+                    myprocdef, line = self.extract_decay_chain_process(line)
+            except self.MadGraphCmdError as error:
+                print "Empty or wrong format process, please try again. Error:\n" \
+                      + str(error)
+                myprocdef = None
+
+            if myprocdef:
+
+                cpu_time1 = time.time()
+
+                myproc = diagram_generation.MultiProcess(myprocdef)
+
+                for amp in myproc.get('amplitudes'):
+                    if amp not in self.__curr_amps:
+                        self.__curr_amps.append(amp)
+                    else:
+                        print "Warning: Already in processes:"
+                        print amp.nice_string_processes()
+
+                cpu_time2 = time.time()
+
+                nprocs = len(myproc.get('amplitudes'))
+                ndiags = sum([amp.get_number_of_diagrams() for \
+                                  amp in myproc.get('amplitudes')])
+                print "%i processes with %i diagrams generated in %0.3f s" % \
+                      (nprocs, ndiags, (cpu_time2 - cpu_time1))
+                ndiags = sum([amp.get_number_of_diagrams() for \
+                                  amp in self.__curr_amps])
+                print "Total: %i processes with %i diagrams" % \
+                      (len(self.__curr_amps), ndiags)                
+            else:
+                print "Empty or wrong format process, please try again."
+
+
+    def complete_add(self, text, line, begidx, endidx):
+        "Complete the add command"
+
+        # Format
+        if len(self.split_arg(line[0:begidx])) == 1:
+            return self.list_completion(text, self.__add_opts)
+
+    # Export a matrix element
     def do_export(self, line):
         """Export a generated amplitude to file"""
 
@@ -584,7 +795,7 @@
             self.help_export()
             return False
 
-        if not filter(lambda amp: amp.get("diagrams"), self.__curr_amps):
+        if not self.__curr_amps:
             print "No process generated, please generate a process!"
             return False
 
@@ -594,7 +805,7 @@
         ndiags, cpu_time = generate_matrix_elements(self)
         calls = 0
         path = args[1]
-        
+
         if args[0] == 'v4standalone':
             for me in self.__curr_matrix_elements.get('matrix_elements'):
                 filename = os.path.join(path, 'matrix_' + \
@@ -615,13 +826,26 @@
                         export_v4.generate_subprocess_directory_v4_standalone(\
                             me, self.__curr_fortran_model, path)
 
+        if args[0] == 'v4madevent':
+            for me in self.__curr_matrix_elements.get('matrix_elements'):
+                calls = calls + \
+                        export_v4.generate_subprocess_directory_v4_madevent(\
+                            me, self.__curr_fortran_model, path)
+                
         print ("Generated helas calls for %d subprocesses " + \
               "(%d diagrams) in %0.3f s") % \
               (len(self.__curr_matrix_elements.get('matrix_elements')),
                ndiags, cpu_time)
-        
+
         print "Wrote %d helas calls" % calls
 
+        # Replace the amplitudes with the actual amplitudes from the
+        # matrix elements, which allows proper diagram drawing also of
+        # decay chain processes
+        self.__curr_amps = diagram_generation.AmplitudeList(\
+               [me.get('base_amplitude') for me in \
+                self.__curr_matrix_elements.get('matrix_elements')])
+
     def complete_export(self, text, line, begidx, endidx):
         "Complete the export command"
 
@@ -683,7 +907,7 @@
             self.help_draw()
             return False
         
-        if not filter(lambda amp: amp.get("diagrams"), self.__curr_amps):
+        if not self.__curr_amps:
             print "No process generated, please generate a process!"
             return False
 
@@ -700,9 +924,15 @@
                     print "invalid syntax: '%s'. Please try again" % data
                     self.help_draw()
                     return False
-                option.set(key,value)
+                option.set(key, value)
+
+        # Collect amplitudes
+        amplitudes = diagram_generation.AmplitudeList()
 
         for amp in self.__curr_amps:
+            amplitudes.extend(amp.get_amplitudes())            
+
+        for amp in amplitudes:
             filename = os.path.join(args[0], 'diagrams_' + \
                                     amp.get('process').shell_string() + ".eps")
             plot = draw.MultiEpsDiagramDrawer(amp['diagrams'],
@@ -745,7 +975,21 @@
 
         print "syntax: generate INITIAL STATE > REQ S-CHANNEL > FINAL STATE $ EXCL S-CHANNEL / FORBIDDEN PARTICLES COUP1=ORDER1 COUP2=ORDER2"
         print "-- generate diagrams for a given process"
-        print "   Example: u d~ > w+ > m+ vm g $ a / z h QED=3 QCD=0"
+        print "   Example: u d~ > w+ > m+ vm g $ a / z h QED=3 QCD=0 @1"
+        print "Decay chain syntax:"
+        print "   core process, decay1, (decay2, (decay3, ...)), ...  etc"
+        print "   Example: g g > t~ t @2, (t~ > W- b~, W- > e- ve~), t > W+ b"
+        print "   Note that identical particles will all be decayed"
+
+    def help_add(self):
+
+        print "syntax: add process INITIAL STATE > REQ S-CHANNEL > FINAL STATE $ EXCL S-CHANNEL / FORBIDDEN PARTICLES COUP1=ORDER1 COUP2=ORDER2"
+        print "-- generate diagrams for a process and add to existing processes"
+        print "   Syntax example: u d~ > w+ > m+ vm g $ a / z h QED=3 QCD=0 @1"
+        print "Decay chain syntax:"
+        print "   core process, decay1, (decay2, (decay3, ...)), ...  etc"
+        print "   Example: g g > t~ t @2, (t~ > W- b~, W- > e- ve~), t > W+ b"
+        print "   Note that identical particles will all be decayed"
 
     def help_define(self):
         print "syntax: define multipart_name [ part_name_list ]"
@@ -757,7 +1001,11 @@
               " FILEPATH"
         print """-- export matrix elements. For v4standalone, the resulting
         files will be FILEPATH/matrix_\"process_string\".f. For v4sa_dirs,
-        the result is a set of complete MG4 Standalone process directories."""
+        the result is a set of complete MG4 Standalone process directories.
+        For v4madevent, the path needs to be to a MadEvent SubProcesses
+        directory, and the result is the Pxxx directories (including the
+        diagram .ps and .jpg files) for the subprocesses as well as a
+        correctly generated subproc.mg file."""
 
     def help_draw(self):
         print "syntax: draw FILEPATH [option=value]"

=== modified file 'madgraph/iolibs/drawing.py'
--- madgraph/iolibs/drawing.py	2010-01-29 06:40:30 +0000
+++ madgraph/iolibs/drawing.py	2010-05-07 07:20:46 +0000
@@ -86,12 +86,12 @@
         return
 
     def def_end_point(self, vertex):
-        """-Re-Define the starting point of the line."""
+        """-Re-Define the starting point of the line. with check"""
 
         if not isinstance(vertex, VertexPoint):
             raise self.FeynmanLineError, 'The end point should be a ' + \
                  'Vertex_Point object'
-
+                 
         self.end = vertex
         vertex.add_line(self)
         return
@@ -741,7 +741,7 @@
                 vertex_point = VertexPoint(vertex)
                 self.vertexList.append(vertex_point)
                 # If initial state particle, we will need to flip begin-end
-                if line.get('state') == 'initial':
+                if line.get('state') == False:
                     if line.start:
                         line.inverse_begin_end()
                     line.def_begin_point(vertex_point)
@@ -797,7 +797,7 @@
         2) Add this vertex in vertexList of the diagram
         3) Update vertex.line list. (first update the leg into line if needed)
         4) assign line.start[end] to this vertex. (in end if start is already
-                assigned to another vertex). the start-end will be flipped later
+                assigned to another vertex). the start-end will be flip later
                 if needed."""
 
         #1) Extend to a vertex point
@@ -875,20 +875,32 @@
         #consequence we replace in line2 the common vertex by the second vertex 
         #present in line1, such that the new line2 is the sum of the two 
         #previous lines.
+        to_del = pos1
         if line1.end is line2.start:
             line2.def_begin_point(line1.start)
         else:
-            line2.def_end_point(line1.end)
+            if line1.end:
+                # Standard case
+                line2.def_end_point(line1.end)
+            else:
+                # line1 is a initial/final line. We need to keep her status
+                #so we will delete line2 instead of line1
+                to_del = pos2
+                line1.def_begin_point(line2.end)
         # Remove line completely
-        line1.start.remove_line(line1)
-        del self.lineList[pos1]
+        if to_del == pos1:
+            line1.start.remove_line(line1)
+            del self.lineList[pos1]
+        else:
+            line2.start.remove_line(line2)
+            del self.lineList[pos2]
         # Remove last_vertex
         self.vertexList.remove(last_vertex)
 
     def deal_last_line(self, last_line):
         """The line of the last vertex breaks the rules that line before
         '>' exist previously and the one after don't. The last one can also
-        already exist and for the one befor the '>' sometimes they arrive 
+        already exist and for the one before the '>' sometimes they arrive 
         with a second object which is equivalent to another one but not 
         the same object. discover those case and treat this properly."""
 
@@ -898,7 +910,6 @@
             id1 = self.find_leg_id(last_line)
             # Find if they are a second call to this line
             id2 = self.find_leg_id(last_line, end=len(self._treated_legs) - id1)
-
             if id2 is not None:
                 # Line is duplicate in linelist => remove this duplication
                 line = self.lineList[id2]
@@ -920,7 +931,7 @@
                     else:
                         line.def_start_point(last_line.start)
                     # Remove last_line from the vertex    
-                    last_line.begin.remove_line(last_line)
+                    last_line.start.remove_line(last_line)
 
                 # Remove last_line
                 self.lineList.remove(last_line)
@@ -977,7 +988,7 @@
                 continue
             # Check if T-channel or not. Note that T-channel tag is wrongly 
             #define if only one particle in initial state.
-            if line.get('state') == 'initial':
+            if line.get('state') == False:
                 # This is T vertex. => level is 1
                 line.end.def_level(1)
             else:
@@ -1016,7 +1027,7 @@
         evolution direction (which will be the wrong one at the next step)."""
 
         for line in t_vertex.line:
-            if line.get('state') == 'initial' and line.start is t_vertex:
+            if line.get('state') == False and line.start is t_vertex:
                 return line.end
 
     def find_vertex_at_level(self, previous_level):
@@ -1035,7 +1046,7 @@
                 continue
 
             for line in vertex.line:
-                if line.start is vertex and line.get('state') == 'final':
+                if line.start is vertex and line.get('state') == True:
                         vertex_at_level.append(line.end)
 
         return vertex_at_level
@@ -1233,8 +1244,8 @@
         This occur for 1>X diagram where T-channel are wrongly define."""
         
         for line in self.lineList:
-            if line.get('state') == 'initial':
-                line.set('state','final')
+            if line.get('state') == False:
+                line.set('state', True)
 
 
     def solve_line_direction(self):
@@ -1246,7 +1257,7 @@
 
         # Use the basic rules. Assigns correctly but for T-channel
         for line in self.lineList:
-            if line.get('state') == 'final':
+            if line.get('state') == True:
                 line.define_line_orientation()
         # The define line orientation use level information and in consequence 
         #fails on T-Channel. So in consequence we still have to fix T-channel
@@ -1269,7 +1280,7 @@
             for line in t_vertex.line:
 
                 # Identify the next T-channel particles
-                if line.get('state') == 'initial' and t_old is not line and \
+                if line.get('state') == False and t_old is not line and \
                     line.start is t_vertex:
                     t_next = line
                     
@@ -1320,7 +1331,7 @@
             if line.is_external():
                 # Check the size of final particles to restrict to the max_size
                 #constraints.
-                if line.get('state') == 'initial' or not line.is_external():
+                if line.get('state') == False or not line.is_external():
                     continue 
                 size = line.get_length() * self.max_level
                 if size > finalsize:
@@ -1874,4 +1885,4 @@
         """Convert the value in a number"""
         
         return float(value)
-        
\ No newline at end of file
+        

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2010-02-12 12:53:59 +0000
+++ madgraph/iolibs/export_v4.py	2010-05-07 07:20:46 +0000
@@ -15,17 +15,18 @@