~sebastien-wertz/mg5amcnlo/standalone_cpp_mem

################################################################################
#
# Copyright (c) 2009 The MadGraph Development team and Contributors
#
# This file is a part of the MadGraph 5 project, an application which 
# automatically generates Feynman diagrams and matrix elements for arbitrary
# high-energy processes in the Standard Model and beyond.
#
# It is subject to the MadGraph license which should accompany this 
# distribution.
#
# For more information, please visit: http://madgraph.phys.ucl.ac.be
#
################################################################################
"""Definition for the objects used in the decay module.
   DecayParticle: this object contains all the decay related properties
                  including all the decay vertices and decay channels.
                  This object also has the 'is_stable' label to denote
                  wether this particle is stable.
   DecayParticleList: this helper class will help to turn Particle 
                      in base_objects into DecayParticle.
   DecayModel: This model contains DecayParticle. Particle in base_objects
               will be automatically converted into DecayParticle either 
               during the initialization or when a ParticleList is set as
               the particles of DecayModel through set function.
               This model can search all the decay_vertexlist for all
               its particles at one time. The function 'find_stable_particles'
               and 'find_decay_groups_general' will find stable particles
               automatically from interaction and mass without any input.
   Channel: A diagram object specialized for decay process. This includes
            several helper functions for channel generations and the functions
            to calculate the approximate decay width.
   ChannelList: A list of channels.
   
   Users can run DecayParticle.find_channels(final_state_number, model)
   to get the channels with final state number they request. Or they
   may run DecayModel.find_all_channels(final_state_number) to search
   the channels for all particles in the given model."""
   
from __future__ import division

import array
import cmath
import copy
import itertools
import logging
import math
import os
import re
import sys
import time

import madgraph.core.base_objects as base_objects
import madgraph.core.diagram_generation as diagram_generation
import madgraph.core.color_amp as color_amp
import madgraph.core.color_algebra as color
import madgraph.core.helas_objects as helas_objects
import madgraph.core.diagram_generation as diagram_generation

import madgraph.various.misc as misc

import models.import_ufo as import_ufo
from madgraph import MadGraph5Error, MG5DIR

import models.model_reader as model_reader

ZERO = 0
#===============================================================================
# Logger for decay_module
#===============================================================================

logger = logging.getLogger('madevent.decay_objects')


#===============================================================================
# DecayParticle
#===============================================================================
class DecayParticle(base_objects.Particle):
    """DecayParticle includes the decay vertices and channels for this
       particle. The label 'is_stable' specifies whether the particle
       is stable. The function find_channels will launch all necessary
       functions it needs. If users try to find channels with more final
       state particles, they can run find_channels_nextlevel to get the
       channels in next level.
    """
    sorted_keys = ['name', 'antiname', 'spin', 'color',
                   'charge', 'mass', 'width', 'pdg_code',
                   'line',
                   'is_part', 'self_antipart', 'is_stable',
                   'decay_vertexlist', 'decay_channels', 'apx_decaywidth',
                   'apx_decaywidth_err', '2body_massdiff',
                   'radiations'
                  ]


    def __init__(self, init_dict={}, force=False):
        """Creates a new particle object. If a dictionary is given, tries to 
        use it to give values to properties.
        A repeated assignment is to avoid error of inconsistent pdg_code and
        initial particle id of vertex"""

        dict.__init__(self)
        self.default_setup()

        assert isinstance(init_dict, dict), \
                            "Argument %s is not a dictionary" % repr(init_dict)

        #To avoid the pdg_code remain 0 and then induce the error when
        #set the vertexlist
        try:
            pid = init_dict['pdg_code']
            self.set('pdg_code', pid)
        except KeyError:
            pass

        for item in init_dict.keys():
            self.set(item, init_dict[item], force)
            
            
    def default_setup(self):
        """Default values for all properties"""
        
        super(DecayParticle, self).default_setup()

        self['is_stable'] = False
        #log of the find_vertexlist history
        self['vertexlist_found'] = False
        self['max_vertexorder'] = 0

        # The decay_vertexlist is a dictionary with vertex list as items and
        # final state particles and on shell condition as keys.
        # decay_channels corresponds to one diagram for each channel, 
        # while decay_amplitudes are a series of diagrams with the same
        # initial and final states.

        self['decay_vertexlist'] = {}
        self['decay_channels'] = {}
        self.decay_amplitudes = {}
        self['apx_decaywidth'] = 0.
        self['apx_decaywidth_err'] = 0.
        self['2body_massdiff'] = 0.
        self['radiations']=[]
        
        
    def get(self, name):
        """ Evaluate some special properties first if the user request. """

        if name == 'apx_decaywidth' \
                and not self[name] \
                and not self['is_stable']:
            self.update_decay_attributes(True, False, True)
            return self[name]
        elif name == 'apx_decaywidth_err' and not self[name]:
            self.estimate_width_error()
            return self[name]
        elif name == 'decay_amplitudes':
            return self.decay_amplitudes
        else:
            # call the mother routine
            return DecayParticle.__bases__[0].get(self, name)

    def set(self, name, value, *args,**opts):
        """assign a value"""

        if name == 'decay_amplitudes':
            self.decay_amplitudes = value
        else:
            out = super(DecayParticle, self).set(name, value, *args, **opts)
            return out
    def check_vertex_condition(self, partnum, onshell, 
                              value = base_objects.VertexList(), model = {}):
        """Check the validity of decay condition, including,
           partnum: final state particle number,
           onshell: on-shell condition,
           value  : the assign vertexlist
           model  : the specific model"""

        #Check if partnum is an integer.
        if not isinstance(partnum, int):
            raise self.PhysicsObjectError, \
                "Final particle number %s must be an integer." % str(partnum)

        #Check if onshell condition is Boolean number.
        if not isinstance(onshell, bool):
            raise self.PhysicsObjectError, \
                "%s must be a Boolean number" % str(onshell)
                
        #Check if the value is a Vertexlist(in base_objects) or a list of vertex
        if not isinstance(value, base_objects.VertexList):
            raise self.PhysicsObjectError, \
                "%s must be VertexList type." % str(value)
                    
        #Check if the model is a valid object.
        if not (isinstance(model, base_objects.Model) or model == {}):
            raise self.PhysicsObjectError, \
                "%s must be a Model" % str(model)
        elif model:
            #Check if the mother particle is in the 'model'
            if not (self.get_pdg_code() in model.get('particle_dict').keys()):
                raise self.PhysicsObjectError, \
                    "The model, %s, does not contain particle %s." \
                    %(model.get('name'), self.get_name())

                            
    def check_vertexlist(self, partnum, onshell, value, model = {}):
        """Check if the all the vertex in the vertexlist satisfy the following
           conditions. If so, return true; if not, raise error messages.

           1. There is an appropriate leg for initial particle.
           2. The number of final particles equals to partnum.
           3. If model is not None, check the onshell condition and
              the initial particle id is the same as calling particle.
        """
        #Check the validity of arguments first
        self.check_vertex_condition(partnum, onshell, value, model)
       
        #Determine the number of final particles.
        #Find all the possible initial particle(s).
        #Check onshell condition if the model is given.
        if model:
            if (abs(eval(self.get('mass')) == 0.)) and (len(value) != 0):
                raise self.PhysicsObjectError, \
                    "Massless particle %s cannot decay." % self['name']

        for vert in value:
            # Reset the number of initial/final particles,
            # initial particle id, and total and initial mass
            num_ini = 0
            radiation = False
            num_final = 0
                
            if model:
                # Calculate the total mass
                total_mass = sum([abs(eval(model.get_particle(l['id']).get('mass'))) for l in vert['legs']])
                ini_mass = abs(eval(self.get('mass')))
                
                # Check the onshell condition
                if (ini_mass.real > (total_mass.real - ini_mass.real))!=onshell:
                    raise self.PhysicsObjectError, \
                        "The on-shell condition is not satisfied."

            for leg in vert.get('legs'):
                # Check if all legs are label by true
                if not leg.get('state'):
                    raise self.PhysicsObjectError, \
                        "The state of leg should all be true"

                # Identify the initial particle
                if leg.get('id') == self.get_pdg_code():
                    # Double anti particle is also radiation
                    if num_ini == 1:
                        radiation = True
                    num_ini = 1
                elif leg.get('id') == self.get_anti_pdg_code() and \
                        not self.get('self_antipart'):
                    radiation = True            

            # Calculate the final particle number
            num_final = len(vert.get('legs'))-num_ini

            # Check the number of final particles is the same as partnum
            if num_final != partnum:
                raise self.PhysicsObjectError, \
                    "The vertex is a %s -body decay, not a %s -body one."\
                    % (str(num_final), str(partnum))

            # Check if there is any appropriate leg as initial particle.
            if num_ini == 0:
                raise self.PhysicsObjectError, \
                    "There is no leg satisfied the mother particle %s"\
                    % str(self.get_pdg_code())

            # Check if the vertex is radiation
            if radiation:
                raise self.PhysicsObjectError, \
                    "The vertex is radiactive for mother particle %s"\
                    % str(self.get_pdg_code())

        return True

    def check_channels(self, partnum, onshell, value = [], model = {}):
        """Check the validity of decay channel condition, including,
           partnum: final state particle number,
           onshell: on-shell condition,
           value  : the assign channel list, all channels in it must
                    be consistent to the given partnum and onshell.
           model  : the specific model."""

        # Check if partnum is an integer.
        if not isinstance(partnum, int):
            raise self.PhysicsObjectError, \
                "Final particle number %s must be an integer." % str(partnum)
        
        # Check if onshell condition is Boolean number.
        if not isinstance(onshell, bool):
            raise self.PhysicsObjectError, \
                "%s must be a Boolean number" % str(onshell)
                
        # Check if the value is a ChannelList
        if (not isinstance(value, ChannelList) and value):
            raise self.PhysicsObjectError, \
                "%s must be ChannelList type." % str(value)
                

        # Check if the partnum is correct for all channels in value
        if any(ch for ch in value if \
                   len(ch.get_final_legs()) != partnum):
            raise self.PhysicsObjectError, \
                "The final particle number of channel should be %d."\
                % partnum
        
        # Check if the initial particle in all channels are as self.
        if any(ch for ch in value if \
                   abs(ch.get_anti_initial_id()) != abs(self.get('pdg_code'))):
            raise self.PhysicsObjectError, \
                "The initial particle is not %d or its antipart." \
                % self.get('pdg_code')

        # Check if the onshell condition is right
        if not (isinstance(model, base_objects.Model) or model == {}):
            raise self.PhysicsObjectError, \
                "%s must be a Model" % str(model)
        elif model:
            # Check if the mother particle is in the 'model'
            if not (self.get_pdg_code() in model.get('particle_dict').keys()):
                raise self.PhysicsObjectError, \
                    "The model, %s, does not contain particle %s." \
                    %(model.get('name'), self.get_name())
            if any([ch for ch in value if onshell != ch.get_onshell(model)]):
                raise self.PhysicsObjectError, \
                    "The onshell condition is not consistent with the model."
        return True


    def check_amplitudes(self, partnum, value):
        """Check the validity of DecayAmplitudes condition, including,
           partnum: final state particle number,
           value  : the assign amplitudelist, all amplitudes in it must
                    be consistent to the processess they hold.
        """
        # Check if partnum is an integer.
        if not isinstance(partnum, int):
            raise self.PhysicsObjectError, \
                "Final particle number %s must be an integer." % str(partnum)
        
        # Check if the value is a DecayAmplitudeList
        if (not isinstance(value, DecayAmplitudeList) and value):
            raise self.PhysicsObjectError, \
                "%s must be DecayAmplitudeList type." % str(value)
                
        return True

    def filter(self, name, value):
        """Filter for valid DecayParticle vertexlist."""
        
        if name == 'decay_vertexlist' or name == 'decay_channels':
            #Value must be a dictionary.
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                    "Decay_vertexlist or decay_channels %s must be a dictionary." % str(value)

            # key must be two element tuple
            for key, item in value.items():
                if not isinstance(key, tuple):
                    raise self.PhysicsObjectError,\
                        "Key %s must be a tuple." % str(key)
                
                if len(key) != 2:
                    raise self.PhysicsObjectError,\
                        "Key %s must have two elements." % str(key)
                
                if name == 'decay_vertexlist':
                    self.check_vertexlist(key[0], key[1], item)
                if name == 'decay_channels':
                    self.check_channels(key[0], key[1], item)          

        if name == 'decay_amplitudes':

            #Value must be a dictionary.
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                    "Decay_amplitudes %s must be a dictionary." % str(value)

            # For each key and item, check them with check_amplitudes
            for key, item in value.items():                
                self.check_amplitudes(key, item)
                    
        if name == 'vertexlist_found' or name == 'is_stable':
            if not isinstance(value, bool):
                raise self.PhysicsObjectError, \
                    "Propery %s should be Boolean type." % name

        if name == 'max_vertexorder':
            if not isinstance(value, int):
                raise self.PhysicsObjectError, \
                    "Property %s should be int type." % name

        # Check apx_decaywidth and apx_decaywidth_err
        if name == 'apx_decaywidth' or name == 'apx_decaywidth_err' \
                or name == '2body_massdiff':
            if not isinstance(value, float) and not isinstance(value, int):
                raise self.PhysicsObjectError, \
                    "Property %s must be float type." % str(value)

        super(DecayParticle, self).filter(name, value)

        return True

    def reset_decay_attributes(self, reset_width, reset_err, reset_br):
        """ Depend on the given arguments, 
            reset the apx_decaywidth, apx_decaywidth_err, and
            branching ratio of the amplitudes in this particle.
            It is necessary when the channels are changed, 
            e.g. find next level."""

        # Reset decay width
        if reset_width:
            self['apx_decaywidth'] = 0.

        # Reset err
        if reset_err:
            self['apx_decaywidth_err'] = 0.
        
        # Reset the branching ratio inside amplitudes
        if reset_br:
            for n, amplist in self.decay_amplitudes.items():
                for amp in amplist:
                    amp.reset_width_br()

    def update_decay_attributes(self, reset_width, reset_err, reset_br, model=None):
        """ This function will update the attributes related to decay width,
            including total width, branching ratios (of each amplitudes),
            and error of width.
            The arguments specify which attributes needs to be updated.
            Note that the width of amplitudes will not be changed.
            If the apx_decaywidth_err needs to be updated, the model
            must be provided!. """
        
        # Reset the related properties
        self.reset_decay_attributes(reset_width, reset_err, reset_br)

        # Update the total width first.
        # (So the decaywidth_err and branching ratios can be calculated with
        # the new width.)
        if reset_width:
            for n, amplist in self.decay_amplitudes.items():
                for amp in amplist:
                    # Do not calculate the br in this moment
                    self['apx_decaywidth'] += amp.get('apx_decaywidth')

        # Update the apx_decaywidth_err
        if reset_err:
            self.estimate_width_error(model)

        # Update the branching ratio in the end with the updated total width
        if reset_br:
            for n, amplist in self.decay_amplitudes.items():
                for amp in amplist:
                    # Reset the br first, so the get function will recalculate
                    # br automatically.
                    amp.get('apx_br')


    def estimate_width_error(self, model=None):
        """ This function will estimate the width error from the highest order
            off-shell channels.
            If model is provided, the apx_decaywidth_err of each channel will
            be (re)calculated. """

        if self['apx_decaywidth']:
            final_level = max([k[0] for k, i in self['decay_channels'].items()])

            # Do not recalculate the apx_decaywidth_nextlevel here
            err = sum([c.get('apx_decaywidth_nextlevel', model) \
                           for c in self.get_channels(final_level, False)])/ \
                           self['apx_decaywidth']

        elif self.get('is_stable'):
            err = 0.
        else:
            err = 1.

        self['apx_decaywidth_err'] = err

        return err

    def decaytable_string(self, format='normal'):
        """ Output the string for the decay table.
            If format is 'full', all the channels in the process will be
            printed."""        

        seperator = str('#'*80)
        output = '\n'+seperator
        output += str('\n#\n#\tPDG\t\tWIDTH\t\tERROR\n')
        output += str('DECAY\t%8d\t%.5e     %.3e  #%s decay\n') \
            %(self.get('pdg_code'), 
              self.get('apx_decaywidth'),
              self.get('apx_decaywidth_err'),
              self.get('name'))
        output += seperator
        # Write the decay table from 2-body decay.
        n = 2
        while n:
            if n in self.get('decay_amplitudes').keys():
                # Do not print empty amplitudes
                if len(self.get_amplitudes(n)):
                    # Titie line
                    output += '\n#\tBR\tNDA       '
                    # ID (ID1, ID2, ID3...)
                    output += '        '.join(['ID%d' %(i+1) for i in range(n)])
                    if format == 'cmp':
                        output += '\tratio'
                    # main content
                    output += '\n%s' \
                        % self.get_amplitudes(n).decaytable_string(format)

                # Increase n to nextlevel    
                n += 1
            else:
                break

        return output+'\n#\n#'
        

    def get_vertexlist(self, partnum ,onshell):
        """Return the n-body on/offshell decay vertexlist. 
           If there is no such vertex, return None.
        
           partnum = number of final particles.
           If onshell=false, return the on-shell list and vice versa.
        """
        #check the validity of arguments
        self.check_vertex_condition(partnum, onshell)
        
        return self.get('decay_vertexlist').get((partnum, onshell), [])

    def set_vertexlist(self, partnum ,onshell, value, model = {}):
        """Set the n_body_decay_vertexlist,
           partnum: n, 
           onshell: True for on-shell decay, and False for off-shell
           value: the decay_vertexlist that is tried to assign.
           model: the underlying model for vertexlist
                  Use to check the correctness of on-shell condition.
        """
        #Check the vertexlist by check_vertexlist
        #Error is raised (by check_vertexlist) if value is not valid
        if self.check_vertexlist(partnum, onshell, value, model):
            self['decay_vertexlist'][(partnum, onshell)] = value

    def get_max_vertexorder(self):
        """ Get the max vertex order of this particle"""
        # Do not include keys without vertexlist in it
        # Both onshell and offshell are consider
        if not self['vertexlist_found']:
            logger.warning("The vertexlist of this particle has not been searched."+"Try find_vertexlist first.")
            return

        vertnum_list = [k[0] for k in self['decay_vertexlist'].keys() \
             if self['decay_vertexlist'][k]]
        if vertnum_list:
            self['max_vertexorder'] = max(vertnum_list)
        else:
            self['max_vertexorder'] = 0

        return self['max_vertexorder']

    # OBSOLETE function. It is recommended to run the find_vertexlist in
    # DecayModel object.
    def find_vertexlist(self, model, option=False):
        """Find the possible decay channel to decay,
           for both on-shell and off-shell.
           If option=False (default), 
           do not rewrite the VertexList if it exists.
           If option=True, rewrite the VertexList anyway.
        """
        
        #Raise error if self is not in model.
        if not (self.get_pdg_code() in model.get('particle_dict').keys()):
            raise self.PhysicsObjectError, \
                    "The parent particle %s is not in the model %s." \
                        % (self.get('name'), model.get('name'))

        #Raise error if option is not Boolean value
        if not isinstance(option, bool):
            raise self.PhysicsObjectError, \
                    "The option %s must be True or False." % str(option)
        
        #If 'vertexlist_found' is true and option is false,
        #no action is proceed.
        if self['vertexlist_found'] and not option:
            return 'The vertexlist has been setup.', \
                'No action proceeds because of False option.'

        # Reset the decay vertex before finding
        self['decay_vertexlist'] = {(2, False) : base_objects.VertexList(),
                                    (2, True)  : base_objects.VertexList(),
                                    (3, False) : base_objects.VertexList(),
                                    (3, True)  : base_objects.VertexList()}

        #Set the vertexlist_found at the end
        self['vertexlist_found'] = True

        # Do not include the massless and stable particle
        model.get('stable_particles')
        if self.get('is_stable'):
            return

        #Go through each interaction...
        for temp_int in model.get('interactions'):
            #Save the particle dictionary (pdg_code & anti_pdg_code to particle)
            partlist = temp_int.get('particles')

            #The final particle number = total particle -1
            partnum = len(partlist)-1
            #Allow only 2 and 3 body decay
            if partnum > 3:
                continue

            #Check if the interaction contains mother particle
            if model.get_particle(self.get_anti_pdg_code()) in partlist:
                #Exclude radiation
                part_id_list = [p.get('pdg_code') for p in partlist]
                if (part_id_list.count(self.get('pdg_code'))) > 1:
                    continue

                total_mass = 0
                ini_mass = abs(eval(self.get('mass')))
                vert = base_objects.Vertex()
                legs = base_objects.LegList()

                # Setup all the legs and find final_mass
                for part in partlist:
                    legs.append(base_objects.Leg({'id': part.get_pdg_code()}))
                    total_mass += abs(eval(part.get('mass')))
                    #Initial particle has not been found: ini_found = True
                    if (part == model.get_particle(self.get_anti_pdg_code())):
                        ini_leg = legs.pop()
                        ini_leg.set('id', self.get_pdg_code())
                    
                #Sort the outgoing leglist for comparison sake (removable!)
                legs.sort(key= lambda l: l['id'], reverse=True)
                # Append the initial leg
                legs.append(ini_leg)

                vert.set('id', temp_int.get('id'))
                vert.set('legs', legs)
                temp_vertlist = base_objects.VertexList([vert])

                #Check validity of vertex (removable)
                """self.check_vertexlist(partnum,
                ini_mass > final_mass,
                temp_vertlist, model)"""

                #Append current vert to vertexlist
                try:
                    self['decay_vertexlist'][(partnum, \
                                            ini_mass > (total_mass-ini_mass))].\
                                            append(vert)
                except KeyError:
                    self['decay_vertexlist'][(partnum, \
                                            ini_mass > (total_mass-ini_mass))] \
                                            = base_objects.VertexList([vert])

        

    def get_channels(self, partnum ,onshell):
        """Return the n-body decay channels.
           partnum = n.
           If onshell=false, return the on-shell list and vice versa.
        """
        #check the validity of arguments
        self.check_channels(partnum, onshell)
        return self.get('decay_channels')[(partnum, onshell)]

    def set_channels(self, partnum ,onshell, value, model = {}):
        """Set the n_body_decay_vertexlist, value is overloading to both
           ChannelList and list of Channels (auto-transformation will proceed)
           partnum: n, 
           onshell: True for on-shell decay, and False for off-shell
           value: the decay_vertexlist that is tried to assign.
           model: the underlying model for vertexlist
                  Use to check the correctness of on-shell condition.
        """
        #Check the vertexlist by check_vertexlist
        #Error is raised (by check_vertexlist) if value is not valid
        if isinstance(value, ChannelList):
            if self.check_channels(partnum, onshell, value, model):
                self['decay_channels'][(partnum, onshell)] = value
        elif isinstance(value, list) and \
                all([isinstance(c, Channel) for c in value]):
            value_transform = ChannelList(value)
            if self.check_channels(partnum, onshell, value_transform, model):
                self['decay_channels'][(partnum, onshell)] = value_transform
        else:
            raise self.PhysicsObjectError, \
                "The input must be a list of diagrams."

    def get_max_level(self):
        """ Get the max channel level that the particle have so far. """
        # Turn off the logger in get_amplitude temporarily
        
        # Initial value
        n = 2
        # Look at the amplitudes or channels to find the max_level
        while self.get_amplitudes(n) or ((n,False) in self['decay_channels'].keys()):
            n += 1

        # n is the failed value, return n-1.
        return (n-1)

    def get_amplitude(self, final_ids):
        """Return the amplitude with the given final pids.
           If no suitable amplitude is found, retun None.
        """
        if not isinstance(final_ids, list):
            raise self.PhysicsObjectError,\
                "The final particle ids %s must be a list of integer." \
                %str(final_ids)

        if any([not isinstance(i, int) for i in final_ids]):
            raise self.PhysicsObjectError,\
                "The final particle ids %s must be a list of integer." \
                %str(final_ids)

        # Sort the given id list first
        final_ids.sort()
        # Search in the amplitude list of given final particle number
        if self.get_amplitudes(len(final_ids)):
            for amp in self.get_amplitudes(len(final_ids)):
                ref_list = [l['id'] \
                                for l in amp.get('process').get_final_legs()]
                ref_list.sort()
                if ref_list == final_ids:
                    return amp
                
        return None

    def get_amplitudes(self, partnum):
        """Return the n-body decay amplitudes.
           partnum = n.
           If the amplitudes do not exist, return none.
        """
        #check the validity of arguments
        if not isinstance(partnum, int):
            raise self.PhysicsObjectError, \
                "The particle number %s must be an integer."  %str(partnum)

        try:
            return self.get('decay_amplitudes')[partnum]
        except KeyError:
            logger.info('The amplitudes of %d for final particle number % d do not exist' % (self['pdg_code'],partnum))
            return []

    def set_amplitudes(self, partnum, value):
        """Set the n_body decay_amplitudes.
           partnum: n, 
           value: the decay_amplitudes that is tried to assign.
        """

        #Check the value by check_amplitudes
        if isinstance(value, DecayAmplitudeList):
            if self.check_amplitudes(partnum, value):
                self.decay_amplitudes[partnum] = value
        elif isinstance(value, list) and \
                all([isinstance(a, DecayAmplitude) for a in value]):
            new_value = DecayAmplitudeList(value)
            if self.check_amplitudes(partnum, new_value):
                self.decay_amplitudes[partnum] = new_value
        else:
            raise self.PhysicsObjectError, \
                "The input must be a list of decay amplitudes."
        
              
    def find_channels(self, partnum, model, min_br=0):
        """ Function for finding decay channels up to the final particle
            number given by max_partnum.
            The branching ratios and err are recalculated in the end.
            Algorithm:
            1. Any channel must be a. only one vertex, b. an existing channel
               plus one vertex.
            2. Given the maxima channel order, the program start looking for
               2-body decay channels until the channels with the given order.
            3. For each channel order,
               a. First looking for any single vertex with such order and 
                  construct the channel. Setup the has_idpart property.
               b. Then finding the channel from off-shell sub-level channels.
                  Note that any further decay of on-shell sub-level channel
                  is not a new channel. For each sub-level channel, go through
                  the final legs and connect vertex with the right order
                  (with the aid of connect_channel_vertex function).
               c. If the new channel does not exist in 'decay_channels',
                  then if the new channel has no identical particle, append it.
                  If the new channel has identical particle, check if it is not
                  equivalent with the existing channels. Then append it.
         """

        # Check validity of argument
        if not isinstance(partnum, int):
            raise self.PhysicsObjectError, \
                "Max final particle number %s should be integer." % str(partnum)
        if not isinstance(model, DecayModel):
            raise self.PhysicsObjectError, \
                "The second argument %s should be a DecayModel object." \
                % str(model)            
        if not self in model['particles']:
            raise self.PhysicsObjectError, \
                "The model %s does not contain particle %s" \
                % (model.get('name'), self.get('name'))

        # If vertexlist has not been found before, run model.find_vertexlist
        if not model['vertexlist_found']:
            logger.info("Process the model to add decay informations.")
            model.find_vertexlist()

        # Find stable particles of this model
        model.get('stable_particles')

        # If the channel list exist, return.
        if (partnum, True) in self['decay_channels'].keys() or \
                (partnum, False) in self['decay_channels'].keys():
            logger.info("Particle %s has found channels in %d-body level. " \
                            % (self['name'], partnum) +\
                        "No channel search will not proceed." )
            return

        # Set empty decay_channels for stable particles
        if self.get('is_stable'):
            for num in range(2, partnum+1):
                self['decay_channels'][(num, True)] = ChannelList()
                self['decay_channels'][(num, False)] = ChannelList()
            logger.info("Particle %s is stable. " % self['name'] +
                            "No channel search will not proceed." )
            return

        # Running the coupling constants
        model.running_externals(abs(eval(self.get('mass'))))
        model.running_internals()

        self['apx_decaywidth'] = 0.
        # Find channels from 2-body decay to partnum-body decay channels.
        for clevel in range(2, partnum+1):
            self.find_channels_nextlevel(model, min_br)

        # Update decay width err and branching ratios, but not the total width
        self.update_decay_attributes(False, True, True)

    def find_channels_nextlevel(self, model, min_br=0):
        """ Find channels from all the existing lower level channels to 
            the channels with one more final particle. It is implemented when
            the channels lower than clevel are all found. The user could
            first setup channels by find_channels and then use this function
            if they want to have high level decay.
            NOTE that the width of channels are added sucessively into mother
            particle during the search.
            NO calculation of branching ratios and width err are performed!
            For the search of two-body decay channels, the 2body_massdiff
            is recorded in the end."""

        # Find the next channel level
        try:
            clevel = max([key[0] for key in self['decay_channels'].keys()])+1
        except ValueError:
            clevel = 2

        # Initialize the item in dictionary
        self['decay_channels'][(clevel, True)] = ChannelList()
        self['decay_channels'][(clevel, False)] = ChannelList()

        # Return if this particle is stable
        if self.get('is_stable'):
            logger.info("Particle %s is stable." %self['name'] +\
                            "No channel search will not proceed.")
            return
          
        connect_channel_vertex = self.connect_channel_vertex
        check_repeat = self.check_repeat

        # If there is a vertex in clevel, construct it
        if (clevel, True) in self['decay_vertexlist'].keys() or \
                (clevel, False) in self['decay_vertexlist'].keys():
            for vert in (self.get_vertexlist(clevel, True) + \
                             self.get_vertexlist(clevel, False)):
                
                if len(vert.get('legs')) >3:
                    if not model.keep_Npoint(vert, self):
                        continue
                
                temp_channel = Channel()
                temp_vert = copy.deepcopy(vert)

                # Set the leg number (starting from 2)
                [l.set('number', i) \
                     for i, l in enumerate(temp_vert.get('legs'), start=2)]

                # The final one has 'number' as 1, pid as anti pid of mother
                # (to be consistent with trunk)
                temp_vert.get('legs')[-1]['number'] = 1
                temp_vert.get('legs')[-1]['state'] = False
                temp_vert.get('legs')[-1]['id'] = self.get_anti_pdg_code()
                temp_channel['vertices'].append(temp_vert)

                # Run initial setups for the new channel
                temp_channel.initial_setups(model, True)

                # Setup the 'has_idpart' property
                if Channel.check_idlegs(temp_vert):
                    temp_channel['has_idpart'] = True


                # Add width to total width if onshell
                if temp_channel.get_onshell(model):
                    if temp_channel.has_goldstone(model):
                        continue
                    self['apx_decaywidth'] += temp_channel.get_apx_decaywidth(model)                    

                # Append this channel after all the setups.
                self.get_channels(clevel, temp_channel.get_onshell(model)).\
                    append(temp_channel)
        
        start = time.time()
        repeat_time = 0
        # Go through sub-channels and try to add vertex to reach partnum
        for sub_clevel in range(max((clevel - model.get_max_vertexorder()+1),2),
 clevel):
            # The vertex level that should be combine with sub_clevel
            vlevel = clevel - sub_clevel+1
            # Go through each 'off-shell' channel in the given sub_clevel.
            # Decay of on-shell channel is not a new channel.
            tot = len(self.get_channels(sub_clevel, False))
            for nb,sub_c in enumerate(self.get_channels(sub_clevel, False)):
                if (nb + 1) % 100 == 0:
                    logger.info('%i / %i: %ss' % (nb+1, tot, int(time.time()-start)))
                # Scan each leg to see if there is any appropriate vertex
                for index, leg in enumerate(sub_c.get_final_legs()):

                    # Get the particle even for anti-particle leg.
                    inter_part = model.get_particle(abs(leg['id']))

                    # If this inter_part is stable, do not attach vertices to it
                    if inter_part.get('is_stable'):
                        continue
                    # Get the vertexlist in vlevel
                    # Both on-shell and off-shell vertex 
                    # should be considered.
                    vlist_a = inter_part.get_vertexlist(vlevel, True)
                    vlist_b = inter_part.get_vertexlist(vlevel, False)

                    minv_max = eval(self['mass']) - \
                                    sum([eval(model.get_particle(abs(l['id']))['mass'])
                                         for l in sub_c.get_final_legs() if l!=leg])
                    
                    allow_qcd=True # if two colored particle are lower than the pion. 
                                   # those computation are meaningless.
                    if 0 < minv_max.real < 0.100:
                        logger.warning("WARNING: Mass gap lower than pion mass for decay of %s "
                                                             % self['pdg_code'])
                        logger.warning("decay into colored particle will be remove.")
                        allow_qcd=False
                        if model.get_particle(abs(leg['id']))['color'] != 1:
                            continue
                        
                        
                    # Find appropriate vertex
                    for vert in (vlist_a + vlist_b):
                        # Connect sub_channel to the vertex
                        # the connect_channel_vertex will
                        # inherit the 'has_idpart' from sub_c

                        if not allow_qcd:
                            nb_colored = sum([1 for l in vert['legs'] if \
                                     model.get_particle(abs(l['id']))['color'] != 1])
                            if nb_colored >=1:
                                continue

                        temp_c = self.connect_channel_vertex(sub_c, index, 
                                                             vert, model)
                        temp_c_o = temp_c.get_onshell(model)
                        
                        if temp_c_o and temp_c.has_goldstone(model):
                            continue
                        
                        # Append this channel if it is new
                        rstart = time.time()
                        status = self.check_repeat(clevel, temp_c_o, temp_c)
                        repeat_time += time.time() - rstart
                        if not status:

                            # Check gauge dependence
                            if not temp_c.check_gauge_dependence(model):
                                continue
                            # Calculate the width if onshell
                            # Add to the apx_decaywidth of mother particle
                            if temp_c_o:
                                self['apx_decaywidth'] += temp_c.\
                                    get_apx_decaywidth(model)
                            
                            self.get_channels(clevel, temp_c_o).append(temp_c)
            
            # add valid 3-body decay which start from a valid 2-body decay but 
            #where the sub-decay is forbidden S1 > S2 (S2 > S3 S3)
            #with 2* MS2 < MS1 but 2*MS3> MS2 so S2 is not onshell!
            tot = len(self.get_channels(sub_clevel, True))
            for nb,sub_c in enumerate(self.get_channels(sub_clevel, True)):
                if (nb + 1) % 100 == 0:
                    logger.info('%i / %i: %ss' % (nb+1, tot, int(time.time()-start)))
                # Scan each leg to see if there is any appropriate vertex
                for index, leg in enumerate(sub_c.get_final_legs()):

                    # Get the particle even for anti-particle leg.
                    inter_part = model.get_particle(abs(leg['id']))

                    # If this inter_part is stable, do not attach vertices to it
                    if inter_part.get('is_stable'):
                        continue
                    # Get the vertexlist in vlevel
                    # ONLY off-shell vertex should be considered.
                    vlist_b = inter_part.get_vertexlist(vlevel, False)


                    # Find appropriate vertex
                    for vert in vlist_b:
                        # Connect sub_channel to the vertex
                        # the connect_channel_vertex will
                        # inherit the 'has_idpart' from sub_c

                        temp_c = self.connect_channel_vertex(sub_c, index, 
                                                             vert, model)
                        temp_c_o = temp_c.get_onshell(model)
                        if not temp_c_o:
                            # if offshell, this starts to be extremelly tricky
                            #and time consuming to found all the fake onshell
                            #decay. We therefore neglect all of them. This starts 
                            #to be problematic for 4 body decay (for valid 2-body)
                            #so this should be fine.
                            continue
                        elif temp_c.has_goldstone(model):
                            continue
                        # Append this channel if it is new
                        rstart = time.time()
                        status = self.check_repeat(clevel, temp_c_o, temp_c)
                        repeat_time += time.time() - rstart
                        if not status:
                            # Check gauge dependence
                            if not temp_c.check_gauge_dependence(model):
                                continue
                            # Calculate the width if onshell
                            # Add to the apx_decaywidth of mother particle
                            self['apx_decaywidth'] += temp_c.\
                                    get_apx_decaywidth(model)
                            
                            self.get_channels(clevel, temp_c_o).append(temp_c)               

            
            if hasattr(self, 'check_repeat_tag'):          
                del self.check_repeat_tag

        # For two-body decay, record the maximal mass difference
        if clevel == 2:
            for channel in self.get_channels(2, True):
                mass_diff = abs(eval(self['mass'])) - sum(channel.get('final_mass_list'))
                if mass_diff > self['2body_massdiff']:
                    self.set('2body_massdiff', mass_diff)


        ## Sort the channels by their width
        #self.get_channels(clevel, False).sort(channelcmp_width)


        # Group channels into amplitudes
        self.group_channels_2_amplitudes(clevel, model, min_br)
        

    def connect_channel_vertex(self, sub_channel, index, vertex, model):
        """ Helper function to connect a vertex to one of the legs 
            in the channel. The argument 'index' specified the position of
            the leg in sub_channel final_legs which will be connected with 
            vertex. If leg is for anti-particle, the vertex will be transform
            into anti-part with minus vertex id."""

        # Copy the vertex to prevent the change of leg number
        new_vertex = copy.deepcopy(vertex)
        
        # Setup the final leg number that is used in the channel.
        leg_num = max([l['number'] for l in sub_channel.get_final_legs()])


        # Find the correct vertex id for decays of antiparticles,
        # using dict in model.conj_int_dict
        # particle legs must be in front of anti-particle legs
        if sub_channel.get_final_legs()[index]['id'] < 0:
            new_vertex['id'] = model['conj_int_dict'][vertex['id']]
            for leg in new_vertex['legs']:
                leg['id']  = model.get_particle(leg['id']).get_anti_pdg_code()

            new_vertex['legs'][:-1] = sorted(new_vertex['legs'][:-1],
                                             key=lambda leg: leg['id'],
                                             reverse=True)


        # Legs continue the number of the final one in sub_channel,
        # except the first and the last one ( these two should be connected.)
        for leg in new_vertex['legs'][1:len(new_vertex['legs'])-1]:
            # For the first and final legs, this loop is useless
            # only for reverse leg id in case of mother leg is anti-particle.
            leg_num += 1
            leg['number'] = leg_num


        # Assign correct number to each leg in the vertex.
        # The first leg follows the mother leg.
        new_vertex['legs'][-1]['number'] = \
            sub_channel.get_final_legs()[index]['number']
        new_vertex['legs'][0]['number'] = new_vertex['legs'][-1]['number']
        
        # Combining vertex with channel
        new_channel = Channel()

        # New vertex is first
        new_channel['vertices'].append(new_vertex)

        # Then extend the vertices of the old channel
        # deepcopy is necessary to get new legs        
        new_channel['vertices'].extend(copy.deepcopy(sub_channel['vertices']))



        # Setup properties of new_channel (This slows the time)
        new_channel.initial_setups(model, True)

        # The 'has_idpart' property of descendent of a channel 
        # must derive from the mother channel and the vertex
        # (but 'id_part_list' will change and should not be inherited)
        new_channel['has_idpart'] = (sub_channel['has_idpart'] or \
                                         bool(Channel.check_idlegs(vertex)))

        return new_channel


    def check_repeat(self, clevel, onshell, channel):
        """ Check whether there is any equivalent channels with the given 
            channel. Use the check_channels_equiv function."""

        # old method (much slower), keep-it for comparison
        #previous_answer = any([Channel.check_channels_equiv(other_c, channel)\
        #               for other_c in self.get_channels(clevel, onshell)
        #               if abs(sum(other_c['final_mass_list'])-
        #                 sum(channel['final_mass_list']))< 0.01])

        # To optimize the check, check if the final_mass is the same first
        # this will significantly reduce the number of call to
        # check_channels_equiv
        # further speed improve, store the tag and avoid to recompute the final mass.
        if not hasattr(self, 'check_repeat_tag'):
            self.check_repeat_tag = {(clevel,onshell): {}}
        elif not (clevel,onshell) in  self.check_repeat_tag.keys():
            self.check_repeat_tag[(clevel,onshell)] =  {}
      
        tag = channel.get('tag')
        mass = 100 * sum(channel['final_mass_list']) // 1 #avoid numerical inacuracy
        if mass in self.check_repeat_tag[(clevel,onshell)]:
            repeat = any(tag == other_tag 
                     for other_tag in self.check_repeat_tag[(clevel,onshell)][mass])
        else:
            repeat = False
        
        #add the information to the dictionary to fasten following calls
        if not repeat:
            if mass in self.check_repeat_tag[(clevel,onshell)]:
                self.check_repeat_tag[(clevel,onshell)][mass].append(tag)
            else:
                self.check_repeat_tag[(clevel,onshell)][mass] = [tag]
        return repeat
                  
    def group_channels_2_amplitudes(self, clevel, model, min_br=0):
        """ After the channels is found, combining channels with the same 
            final states into amplitudes.
            NO CALCULATION of branching ratio at this stage!
            clevel: the number of final particles
            model: the model use in find_channels
        """

        if not isinstance(clevel, int):
            raise self.PhysicsObjectError, \
                "The channel level %s must be an integer." % str(clevel)

        if not isinstance(model, DecayModel):
            raise self.PhysicsObjectError, \
                "The model must be an DecayModel object."

        # Reset the value of decay_amplitudes
        self.set_amplitudes(clevel, DecayAmplitudeList())
        # Sort the order of onshell channels according to their final mass list.
        self.get_channels(clevel, True).sort(channelcmp_final)

        total_width = self.get('apx_decaywidth')
        for channel in self.get_channels(clevel, True):
            found = False
            # Record the final particle id.
            final_pid = sorted([l.get('id') for l in channel.get_final_legs()])
        
            # Check if there is a amplitude for it. Since the channels with 
            # the similar final states are put together. Use reversed order
            # of loop.
            for amplt in reversed(self.decay_amplitudes[clevel]):
                # Do not include the first leg (initial id)
                if sorted([l.get('id') for l in amplt['process']['legs'][1:]])\
                        == final_pid:
                    amplt.add_std_diagram(channel)
                    found = True
                    break

            # If no amplitude is satisfied, initiate a new one.
            if not found:
                #print 'start : ', channel.nice_string()
                self.get_amplitudes(clevel).append(DecayAmplitude(channel,
                                                                  model))
        if min_br:
            total_width = self.get('apx_decaywidth')
        # calculate the apx_decaywidth and for every amplitude.
        # apr_br WILL NOT be calculated!
        for amp in list(self.get_amplitudes(clevel)):
            approx_width = amp.get('apx_decaywidth')
            if min_br:
                if approx_width or total_width:
                    br = approx_width / total_width
                else:
                    br = 0
                if br.real < min_br:
                    self.decay_amplitudes[clevel].remove(amp)
                
        self.get_amplitudes(clevel).sort(amplitudecmp_width)




    # This helper function is obselete in current algorithm...
    def generate_configlist(self, channel, partnum, model):
        """ Helper function to generate all the configuration to add
            vertetices to channel to create a channel with final 
            particle number as partnum"""

        current_num = len(channel.get_final_legs())
        configlist = []
        limit_list = [model.get_particle(abs(l.get('id'))).get_max_vertexorder() for l in channel.get_final_legs()]

        for leg_position in range(current_num):
            if limit_list[leg_position] >= 2:
                ini_config = [0]* current_num
                ini_config[leg_position] = 1
                configlist.append(ini_config)

        for i in range(partnum - current_num -1):
            new_configlist = []
            # Add particle one by one to each config
            for config in configlist:
                # Add particle to each leg if it does not exceed limit_list
                for leg_position in range(current_num):
                    # the current config + new particle*1 + mother 
                    # <= max_vertexorder
                    if config[leg_position] + 2 <= limit_list[leg_position]:
                        temp_config = copy.deepcopy(config)
                        temp_config[leg_position] += 1
                        if not temp_config in new_configlist:
                            new_configlist.append(temp_config)
            if not new_configlist:
                break
            else:
                configlist = new_configlist

        # Change the format consistent with max_vertexorder
        configlist = [[l+1 for l in config] for config in configlist]
        return configlist
                                            
                    
#===============================================================================
# DecayParticleList
#===============================================================================
class DecayParticleList(base_objects.ParticleList):
    """A class to store list of DecayParticle, Particle is also a valid
       element, but will automatically convert to DecayParticle"""

    def __init__(self, init_list=None, force=False):
        """Creates a new particle list object. If a list of physics 
        object is given, add them."""

        list.__init__(self)

        if init_list is not None:
            for object in init_list:
                self.append(object, force)

    def append(self, object, force=False):
        """Append DecayParticle, even if object is Particle"""

        if not force:
            assert self.is_valid_element(object), \
                "Object %s is not a valid object for the current list" %repr(object)

        if isinstance(object, DecayParticle):
            list.append(self, object)
        else:
            list.append(self, DecayParticle(object, force))

    def generate_dict(self):
        """Generate a dictionary from particle id to particle.
        Include antiparticles.
        """

        particle_dict = {}

        for particle in self:
            particle_dict[particle.get('pdg_code')] = particle
            if not particle.get('self_antipart'):
                antipart = copy.copy(particle)
                antipart.set('is_part', False)
                particle_dict[antipart.get_pdg_code()] = antipart

        return particle_dict
    
#===============================================================================
# DecayModel: Model object that is used in this module
#===============================================================================
class DecayModel(model_reader.ModelReader):
    """DecayModel object is able construct the decay vertices for
       all its particles by find_vertexlist. When the user try to get stable
       particles, it will find all stable particles automatically according to
       the particle mass and interactions by find_stable_particles functions.
       The run of find_channels uses the function of DecayParticle. 
       Note that Particle objects will be converted into DecayParticle 
       either during the initialization or when the set function is used.
    """
    sorted_keys = ['name', 'particles', 'parameters', 'interactions', 
                   'couplings', 'lorentz', 
                   'stable_particles', 'vertexlist_found',
                   'reduced_interactions', 'decay_groups', 'max_vertexorder',
                   'decaywidth_list', 
                   'helascalls',
                   'ab_model', 'abmodel_generated', 'coupling_dict','parameter_dict'
                  ]

    def __init__(self, init_dict = {}, force=False):
        """Reset the particle_dict so that items in it is 
           of DecayParitcle type"""

        dict.__init__(self)
        self.default_setup()

        assert isinstance(init_dict, dict), \
            "Argument %s is not a dictionary" % repr(init_dict)

        # Must set particles first so that the particle_dict
        # can point to DecayParticle
        # Futhermore, the set of interactions can have correct particle_dict
        if 'particles' in init_dict.keys():
            p =  init_dict['particles'][0]
            self.set('particles', init_dict['particles'], force)

        self['particle_dict'] = {}
        self.get('particle_dict')
        # Do not copy the particle_dict, it may be old version and point only
        # to Particle rather than DecayParticle.
        for item in init_dict.keys():
            if item != 'particles' and item != 'particle_dict':
                try:
                    self.set(item, init_dict[item], force)
                except:
                    pass

        
    def default_setup(self):
        """The particles is changed to ParticleList"""
        super(DecayModel, self).default_setup()
        self['particles'] = DecayParticleList()
        # Other properties
        self['vertexlist_found'] = False
        self['max_vertexorder'] = 0
        self['conj_int_dict'] = {}

        self['helascalls'] = {}
        self['decay_groups'] = []
        self['reduced_interactions'] = []
        self['stable_particles'] = []
        # Width from the value of param_card.
        self['decaywidth_list'] = {}
        # Properties for abstract model
        self['ab_model'] = AbstractModel()
        self['abmodel_generated'] = False
        self['invalid_Npoint'] = []
        

    def get_sorted_keys(self):
        return self.sorted_keys

    def get(self, name):
        """ Evaluate some special properties first if the user request. """

        if name == 'stable_particles' and not self['stable_particles']:
            self.find_stable_particles()
            return self['stable_particles']
        # reduced_interactions might be empty, cannot judge the evaluation is
        # done or not by it.
        elif (name == 'decay_groups' or name == 'reduced_interactions') and \
                not self['decay_groups']:
            self.find_decay_groups_general()
            return self[name]
        elif name == 'max_vertexorder' and self['max_vertexorder'] == 0:
            self.get_max_vertexorder()
            return self['max_vertexorder']
        else:
            # call the mother routine
            return DecayModel.__bases__[0].get(self, name)
        
    def filter(self, name, value):
        if name == 'vertexlist_found':
            if not isinstance(value, bool):
                raise self.PhysicsObjectError, \
                    "Property %s should be bool type." % name
        if name == 'max_vertexorder':
            if not isinstance(value, int):
                raise self.PhysicsObjectError,\
                    "Property %s should be int type." % name
        if name == 'stable_particles' or name == 'decay_groups':
            if not isinstance(value, list):
                raise self.PhysicsObjectError,\
                    "Property %s should be a list contains several particle list." % name
            for plist in value:                
                if not isinstance(plist, list):
                    raise self.PhysicsObjectError,\
                    "Property %s should be a list contains several particle list." % name
                for p in plist:
                    if not isinstance(p, DecayParticle):
                        raise self.PhysicsObjectError,\
                            "Property %s should be a list contains several particle list." % name

        super(DecayModel, self).filter(name, value)
        
        return True
            
        
    def set(self, name, value, force=False):
        """Change the Particle into DecayParticle"""
        #Record the validity of set by mother routine
        return_value = super(DecayModel, self).set(name, value, force)
        #Reset the dictionaries
        
        if return_value:
            if name == 'particles':
                #Reset dictionaries and decay related properties.
                self['particle_dict'] = {}
                self['got_majoranas'] = None
                self['vertexlist_found'] = False
                self['max_vertexorder'] = 0
                self['decay_groups'] = []
                self['reduced_interactions'] = []
                self['stable_particles'] = []
                self['decaywidth_list'] = {}

                #Convert to DecayParticleList
                self['particles'] = DecayParticleList(value, force=True)
                #Generate new dictionaries with items are DecayParticle
                self.get('particle_dict')
                self.get('got_majoranas')
            if name == 'interactions':
                # Reset dictionaries and decay related properties.
                self['interaction_dict'] = {}
                self['ref_dict_to1'] = {}
                self['ref_dict_to0'] = {}
                self['vertexlist_found'] = False
                self['max_vertexorder'] = 0
                self['decay_groups'] = []
                self['reduced_interactions'] = []
                self['stable_particles'] = []
                self['decaywidth_list'] = {}

                # Generate interactions with particles are DecayParticleList
                # Get particle from particle_dict so that the particles
                # of interaction is the alias of particles of DecayModel
                for inter in self['interactions']:
                    inter['particles']=DecayParticleList([self.get_particle(part.get_pdg_code()) for part in inter['particles']])
                # Generate new dictionaries
                self.get('interaction_dict')
                self.get('ref_dict_to1')
                self.get('ref_dict_to0')

        return return_value

    def get_max_vertexorder(self):
        """find the maxima vertex order (i.e. decay particle number)"""
        if not self['vertexlist_found']:
            logger.warning("Use find_vertexlist before get_max_vertexorder!")
            return

        # Do not include key without any vertexlist in it
        self['max_vertexorder'] = max(sum([[k[0] \
                                            for k in \
                                            p.get('decay_vertexlist').keys() \
                                            if p.get('decay_vertexlist')[k]] \
                                            for p in self.get('particles')], [])
                                      )
        return self['max_vertexorder']

    def get_helascalls(self, number):
        """ Get the helascalls with the number of final particles. Return
            empty list if such helascalls does not exist."""

        return self.get('helascalls').get(number, [])

    def add_helascalls(self, number, diagram):
        """ Add the helascalls with the number of final particles. 
            Use deepcopy of the diagram."""
        
        # diagram should have std_diagram already
        assert diagram['std_diagram']

        std_channel = Channel()
        # Replace the vertices with standard expression 
        # written by IdentifyHelasTagHelasTag.
        std_channel['vertices'] = copy.deepcopy(diagram['std_diagram']['vertices'])
        std_channel['helastag'] = diagram['helastag']


        try:
            self['helascalls'][number].append(std_channel)
        except KeyError:
            self['helascalls'][number] = base_objects.DiagramList([std_channel])


    def collect_helascalls(self, part, number):
        """ Collect diagrams from the ones in the given particle 
            with given final state number that are not in 'helascalls'."""
            
        
        if not isinstance(number, int) or not isinstance(part, DecayParticle):
            raise self.PhysicsObjectError, \
                "Wrong argument types."

        amps = part.get_amplitudes(number)
        diagrams = sum([[d for d in amp['diagrams']] for amp in amps],[])
        calls = self.get_helascalls(number)

        for d in diagrams:
            helastag, std_d = d.get_helas_properties(self)

            found = False
            for num, call in enumerate(calls):
                if not call['helastag']:
                    raise self.PhysicsObjectError, \
                        "The helastag of helascalls should exist."
                if call['helastag'] == helastag:
                    found = True
                    d['helas_number'] = num
                    break

            if not found:
                d['helas_number'] = len(calls)      
                self.add_helascalls(number, d)
                # for empty calls, update the calls to direct to the actual
                # list rather than a fake empty list [].
                if not calls:
                    calls = self.get_helascalls(number)


    def generate_abstract_model(self, force=False):
        """ Generate the abstract particles in AbstractModel."""

        logger.info("Generating the abstract model...")

        self.get('stable_particles')
        # If vertexlist has not been found before, run model.find_vertexlist
        if not self['vertexlist_found']:
            logger.info("Vertexlist of this model has not been searched."+ \
                "Automatically run the model.find_vertexlist()")
            self.find_vertexlist()

        # Add all particles, including stable ones 
        # (could appear in final states)
        self['ab_model'].setup_particles(self.get('particles'), force)

        # Add interactions, except ones with all stable particles or ones
        # with radiation process.
        self['ab_model'].setup_interactions(self.get('interactions'), 
                                            self['conj_int_dict'],
                                            force)

    def generate_abstract_amplitudes(self, part, clevel):
        """ Generate the abstract amplitudes in AbstractModel."""

        # Skip search if this particle is stable
        if part.get('is_stable'):
            logger.info("Particle %s is stable." %part['name'] +\
                            "No abstract amplitude will be generated.")
            return
        
        logger.info("Generating the abstract amplitudes of %s %d-body decays..." % (part['name'], clevel))

        self['ab_model'].generate_ab_amplitudes(part.get_amplitudes(clevel))

            
    def find_vertexlist(self, force=False):
        """ Check whether the interaction is able to decay from mother_part.
            Set the 'decay_vertexlist' of the corresponding particles.
            Utilize in finding all the decay table of the whole model.
            force: set vertices even though it's already in the vertex list.
            Note: only construct VertexList with particle as mother,
                  fix the antiparticle vertices by conj_int_dict.
        """

        # Return if self['vertexlist_found'] is True.\
        if self['vertexlist_found']:
            logger.info("Vertexlist has been searched before.")
            return
        # Find the stable particles of this model and do not assign decay vertex
        # to them.
        self.get('stable_particles')
        
        # Valid initial particle list
        ini_list = []

        for part in self.get('particles'):
            if not part.get('is_stable'):
                #All valid initial particles (mass != 0 and is_part == True)
                ini_list.append(part.get_pdg_code())


        #Prepare the vertexlist        
        for inter in self['interactions']:

            if 'type' in inter and inter['type'] != 'base':
                continue #need only the LO interactions.
            #Calculate the particle number
            partnum = len(inter['particles']) - 1

            temp_legs = base_objects.LegList()
            total_mass = 0
            validity = False

            for num, part in enumerate(inter['particles']):

                #Check if the interaction contains valid initial particle
                if part.get_anti_pdg_code() in ini_list:
                    validity = True

                #Create the original legs
                temp_legs.append(base_objects.Leg({'id':part.get_pdg_code()}))
                total_mass += abs(eval(part.get('mass')))
            
            #Exclude interaction without valid initial particle
            if not validity:
                continue
                
            for num, part in enumerate(inter['particles']):

                # Get anti_pdg_code (pid for incoming particle)
                pid = part.get_anti_pdg_code()

                # Exclude illegal initial particle 
                # (including unstable antiparticles)
                if not pid in ini_list:
                    #print 'bypass', pid, 'not in', ini_list 
                    continue

                # Exclude initial particle (or its own antiparticle) 
                # appears in final particles
                # e.g. radiation, t > t h, and p > p~ etc.                
                pid_list = [p.get_pdg_code() for p in inter.get('particles')]

                # N.B. can't use #pid + #antipid here, which fail
                # for neutral particles.
                if pid_list.count(pid)+pid_list.count(-pid) > 1:

                    antipid = part.get_pdg_code()

                    # Simply ignore processes that're not radiation-like.
                    if len(pid_list)>3:
                        #print 'wrong line'
                        continue
                    # If particle is not antiparticle, #pid=1 and #antipid=1
                    if antipid != pid and not (pid_list.count(pid) ==1 and \
                                                   pid_list.count(-pid) ==1):
                        continue
                    # If particle is antiparticle, #pid=2
                    elif antipid == pid and pid_list.count(pid) ==2:
                        continue

                    # For radiation like process, we need to collect them to
                    # check processes which are identical to 
                    # radiation + proper decay


                    # skip the 2nd initial particle (e.g. b b a)
                    if num > [i for i, l in enumerate(pid_list) if l==antipid][0]:
                        continue

                    # Extract the radiation
                    for temp_id in pid_list:
                        if temp_id != pid and temp_id != antipid:
                            self.get_particle(pid)['radiations'].append(temp_id)

                    """# If the "radiation" turn particle into antiparticle
                    # (p > ~p + something)
                    # Include both the "radiation" and "anti-radiation", the
                    # product emitted by anti-particle
                    if pid_list.count(antipid) > 1:
                        antiradiation = [self.get_particle(l).get_anti_pdg_code() for l in radiation]
                        if set(antiradiation) != set(radiation):
                            self.get_particle(pid)['radiative_products'].append(antiradiation)"""
                        
                    continue

                ini_mass = abs(eval(part.get('mass')))
                onshell = ini_mass > (total_mass - ini_mass)

                # Create new legs for the sort later,
                # turn pid of initial particle into anti-one
                temp_legs_new = copy.deepcopy(temp_legs)
                temp_legs_new[num]['id'] = pid

                # Put initial leg in the last 
                # and sort other legs for comparison
                inileg = temp_legs_new.pop(num)

                # IMPORTANT
                # Particle legs must be in front of anti-particle legs
                # for the sake of HELAS output
                temp_legs_new.sort(key=lambda leg: leg['id'],reverse=True)

                temp_legs_new.append(inileg)
                temp_vertex = base_objects.Vertex({'id': inter.get('id'),
                                                   'legs':temp_legs_new})


                # Record the vertex with key = (interaction_id, part_id)
                try:
                    # Raise error if this vertex appears already
                    # (indication of mother appears in final state)
                    if temp_vertex in self.get_particle(pid)['decay_vertexlist'][(partnum, onshell)] and not force:

                        raise self.PhysicsObjectError,\
                            "found duplicated %d-point vertices in %s" \
                            % (partnum+1, self.get_particle(pid)['name'])

                    self.get_particle(pid)['decay_vertexlist'][(\
                            partnum, onshell)].append(temp_vertex)

                # Create VertexList if necessary
                except KeyError:
                    self.get_particle(pid)['decay_vertexlist'][(\
                            partnum, onshell)] = base_objects.VertexList(\
                        [temp_vertex])

        # Set the property vertexlist_found as True and for all particles
        self['vertexlist_found'] = True
        for part in self['particles']:
            part['vertexlist_found'] = True

        # Find complex conjugate interaction pairs
        self.find_conjugate_dict()

        # Check gauge dependence
        self.gauge_dependence_helper()


    def find_conjugate_dict(self):
        """ Find the dictionary connecting an interaction to its complex 
            conjuagate. """

        # Setup the conj_int_dict for interaction id to the id of 
        # its complex conjugate.
        interlist = [(i['id'], [p.get_anti_pdg_code() for p in i['particles']])\
                         for i in self['interactions']]

        for inter in self['interactions']:
            pids = [p.get_pdg_code() for p in inter['particles']]
            conj_pids = [p.get_anti_pdg_code() for p in inter['particles']]

            # Skip searching if the interaction is self-conjugate
            if sorted(pids) == sorted(conj_pids):
                continue

            found = False
            for item in interlist:
                if sorted(item[1]) == sorted(pids) and \
                        inter['id'] != item[0]:

                    self['conj_int_dict'][inter['id']] = item[0]
                    found = True
                    break

            if not found:
                logger.warning('No hermitian conjugate is found for interaction #%d' %inter['id'])


        #fdata = open(os.path.join(MG5DIR, 'models', self['name'], 'vertexlist_dict.dat'), 'w')
        #fdata.write(str(vertexlist_dict))
        #fdata.close()

    @misc.mute_logger(names=['madgraph.diagram_generation'], levels=[40])
    def keep_Npoint(self, vertex, initpart):
        """ """
        get_mass = lambda x: eval(self.get_particle(x['id']).get('mass'))
        init_mass = eval(initpart.get('mass'))

        if vertex['id'] in self['invalid_Npoint']:
            return False

        interaction = self.get('interaction_dict')[vertex['id']]
        decay_parts = [p for p in interaction['particles']]
        
        if len([1 for p in decay_parts if abs(p['pdg_code'])==abs(initpart['pdg_code'])]) >1:
            self['invalid_Npoint'].append(vertex['id'])
            return False

        # create a process:
        process = base_objects.ProcessDefinition()
        process['model'] = self
        process['orders'] = interaction['orders']
        for order in self.get('coupling_orders'):
            if order not in interaction['orders']:
                process['orders'][order] = 0
        leglist = base_objects.MultiLegList()
        for i,particle in enumerate(interaction['particles']):
            newleg = base_objects.MultiLeg({#'number': i+1,
                                       'ids': [particle['pdg_code'] * (1 if particle['is_part'] else -1)],
                                       'state': True
                                       })
            if abs(newleg['ids'][0]) == abs(initpart['pdg_code']):
                leglist.insert(0, newleg)
            else:
                leglist.append(newleg)
            
            
            
        process['legs'] = leglist
        
        myprocdef = base_objects.ProcessDefinitionList()
        myprocdef.append(process) 
        myproc = diagram_generation.MultiProcess(myprocdef, optimize=False)

        one_offshell = False
        one_subdiag = False
        for amp in myproc['amplitudes']:
            for proc in amp['diagrams']:
                # check if they are substructure
                if not one_subdiag and len(proc['vertices']) >1:
                    one_subdiag = True
                else:
                    continue
                
                #check that all substructure are valid
                #remove if any radiation and two times the same particle in a vertex
                for v in proc['vertices']:
                    if any([get_mass(l)==0 for l in v.get('legs')]):
                        self['invalid_Npoint'].append(vertex['id'])
                        return False

                    ids = set(abs(l['id']) for l in v.get('legs'))
                    if len(ids) != len(vertex.get('legs')):
                        self['invalid_Npoint'].append(vertex['id'])
                        return False

                # check onshell/offshell status                
                prev_mass = 0
                for v in proc['vertices'][:-1]:
                    propa =  v.get('legs')[-1]
                    if propa['number'] == 1:
                        other_mass = sum([get_mass(p) for p in v.get('legs')[1:]])
                        if init_mass.real > other_mass.real + prev_mass:
                            # update the mass for the nex propagator
                            prev_mass = sum([get_mass(p) for p in v.get('legs')[1:-1]])
                            # special for 3 body to check for fake onshell
                            if len(vertex['legs']) == 4:
                                propa_mass = get_mass(v.get('legs')[-1])
                                other_v = proc['vertices'][-1]
                                end_mass = sum([get_mass(p) for p in other_v['legs']
                                                if p['number']!=1])
                                if propa_mass.real < end_mass.real:
                                    one_offshell = True
                                    continue
                                
                        else:
                            one_offshell = True
                            continue
                                
                                
                            
                            
                
        
        if not one_subdiag:
            return True
        elif one_offshell:
            return True
        else:
            return False                
                        

                                        
                     
        



    def gauge_dependence_helper(self):
        """ Check the potential gauge dependence of vertices, i.e.
            3-pt interaction + radiation. 
            Run inside find_vertexlist, after radiations are identified and
            complex-conjugate interactions are found.
            
            For a N (>3) interactions the algorithm is the following.
            1) generate a process object (same particles and same order)
            2) reject the diagram if some of the diagrams have more than one 
            interactions.
            
            REPLACE BY check_Npoint with a check of onshell/offshell in addiation
        """
        return

            

    def color_multiplicity_def(self, colorlist):
        """ Defines the two-body color multiplicity. It is applied in the 
            get_color_multiplicity in channel object.
            colorlist is the a list of two color indices.
            This function will return a list of all possible
            color index of the "mother" particle and the corresponding
            color multiplicities. """
        
        # Raise error if the colorlist is not the right format.
        if not isinstance(colorlist, list):
            raise self.PhysicsObjectError,\
                "The argument must be a list."

        if any([not isinstance(i, int) for i in colorlist]):
            raise self.PhysicsObjectError,\
                "The argument must be a list of integer elements."

        # Sort the colorlist and 
        colorlist.sort()
        color_tuple = tuple(colorlist)
        # The dictionary of color multiplicity
        # The key is the final_color structure and the value is
        # [(ini_color_1, multiplicity_1), (ini_color_2, multiplicity_2), ...]
        color_dict = { \
            # The trivial key and value
            # These define the convention we used.
            (1, 1): [(1, 1)],
            (1, 3): [(3, 1)],
            (1, 8): [(8, 1./2)],
            (1, 6): [(6, 1)],
            (3, 3): [(1, 3), (8, 0.5), (3, 1), (6, 1)],
            # (3, 6) ->  (3, 2 rather than 3), 
            (3, 6): [(3, 2), (8, 3./4)],
            (3, 8): [(3, 0.5*(3-1./3)), (6, 1)],
            (6, 6): [(1, 6), (8, 4./3)],
            (6, 8): [(3, 2), (6, 2-1./3)],
            (8, 8): [(1, 8)],
            # 3-body decay color structure
            # Quick reference
            (3, 3, 8): [(1, 8), (8, 3-1./3)],
            (1, 3, 8): [(3, 3)],
            (1, 3, 3): [(8, 1)],
            (3, 8, 8): [(3, 64./9)]
            }

        return color_dict[color_tuple]





    def read_param_card(self, param_card):
        """Read a param_card and set all parameters and couplings as
        members of this module"""

        self.set_parameters_and_couplings(param_card)
        for param, value in self.get('parameter_dict').items():
            exec("globals()[\'%s\'] = %s" % (param, value))
        for param, value in self.get('coupling_dict').items():
            exec("globals()[\'%s\'] = %s" % (param, value))        

        for particle in self.get('particles'):
            pid = abs(particle['pdg_code'])
            value = self.get('parameter_dict')[particle['width']]
            self['decaywidth_list'][(pid, True)] = float(value.real)
            
        global amZ0, aS
        amZ0 = aS



    def running_externals(self, q, loopnum=2):
        """ Recalculate external parameters at the given scale. """
        
        # Raise error for wrong type of q
        if not isinstance(q, int) and not isinstance(q, long) and \
                not isinstance(q, float):
            raise self.PhysicsObjectError, \
                "The argument %s should be numerical type." %str(q)

        # Declare global value. amZ0 is the alpha_s at Z pole
        global aS, amZ0, mdl_amZ0
        # Define all functions used
        for func in self['functions']:
            exec("def %s(%s):\n   return %s" % (func.name,
                                                ",".join(func.arguments),
                                                func.expr))

        # Setup the alpha_s at different scale
        amt = 0.
        amb = 0.
        amc = 0.
        a_out = 0.

        # Setup parameters
        # MZ, MB are already read in from param_card
    
    
        # get Z mass
        Z = self.get_particle(23)
        if not Z:
            MZ_ref = 91.118
        else:
            MZ_ref = self['parameter_dict'][Z['mass']].real
        
        B = self.get_particle(5)
        if not B:
            MB_ref = 4.7
        else:
            MB_ref = self['parameter_dict'][B['mass']].real
            if not MB_ref:
                MB_ref = 4.7
                
        C = self.get_particle(4)
        if not C:
            MC_ref = 1.42
        else:
            MC_ref = self['parameter_dict'][C['mass']].real
            if not MC_ref:
                MC_ref = 1.42

        # Calculate alpha_s at the scale q
        if q < MB_ref:
            # Running the alpha_s from MZ_ref to MB_ref (fermion_num = 5)
            t = 2 * math.log(MB_ref/MZ_ref)
            amb = self.newton1(t, amZ0, loopnum, 5.)
            if q< MC_ref:
                # Running the alpha_s from MB_ref to MC_ref (fermion_num = 4)
                t = 2.*math.log(MC_ref/MB_ref)
                amc = self.newton1(t, amb, loopnum, 4.)
                try:
                    # Running the alpha_s from MC_ref to q
                    t = 2.*math.log(q/MC_ref)
                    a_out = self.newton1(t, amc, loopnum, 3.)
                except ValueError:
                    a_out = amc
            else:
                # Running the alpha_s from MB_ref to q (fermion_num = 4)
                t = 2.*math.log(q/MB_ref)
                a_out = self.newton1(t, amb, loopnum, 4.)
        else:
            # Running the alpha_s from MZ_ref to MB_ref (fermion_num = 5)
            t = 2 * math.log(q/MZ_ref)
            a_out = self.newton1(t, amZ0, loopnum, 5.)

        # Save the result alpha_s
        aS = a_out


    def newton1(self, t, a_in, loopnum, nf):
        """ Calculate the running strong coupling constant from a_in
            using the t as the energy running factor, 
            loop number given by loopnum, number of fermions by nf."""

        # Setup the accuracy.
        tol = 5e-5

        # Functions used in the beta function
        b0 = {}
        c1 = {}
        c2 = {}
        delc = {}

        for i in range(3,6):
            b0[i] = (11. - 2.*i/3.)/4./math.pi
            c1[i] = (102. - 38./3.*i)/4./math.pi/(11. - 2.*i/3.)
            c2[i] = (2857./2. - 5033.*i/18. + 325* i**2 /54.)\
                /16./math.pi**2/(11. - 2.*i/3.)
            delc[i] = math.sqrt(4*c2[i] - c1[i]**2)

        f2 = lambda x: (1./x + c1[nf] * math.log((c1[nf]*x)/(1. + c1[nf]*x)))
        f3 = lambda x: (1./x + 0.5*c1[nf] * \
                         math.log((c2[nf]* x**2)\
                                      /(1. +c1[nf]*x +c2[nf]* x**2))\
                            -(c1[nf]**2 - 2.*c2[nf])/delc[nf]* \
                            math.atan((2.*c2[nf]*x + c1[nf])/delc[nf]))

        # Return the 1-loop alpha_s
        if loopnum == 1:
            return a_in/(1. + a_in * b0[nf] * t)
        # For higher order correction, setup the initial value of a_out
        else:
            a_out = a_in/(1. + a_in * b0[nf] * t + \
                              c1[nf] * a_in * math.log(1. + a_in * b0[nf] * t))
            if a_out <= 0.:
                a_out = 0.3
            
        # Start the iteration            
        delta = tol +1
        while delta > tol:
            if loopnum == 2:
                f = b0[nf]*t + f2(a_in) - f2(a_out)
                fp = 1./(a_out**2 * (1. + c1[nf]*a_out))
            if loopnum == 3:
                f = b0[nf]*t + f3(a_in) - f3(a_out)
                fp = 1./(a_out**2 * (1. + c1[nf]*a_out + c2[nf]* a_out**2))
            if fp == 0:
                return a_in
            a_out = a_out - f/fp
            delta = abs(f/fp/a_out)

        return a_out

    def running_internals(self):
        """ Recalculate parameters and couplings which depend on
            running external parameters to the given energy scale.
            RUN running_externals before run this function."""

        # Define all functions used
        for func in self['functions']:
            exec("def %s(%s):\n   return %s" % (func.name,
                                                ",".join(func.arguments),
                                                func.expr))

        # External parameters that must be recalculate for different energy
        # scale.
        run_ext_params = ['aS']

        # Extract derived parameters
        derived_parameters = []
        # Take only keys that depend on running external parameters
        ordered_keys = [key for key in self['parameters'].keys() if \
                            key != ('external',) and \
                            any([run_ext in key for run_ext in run_ext_params])]
        # Sort keys, keys depend on fewer number parameters should be
        # evaluated first.
        ordered_keys.sort(key=len)
        for key in ordered_keys:
            derived_parameters += self['parameters'][key]

        # Now calculate derived parameters
        # TO BE IMPLEMENTED use running alpha_s for aS-dependent params
        for param in derived_parameters:
            exec("globals()[\'%s\'] = %s" % (param.name, param.expr))
            if not eval(param.name) and eval(param.name) != 0:
                logger.warning("%s has no expression: %s" % (param.name,
                                                             param.expr))
#             try:
#                 logger.info("Recalculated parameter %s = %f" % \
#                             (param.name, eval(param.name)))
#             except TypeError:
#                 logger.info("Recalculated parameter %s = (%f, %f)" % \
#                             (param.name,\
#                              eval(param.name).real, eval(param.name).imag))
        
        # Extract couplings from couplings that depend on fewer 
        # number of external parameters.
        couplings = []
        # Only take keys that contain running external parameters.
        ordered_keys = [key for key in self['couplings'].keys() if \
                            key != ('external',) and \
                            any([run_ext in key for run_ext in run_ext_params])]
        # Sort keys
        ordered_keys.sort(key=len)
        for key in ordered_keys:
            couplings += self['couplings'][key]

        # Now calculate all couplings
        # TO BE IMPLEMENTED use running alpha_s for aS-dependent couplings
        for coup in couplings:
            exec("globals()[\'%s\'] = %s" % (coup.name, coup.expr))
            if not eval(coup.name) and eval(coup.name) != 0:
                logger.warning("%s has no expression: %s" % (coup.name,
                                                             coup.expr))
            """logger.info("Recalculated coupling %s = (%f, %f)" % \
                        (coup.name,\
                         eval(coup.name).real, eval(coup.name).imag))"""


    def find_decay_groups(self):
        """ Find groups of particles which can decay into each other,
        keeping Standard Model particles outside for now. This allows
        to find particles which are absolutely stable based on their
        interactions.

        Algorithm:

        1. Start with any non-SM particle. Look for all
        interactions which has this particle in them.

        2. Any particles with single-particle interactions with this
        particle and with any number of SM particles are in the same
        decay group.

        3. If any of these particles have decay to only SM
        particles, the complete decay group becomes "sm"
        
        5. Iterate through all particles, to cover all particles and
        interactions.
        """

        self.sm_ids = [1,2,3,4,5,6,11,12,13,14,15,16,21,22,23,24]
        self['decay_groups'] = [[]]

        particles = [p for p in self.get('particles') if \
                     p.get('pdg_code') not in self.sm_ids]

        for particle in particles:
            # Check if particles is already in a decay group
            if particle not in sum(self['decay_groups'], []):
                # Insert particle in new decay group
                self['decay_groups'].append([particle])
                self.find_decay_groups_for_particle(particle)

    def find_decay_groups_for_particle(self, particle):
        """Recursive routine to find decay groups starting from a
        given particle.

        Algorithm:

        1. Pick out all interactions with this particle

        2. For any interaction which is not a radiation (i.e., has
        this particle twice): 

        a. If there is a single non-sm particle in
        the decay, add particle to this decay group. Otherwise, add to
        SM decay group or new decay group.

        b. If there are more than 1 non-sm particles: if all particles
        in decay groups, merge decay groups according to different
        cases:
        2 non-sm particles: either both are in this group, which means
        this is SM, or one is in this group, so the other has to be
        SM, or both are in the same decay group, then this group is SM.
        3 non-sm particles: either 1 is in this group, then the other
        two must be in same group or 2 is in this group, then third
        must also be in this group, or 2 is in the same group, then
        third must be in this group (not yet implemented). No other
        cases can be dealt with.
        4 or more: Not implemented (not phenomenologically interesting)."""
        
        # interactions with this particle which are not radiation
        interactions = [i for i in self.get('interactions') if \
                            particle in i.get('particles') and \
                            i.get('particles').count(particle) == 1 and \
                            (particle.get('self_antipart') or
                             not self.get_particle(particle.get_anti_pdg_code())\
                                 in i.get('particles'))]
                             
        while interactions:
            interaction = interactions.pop(0)
            non_sm_particles = [p for p in interaction.get('particles') \
                                if p != particle and \
                                not p.get('pdg_code') in self.sm_ids and \
                                not (p.get('is_part') and p in \
                                     self['decay_groups'][0] or \
                                     not p.get('is_part') and \
                                     self.get_particle(p.get('pdg_code')) in \
                                     self['decay_groups'][0])]
            group_index = [i for (i, g) in enumerate(self['decay_groups']) \
                           if particle in g][0]

            if len(non_sm_particles) == 0:
                # The decay group of this particle is the SM group
                if group_index > 0:
                    group = self['decay_groups'].pop(group_index)
                    self['decay_groups'][0].extend(group)
                    
            elif len(non_sm_particles) == 1:
                # The other particle should be in my decay group
                particle2 = non_sm_particles[0]
                if not particle2.get('is_part'):
                    particle2 = self.get_particle(particle2.get_anti_pdg_code())
                if particle2 in self['decay_groups'][group_index]:
                    # This particle is already in this decay group,
                    # and has been treated.
                    continue
                elif particle2 in sum(self['decay_groups'], []):
                    # This particle is in a different decay group - merge
                    group_index2 = [i for (i, g) in \
                                    enumerate(self['decay_groups']) \
                                    if particle2 in g][0]
                    group = self['decay_groups'].pop(max(group_index,
                                                      group_index2))
                    self['decay_groups'][min(group_index, group_index2)].\
                                                        extend(group)
                else:
                    # Add particle2 to this decay group
                    self['decay_groups'][group_index].append(particle2)

            elif len(non_sm_particles) > 1:
                # Check if any of the particles are not already in any
                # decay group. If there are any, let another particle
                # take care of this interaction instead, later on.

                non_checked_particles = [p for p in non_sm_particles if \
                                         (p.get('is_part') and not p in \
                                          sum(self['decay_groups'], []) or \
                                          not p.get('is_part') and not \
                                          self.get_particle(\
                                                     p.get_anti_pdg_code()) in \
                                          sum(self['decay_groups'], []))
                                         ]

                if not non_checked_particles:
                    # All particles have been checked. Analyze interaction.

                    if len(non_sm_particles) == 2:
                        # Are any of the particles in my decay group already?
                        this_group_particles = [p for p in non_sm_particles \
                                                if p in self['decay_groups'][\
                                                                   group_index]]
                        if len(this_group_particles) == 2:
                            # There can't be any conserved quantum
                            # number! Should be SM group!
                            group = self['decay_groups'].pop(group_index)
                            self['decay_groups'][0].extend(group)
                            continue
                        elif len(this_group_particles) == 1:
                            # One particle is in the same group as this particle
                            # The other (still non_sm yet) must be SM group.
                            particle2 = [p for p in non_sm_particles \
                                             if p != this_group_particles[0]][0]
                            if not particle2.get('is_part'):
                                particle2 = self.get_particle(particle2.get_anti_pdg_code())

                            group_index2 = [i for (i, g) in \
                                                enumerate(self['decay_groups'])\
                                                if particle2 in g][0]
                            group_2 = self['decay_groups'].pop(group_index2)
                            self['decay_groups'][0].extend(group_2)

                        else:
                            # If the two particles are in another same group,
                            # this particle must be the SM particle.
                            # Transform the 1st non_sm_particle into particle
                            particle1 = non_sm_particles[0]
                            if not particle1.get('is_part'):
                                particle1 = self.get_particle(\
                                    particle1.get_anti_pdg_code())
                            # Find the group of particle1
                            group_index1 = [i for (i, g) in \
                                            enumerate(self['decay_groups']) \
                                            if particle1 in g][0]

                            # If the other non_sm_particle is in the same group
                            # as particle1, try to merge this particle to SM
                            if non_sm_particles[1] in \
                                    self['decay_groups'][group_index1]:
                                if group_index > 0:
                                    group = self['decay_groups'].pop(group_index)
                                    self['decay_groups'][0].extend(group)

                    if len(non_sm_particles) == 3:
                        # Are any of the particles in my decay group already?
                        this_group_particles = [p for p in non_sm_particles \
                                                if p in self['decay_groups'][\
                                                                   group_index]]
                        if len(this_group_particles) == 2:
                            # Also the 3rd particle has to be in this group.
                            # Merge.
                            particle2 = [p for p in non_sm_particles if p not \
                                         in this_group_particles][0]
                            if not particle2.get('is_part'):
                                particle2 = self.get_particle(\
                                                  particle2.get_anti_pdg_code())
                            group_index2 = [i for (i, g) in \
                                            enumerate(self['decay_groups']) \
                                            if particle2 in g][0]
                            group = self['decay_groups'].pop(max(group_index,
                                                              group_index2))
                            self['decay_groups'][min(group_index, group_index2)].\
                                                                extend(group)
                        if len(this_group_particles) == 1:
                            # The other two particles have to be in
                            # the same group
                            other_group_particles = [p for p in \
                                                     non_sm_particles if p not \
                                                     in this_group_particles]
                            particle1 = other_group_particles[0]
                            if not particle1.get('is_part'):
                                particle1 = self.get_particle(\
                                                  particle1.get_anti_pdg_code())
                            group_index1 = [i for (i, g) in \
                                            enumerate(self['decay_groups']) \
                                            if particle1 in g][0]
                            particle2 = other_group_particles[0]
                            if not particle2.get('is_part'):
                                particle2 = self.get_particle(\
                                                  particle2.get_anti_pdg_code())
                            group_index2 = [i for (i, g) in \
                                            enumerate(self['decay_groups']) \
                                            if particle2 in g][0]

                            if group_index1 != group_index2:
                                # Merge groups
                                group = self['decay_groups'].pop(max(group_index1,
                                                                  group_index2))
                                self['decay_groups'][min(group_index1,
                                                      group_index2)].\
                                                                   extend(group)

                        # One more case possible to say something
                        # about: When two of the three particles are
                        # in the same group, the third particle has to
                        # be in the present particle's group. I'm not
                        # doing this case now though.

                    # For cases with number of non-sm particles > 3,
                    # There are also possibilities to say something in
                    # particular situations. Don't implement this now
                    # however.

    def find_decay_groups_general(self):
        """Iteratively find decay groups, valid for vertices in all orders
           SM particle is defined as all MASSLESS particles.
           Algrorithm:
           1. Establish the reduced_interactions
              a. Read non-sm particles only
                 (not in sm_ids and not in decay_groups[0])
              b. If the particle appears in this interaction before, 
                 not only stop read it but also remove the existing one.
              c. If the interaction has only one particle,
                 move this particle to SM-like group and void this interaction.
              d. If the interaction has no particle in it, delete it.   
           2. Iteratively reduce the interaction
              a. If there are two particles in this interaction,
                 they must be in the same group. 
                 And we can delete this interaction since we cannot draw more
                 conclusion from it.
              b. If there are only one particle in this interaction,
                 this particle must be SM-like group
                 And we can delete this interaction since we cannot draw more
                 conclusion from it.
              c. If any two particles in this interaction already belong to the 
                 same group, remove the two particles. Delete particles that 
                 become SM-like as well. If this interaction becomes empty 
                 after these deletions, delete this interaction.
              d. If the iteration does not change the reduced_interaction at all
                 stop the iteration. All the remaining reduced_interaction must
                 contain at least three non SM-like particles. And each of 
                 them belongs to different groups.
           3. If there is any particle that has not been classified,
              this particle is lonely i.e. it does not related to 
              other particles. Add this particle to decay_groups.
        """
        
        # Setup the SM particles and initial decay_groups, reduced_interactions
        self['decay_groups'] = [[]]
        # self['reduced_interactions'] contains keys in 'id' as interaction id,
        # 'particle' as the reduced particle content, and 'groupid_list' as
        # the decay group id list
        self['reduced_interactions'] = []
        sm_ids = []

        # Setup the original 'SM' particles, i.e. particle without mass.
        sm_ids = [p.get('pdg_code') for p in self.get('particles')\
                      if abs(eval(p.get('mass'))) == 0.]
        self['decay_groups'] = [[p for p in self.get('particles')\
                                     if abs(eval(p.get('mass'))) == 0.]]

        #Read the interaction information and setup
        for inter in self.get('interactions'):
            temp_int = {'id':inter.get('id'), 'particles':[]}
            for part in inter['particles']:
                #If this particle is anti-particle, convert it.
                if not part.get('is_part'):
                    part = self.get_particle(part.get_anti_pdg_code())

                #Read this particle if it is not in SM
                if not part.get('pdg_code') in sm_ids and \
                   not part in self['decay_groups'][0]:
                    #If pid is not in the interaction yet, append it
                    if not part in temp_int['particles']:
                        temp_int['particles'].append(part)
                    #If pid is there already, remove it since double particles
                    #is equivalent to none.
                    else:
                        temp_int['particles'].remove(part)

            # If there is only one particle in this interaction, this must in SM
            if len(temp_int['particles']) == 1:
                # Remove this particle and add to decay_groups
                part = temp_int['particles'].pop(0)
                self['decay_groups'][0].append(part)

            # Finally, append only interaction with nonzero particles
            # to reduced_interactions.
            if len(temp_int['particles']):
                self['reduced_interactions'].append(temp_int)
            # So interactions in reduced_interactions are all 
            # with non-zero particles in this stage

        # Now start the iterative interaction reduction
        change = True
        while change:
            change = False
            for inter in self['reduced_interactions']:
                #If only two particles in inter, they are in the same group
                if len(inter['particles']) == 2:
                    #If they are in different groups, merge them.
                    #Interaction is useless.

                    # Case for the particle is in decay_groups
                    if inter['particles'][0] in sum(self['decay_groups'], []):
                        group_index_0 =[i for (i,g) in\
                                        enumerate(self['decay_groups'])\
                                        if inter['particles'][0] in g][0]

                        # If the second one is also in decay_groups, merge them.
                        if inter['particles'][1] in sum(self['decay_groups'], []):
                            if not inter['particles'][1] in \
                                    self['decay_groups'][group_index_0]:
                                group_index_1 =[i for (i,g) in \
                                                enumerate(self['decay_groups'])\
                                                if inter['particles'][1] 
                                                in g][0]
                                # Remove the outer group
                                group_1 = self['decay_groups'].pop(max(\
                                          group_index_0, group_index_1))
                                # Merge with the inner one
                                self['decay_groups'][min(group_index_0, \
                                                 group_index_1)].extend(group_1)
                        # The other one is no in decay_groups yet
                        # Add inter['particles'][1] to the group of 
                        # inter['particles'][0]
                        else:
                            self['decay_groups'][group_index_0].append(
                                inter['particles'][1])
                    # Case for inter['particles'][0] is not in decay_groups yet.
                    else:
                        # If only inter[1] is in decay_groups instead, 
                        # add inter['particles'][0] to its group.
                        if inter['particles'][1] in sum(self['decay_groups'], []):
                            group_index_1 =[i for (i,g) in \
                                            enumerate(self['decay_groups'])\
                                            if inter['particles'][1] in g][0]
                            # Add inter['particles'][0]
                            self['decay_groups'][group_index_1].append(
                                inter['particles'][0])

                        # Both are not in decay_groups
                        # Add both particles to decay_groups
                        else:
                            self['decay_groups'].append(inter['particles'])

                    # No matter merging or not the interaction is useless now. 
                    # Kill it.
                    self['reduced_interactions'].remove(inter)
                    change = True

                # If only one particle in this interaction,
                # this particle must be SM-like group.
                elif len(inter['particles']) == 1:
                    if inter['particles'][0] in sum(self['decay_groups'], []):
                        group_index_1 =[i for (i,g) in \
                                        enumerate(self['decay_groups'])\
                                        if inter['particles'][0] in g][0]
                        # If it is not, merge it with SM.
                        if group_index_1 > 0:
                            self['decay_groups'][0].\
                                extend(self['decay_groups'].pop(group_index_1))

                    # Inter['Particles'][0] not in decay_groups yet, 
                    # add it to SM-like group
                    else:
                        self['decay_groups'][0].extend(inter['particles'])

                    # The interaction is useless now. Kill it.
                    self['reduced_interactions'].remove(inter)
                    change = True
                
                # Case for more than two particles in this interaction.
                # Remove particles with the same group.
                elif len(inter['particles']) > 2:
                    #List to store the id of each particle's decay group
                    group_ids = []
                    # This list is to prevent removing elements during the 
                    # for loop to create errors.
                    # If the value is normal int, the particle in this position 
                    # is valid. Else, it is already removed. 
                    ref_list = range(len(inter['particles']))
                    for i, part in enumerate(inter['particles']):
                        try:
                            group_ids.append([n for (n,g) in \
                                              enumerate(self['decay_groups']) \
                                              if part in g][0])
                        # No group_ids if this particle is not in decay_groups
                        except IndexError:
                            group_ids.append(None)
                            continue
                        
                        # If a particle is SM-like, remove it!
                        # (necessary if some particles turn to SM-like during
                        # the loop then we could reduce the number and decide
                        # groups of the rest particle
                        if group_ids[i] == 0:
                            ref_list[i] = None
                            change = True

                        # See if any valid previous particle has the same group.
                        # If so, both the current one and the previous one
                        # is void
                        for j in range(i):
                            if (group_ids[i] == group_ids[j] and \
                                group_ids[i] != None) and ref_list[j] != None:
                                # Both of the particles is useless for 
                                # the determination of parity
                                ref_list[i] = None
                                ref_list[j] = None
                                change = True
                                break
                    
                    # Remove the particles label with None in ref_list
                    # Remove from the end to prevent errors in list index.
                    for i in range(len(inter['particles'])-1, -1, -1):
                        if ref_list[i] == None:
                            inter['particles'].pop(i)

                    # Remove the interaction if there is no particle in it
                    if not len(inter['particles']):
                        self['reduced_interactions'].remove(inter)

                # Start a new iteration...


        # Check if there is any particle that cannot be classified.
        # Such particle is in the group of its own.
        for part in self.get('particles'):
            if not part in sum(self['decay_groups'], []) and \
                    not part.get('pdg_code') in sm_ids:
                self['decay_groups'].append([part])

        # For conveniences, record the decay_group id in particles of 
        # reduced interactions
        for inter in self['reduced_interactions']:
            inter['groupid_list'] = [[i \
                                          for i, g in \
                                          enumerate(self['decay_groups']) \
                                          if p in g][0] \
                                         for p in inter['particles']]

        return self['decay_groups']

    def find_stable_particles(self):
        """ Find stable particles that are protected by parity conservation
            (massless particle is not included). 
            Algorithm:
            1. Find the lightest massive particle in each group.
            2. For each reduced interaction, the group of the particle which
               lightest mass is greater than the sum of all other particles 
               is not (may not be) stable.
            3. Replace the lightest mass of unstable group as its decay products
            4. Repeat 2., until no replacement can be made.   

        """

        # If self['decay_groups'] is None, find_decay_groups first.
        if not self['decay_groups']:
            self.find_decay_groups_general()

        # The list for the lightest particle for each groups
        # The first element is preserved for the SM-like group.
        stable_candidates = [[]]
        # The list for the lightest mass for each groups
        lightestmass_list = [0.]

        # Set the massless particles into stable_particles
        self['stable_particles'] = [[]]
        for p in self.get('particles'):
            if abs(eval(p.get('mass'))) == 0. :
                p.set('is_stable', True)
                self['stable_particles'][-1].append(p)

        # Find lightest particle in each group.
        # SM-like group is excluded.
        for group in self['decay_groups'][1:]:
            # The stable particles of each group is described by a sublist
            # (take degeneracy into account). Group index is the index in the
            # stable_candidates. Suppose the lightest particle is the 1st one.
            stable_candidates.append([group[0]])
            for part in group[1:]:
                # If the mass is smaller, replace the the list.
                if abs(eval(part.get('mass'))) < \
                        abs(eval(stable_candidates[-1][0].get('mass'))) :
                    stable_candidates[-1] = [part]
                # If degenerate, append current particle to the list.
                elif abs(eval(part.get('mass'))) == \
                        abs(eval(stable_candidates[-1][0].get('mass'))) :
                    stable_candidates[-1].append(part)
            # Record the lightest mass into lightestmass_list
            lightestmass_list.append(abs(eval(stable_candidates[-1][0].get('mass'))))


        # Deal with the reduced interaction
        change = True
        while change:
            change = False
            for inter in self['reduced_interactions']:
                # Find the minial mass for each particle
                masslist = [lightestmass_list[inter['groupid_list'][i]] \
                                for i in range(len(inter['particles']))]
                
                # Replace the minial mass to possible decay products
                for i, m in enumerate(masslist):
                    if 2*m > sum(masslist):
                        # Clear the stable_candidates in this group
                        stable_candidates[inter['groupid_list'][i]] = []
                        # The lightest mass becomes the mass of decay products.
                        lightestmass_list[inter['groupid_list'][i]] = \
                            sum(masslist)-m
                        change = True
                        break

        # Append the resulting stable particles
        for stable_particlelist in stable_candidates:
            if stable_particlelist:
                self['stable_particles'].append(stable_particlelist)

        # Set the is_stable label for particles in the stable_particles
        for p in sum(self['stable_particles'], []):
            p.set('is_stable', True)
            self.get_particle(p.get_anti_pdg_code()).set('is_stable', True)

        # Run the advance find_stable_particles to ensure that
        # all stable particles are found
        self.find_stable_particles_advance()

        return self['stable_particles']

    def find_stable_particles_advance(self):
        """ Find all stable particles. 
            Algorithm:
            1. For each interaction, if one particle has mass larger than 
               the other, than this particle's mass is replaced by 
               the sum of its decay products' masses. 
               The 'is_stable' label of this particle is False.
               
            2. Repeat 1., until no change was made after the whole check.
            3. Particles that have never been labeled as unstable are now
               stable particles.
        """

        # Record the mass of all particles
        # Record whether particles can decay through stable_list
        mass={}
        stable_list = {}
        for part in self.get('particles'):
            mass[part.get('pdg_code')] = abs(eval(part.get('mass')))
            stable_list[part.get('pdg_code')] = True

        # Record minimal mass for avioding round-off error.
        m_min = min(min([m for i,m in mass.items() if m > 0.]), 10**(-6))

        # Start the iteration
        change = True
        while change:
            change = False
            for inter in self['interactions']:
                total_m = sum([mass[p.get('pdg_code')] \
                                   for p in inter['particles']])

                # Skip interaction with total_m = 0
                if total_m == 0.:
                    continue

                # Find possible decay for each particle
                for part in inter['particles']:
                    # If not stable particle yet.
                    if not part.get('is_stable'):
                        # This condition is to prevent round-off error.
                        if (2*mass[part.get('pdg_code')]-total_m) > \
                                10**3*mass[part.get('pdg_code')]*\
                                sys.float_info.epsilon:
                            mass[part.get('pdg_code')] = \
                                total_m - mass[part.get('pdg_code')]
                            part['is_stable'] = False
                            stable_list[part.get('pdg_code')] = False
                            change = True
                            break

        # Record the stable particle
        for part in self.get('particles'):
            if stable_list[part.get('pdg_code')]:               
                part.set('is_stable', True)
                self.get_particle(part.get_anti_pdg_code()).set('is_stable', 
                                                                True)
                if not part in sum(self['stable_particles'], []):
                    self['stable_particles'].append([part])


    def find_channels(self, part, max_partnum, min_br=0):
        """ Function that find channels for a particle.
            Call the function in DecayParticle."""
        part.find_channels(max_partnum, self, min_br)

    def find_all_channels(self, max_partnum, collect_helascalls = True,
                          generate_abstract=False, min_br=0):
        """ Function that find channels for all particles in this model.
            Call the function in DecayParticle.
            It also write a file to compare the decay width from 
            param_card and from the estimation of this module."""

        # If vertexlist has not been found before, run model.find_vertexlist
        if not self['vertexlist_found']:
            logger.info("Vertexlist of this model has not been searched."+ \
                "Automatically run the model.find_vertexlist()")
            self.find_vertexlist()

        # If vertexlist has not been found before, run model.find_vertexlist
        if generate_abstract and not self['abmodel_generated']:
            logger.info("AbstractModel for this model has not been generated."\
                            +"Automatically run the model.generate_abstract_model()")
            self.generate_abstract_model()

        # Find stable particles of this model
        self.get('stable_particles')
        logger.info("Found %s stable particles" % len(self['stable_particles']))

        # Run the width of all particles from 2-body decay so that the 3-body
        # decay could use the width from 2-body decay.
        for part in self.get('particles'):
            # Skip search if this particle is stable
            if part.get('is_stable'):
                continue

            # Recalculating parameters and coupling constants 
            self.running_externals(abs(eval(part.get('mass'))))
            self.running_internals()
            #logger.info("Find 2-body channels of %s" %part.get('name'))
            part.find_channels_nextlevel(self, min_br)

            # Update the total width for the later error estimation
            part.update_decay_attributes(True, False, False)

            if generate_abstract:
                self.generate_abstract_amplitudes(part, 2)
            if collect_helascalls:
                self.collect_helascalls(part, 2)
 
        for part in self.get('particles'):
            # Skip search if this particle is stable
            if part.get('is_stable'):
                continue
                
            # Recalculating parameters and coupling constants 
            self.running_externals(abs(eval(part.get('mass'))))
            self.running_internals()

            if max_partnum > 2:
                # After recalculating the parameters, find the channels to the
                # requested level.
                for clevel in range(3, max_partnum+1):
                    logger.info("Find %d-body channels of %s" %(clevel,
                                                                part.get('name')))
                    part.find_channels_nextlevel(self)

                    if generate_abstract:
                        self.generate_abstract_amplitudes(part, clevel)
                    
                    if collect_helascalls:
                        self.collect_helascalls(part, clevel)

            # update the decay attributes for both max_partnum >2 or == 2.
            # The update should include branching ratios and apx_decaywidth_err
            # So the apx_decaywidth_err(s) are correct even for max_partnum ==2.
            part.update_decay_attributes(True, True, True, self)


    def find_all_channels_smart(self, precision, 
                                collect_helascalls = True,
                                generate_abstract=False):
        """ Function that find channels for all particles in this model.
            Decay channels more than three final particles are searched
            when the precision is not satisfied.
            If generate_abstract = True, the abstract amplitudes will be
            generated."""

        # Raise error if precision is not a float
        if not isinstance(precision, float):
            raise self.PhysicsObjectError, \
                "The precision %s should be float type." % str(precision)

        # If vertexlist has not been found before, run model.find_vertexlist
        if not self['vertexlist_found']:
            logger.info("Vertexlist of this model has not been searched."+ \
                "Automatically run the model.find_vertexlist()")
            self.find_vertexlist()

        # If vertexlist has not been found before, run model.find_vertexlist
        if generate_abstract and not self['abmodel_generated']:
            logger.info("AbstractModel for this model has not been generated."\
                            +"Automatically run the model.generate_abstract_model()")
            self.generate_abstract_model()

        # Find stable particles of this model
        self.get('stable_particles')

        # Run the width of all particles from 2-body decay so that the 3-body
        # decay could use the width from 2-body decay.
        for part in self.get('particles'):
            # Skip search if this particle is stable
            if part.get('is_stable'):
                logger.info("Particle %s is stable." %part['name'] +\
                                "No channel search will not proceed.")
                continue

            # Recalculating parameters and coupling constants 
            self.running_externals(abs(eval(part.get('mass'))))
            self.running_internals()

            logger.info("Find 2-body channels of %s" %part.get('name'))
            part.find_channels_nextlevel(self)
            if generate_abstract:
                self.generate_abstract_amplitudes(part, 2)

            if collect_helascalls:
                self.collect_helascalls(part, 2)

        # Search for higher final particle states, if the precision
        # is not satisfied.
        for part in self.get('particles'):
            # Skip search if this particle is stable
            if part.get('is_stable'):
                continue

            # Update the decaywidth_err
            part.update_decay_attributes(False,True,False, self)

            # If the error (ratio to apx_decaywidth) is larger then precision,
            # find next level channels.
            # Running coupling constants first.
            if part.get('apx_decaywidth_err') > precision:
                self.running_externals(abs(eval(part.get('mass'))))
                self.running_internals()

            clevel = 3
            while part.get('apx_decaywidth_err') > precision:
                logger.info("Find %d-body channels of %s" \
                                %(clevel,
                                  part.get('name')))
                part.find_channels_nextlevel(self)
                if generate_abstract:
                    self.generate_abstract_amplitudes(part, clevel)

                if collect_helascalls:
                    self.collect_helascalls(part, clevel)

                # Note that the width is updated automatically in the
                # find_nextlevel
                part.update_decay_attributes(False,True,False, self)
                clevel += 1

            # Finally, update the branching ratios
            part.update_decay_attributes(False, False, True)         



    def write_summary_decay_table(self, name=''):
        """ Write a table to list the total width of all the particles
            and compare to the value in param_card."""
    
        # Write the result to decaywidth_MODELNAME.dat in 'mg5decay' directory
        path = os.path.join(MG5DIR, 'mg5decay')
        if not name:
            fdata = open(os.path.join(path, 
                                      (self['name']+'_decay_summary.dat')),
                         'w')
            logger.info("\nWrite decay width summary to %s \n" \
                            % str(os.path.join(path,
                                               (self['name']+'_decay_summary.dat'))))

        elif isinstance(name, str):
            fdata = open(os.path.join(path, name),'w')
            logger.info("\nWrite decay width summary to %s \n" \
                            % str(os.path.join(path, name)))

        else:
            raise PhysicsObjectError,\
                "The file name of the decay table must be str." % str(name)

        summary_chart = ''
        summary_chart = (str('# DECAY WIDTH COMPARISON \n') +\
                            str('# model: %s \n' %self['name']) +\
                            str('#'*80 + '\n')+\
                            str('#Particle ID    card value     apprx. value  ratio') +\
                            str('   level    err \n')
                        )

        for part in self.get('particles'):
            # For non-stable particles
            if not part.get('is_stable'):
                # For width available in the param_card.
                try:
                    summary_chart +=(str('#%11d    %.4e     %.4e    %4.2f  %3d        %.2e\n'\
                                            %(part.get('pdg_code'), 
                                              self['decaywidth_list']\
                                                  [(part.get('pdg_code'), True)],
                                              part['apx_decaywidth'],
                                              part['apx_decaywidth']/self['decaywidth_list'][(part.get('pdg_code'), True)],
                                              part.get_max_level(),
                                              part['apx_decaywidth_err']
                                              )))
                # For width not available, do not calculate the ratio.
                except KeyError:
                    summary_chart += (str('#%11d    %.4e     %.4e    %s\n'\
                                             %(part.get('pdg_code'), 
                                               0.,
                                               part['apx_decaywidth'],
                                               'N/A')))
                # For width in param_card is zero.
                except ZeroDivisionError:
                    summary_chart += (str('#%11d    %.4e     %.4e    %s\n'\
                                             %(part.get('pdg_code'), 
                                               0.,
                                               part['apx_decaywidth'],
                                               'N/A')))
            # For stable particles
            else:
                try:
                    if abs(self['decaywidth_list'][(part.get('pdg_code'), True)]) == 0.:
                        ratio = 1
                    else:
                        ratio = 0
                    summary_chart += (str('#%11d    %.4e     %s    %4.2f\n'\
                                             %(part.get('pdg_code'), 
                                               self['decaywidth_list']\
                                                   [(part.get('pdg_code'), True)],
                                               'stable    ',
                                               ratio)))
                except KeyError:
                    summary_chart += (str('#%11d    %.4e     %s    %s\n'\
                                             %(part.get('pdg_code'), 
                                               0.,
                                               'stable    ',
                                               '1'
                                               )))
        # Write the summary_chart into file
        fdata.write(summary_chart)
        fdata.close()


    def write_decay_table(self, mother_card_path, format='normal',name = ''):
        """ Functions that write the decay table of all the particles 
            in this model that including the channel information and 
            branch ratio (call the estimate_width_error automatically 
            in the execution) in a file.
            format:
                normal: write only amplitudes
                cmp: add ratio of decay_width to the value in MG4 param_card
                full: also write the channels in each amplitude."""

        # The list of current allowing formats
        allow_formats = ['normal','full','cmp']

        # Raise error if format is wrong
        assert format in allow_formats, \
                "The format must be \'normal\' or \'full\' or \'cmp\'." \
                % str(name)

        # Write the result to decaywidth_MODELNAME.dat in 'mg5decay' directory
        path = os.path.join(MG5DIR, 'mg5decay')

        if not name:
            if format == 'full':
                fdata = open(os.path.join(path,
                                          (self['name']+'_decaytable_full.dat')),
                             'w')
                logger.info("\nWrite full decay table to %s\n"\
                                %str(os.path.join(path,
                                                  (self['name']+'_decaytable_full.dat'))))
            else:
                fdata = open(os.path.join(path,
                                          (self['name']+'_decaytable.dat')),
                             'w')
                logger.info("\nWrite %s decay table to %s\n"\
                                %(format, 
                                  str(os.path.join(path,
                                          (self['name']+'_decaytable.dat')))))

        elif isinstance(name, str):
            fdata = open(os.path.join(path, name),'w')
            logger.info("\nWrite %s decay table to %s\n"\
                            %(format, 
                              str(os.path.join(path,
                                               name))))

        else:
            raise PhysicsObjectError,\
                "The file name of the decay table must be str." % str(name)

        # Write the param_card used first
        fdata0 = open(mother_card_path, 'r')
        fdata.write(fdata0.read())
        fdata0.close()

        # Write header of the table
        spart = ''
        nonspart = ''
        summary_chart = ''
        seperator = str('#'*80 + '\n')
        fdata.write('\n' + seperator + '#\n'*2 +\
                        str('##    EST. DECAY TABLE    ## \n') +\
                        '#\n'*2 + seperator)

        # Header of summary data
        summary_chart = (str('# DECAY WIDTH COMPARISON \n') +\
                            str('# model: %s \n' %self['name']) +\
                            str('#'*80 + '\n')+\
                            str('#Particle ID    card value     apprx. value  ratio') +\
                            str('   level    err \n')
                        )
        # Header of stable particle output
        spart = ('\n' + seperator + \
                     '# Stable Particles \n'+ \
                     seperator+ \
                     '#%8s    Predicted \n' %'ID')

        # Output order is based on pdg code
        pids = sorted([p.get_pdg_code() for p in self['particles']])

        for pid in pids:
            p = self.get_particle(pid)

            # Write the table only for particles with finite width.
            if p.get('apx_decaywidth'):
                nonspart += p.decaytable_string(format)
                # Try to calculate the ratio in summary_chart
                try:
                    summary_chart +=(str('#%11d    %.4e     %.4e    %4.2f  %3d        %.2e\n'\
                                             %(p.get('pdg_code'), 
                                               self['decaywidth_list']\
                                                  [(p.get('pdg_code'), True)],
                                               p['apx_decaywidth'],
                                               p['apx_decaywidth']/self['decaywidth_list'][(p.get('pdg_code'), True)],
                                               p.get_max_level(),
                                               p['apx_decaywidth_err']
                                              )))
                # For width not available, do not calculate the ratio.
                except KeyError:
                    summary_chart += (str('#%11d    %.4e     %.4e    %s\n'\
                                              %(p.get('pdg_code'), 
                                                0.,
                                                p['apx_decaywidth'],
                                                'N/A')))
                # For width in param_card is zero.
                except ZeroDivisionError:
                    summary_chart += (str('#%11d    %.4e     %.4e    %s\n'\
                                              %(p.get('pdg_code'), 
                                                0.,
                                                p['apx_decaywidth'],
                                                'N/A')))

            else:
                # If width = 0.,
                # see if the stable property is predicted.
                if p.get('is_stable'):
                    spart += str('#%8d    %9s \n' % (p.get('pdg_code'), 'Yes'))
                else:
                    spart += str('#%8d    %9s \n' % (p.get('pdg_code'), 'No'))

                # Try to calculate the ratio if there is reference width
                try:
                    if abs(self['decaywidth_list'][(p.get('pdg_code'), True)]) == 0.:
                        ratio = 1
                    else:
                        ratio = 0
                    summary_chart += (str('#%11d    %.4e     %s    %4.2f\n'\
                                              %(p.get('pdg_code'), 
                                                self['decaywidth_list']\
                                                    [(p.get('pdg_code'), True)],
                                                'stable    ',
                                                ratio)))

                # If no width available, write the ratio as 1
                except KeyError:
                    summary_chart += (str('#%11d    %.4e     %s    %s\n'\
                                             %(p.get('pdg_code'), 
                                               0.,
                                               'stable    ',
                                               '1'
                                               )))
                    
        # Print summary_chart, stable particles, and finally unstable particles
        fdata.write(summary_chart)
        fdata.write(spart)
        fdata.write(nonspart)
        fdata.close()


    def find_nextlevel_ratio(self):
        """ Find the ratio of matrix element square for channels decay to
            next level."""

        pass


    # Helper function to review helas calls
    def write_helas_collection(self, name = ''):
        """ Functions that write the helascalls of the model, attached with
            a list of actual diagrams that belong to the helas call. """


        # Write the result to decaywidth_MODELNAME.dat in 'mg5decay' directory
        path = os.path.join(MG5DIR, 'mg5decay')

        if not name:
            fdata = open(os.path.join(path,
                                      (self['name']+'_helascollection.dat')),
                         'w')
            logger.info("\nWrite %s Helas collection to %s\n"\
                            %(format, 
                              str(os.path.join(path,
                                               (self['name']+'_helascollection.dat')))))

        elif isinstance(name, str):
            fdata = open(os.path.join(path, name),'w')
            logger.info("\nWrite %s Helas collection to %s\n"\
                            %(format, 
                              str(os.path.join(path,
                                               name))))

        else:
            raise PhysicsObjectError,\
                "The file name of the decay table must be str." % str(name)


        collection = ''

        # Output order is based on pdg code
        pids = sorted([p.get_pdg_code() for p in self['particles']])
        
        for key in sorted(self['helascalls'].keys()):
            collection += '%d-body decay: %d calls\n' %(key, len(self.get_helascalls(key)))
            for num, call in enumerate(self.get_helascalls(key)):
                collection += call.nice_string(decay_info=False)+'\n'
                for pid in pids:
                    p = self.get_particle(pid)
                    for amp in p.get_amplitudes(key):
                        for c in amp['diagrams']:
                            if c['helas_number'] == num:
                                collection += "   "+c.nice_string(decay_info=False)+'\n'
            

                    
        # Print collection
        fdata.write(collection)
        fdata.close()


    # Helper Function for reading MG4 param_card
    # And compare with our apx_decaywidth
    def read_MG4_param_card_decay(self, param_card):
        """Read the decay width in MG4 param_card and 
           compare the width with our estimation."""

        if not os.path.isfile(param_card):
            raise MadGraph5Error, \
                "No such file %s" % param_card
    
        # Read in param_card
        logger.info("\nRead MG4 param_card: %s \n" % str(param_card))
        param_lines = open(param_card, 'r').read().split('\n')

        # Define regular expressions
        re_decay = re.compile(\
            "^decay\s+(?P<pid>\d+)\s+(?P<value>-*\d+\.\d+e(\+|-)\d+)\s*")
        re_two_body_decay = re.compile(\
            "^\s+(?P<br>-*\d+\.\d+e(\+|-)\d+)\s+(?P<nda>\d+)\s+(?P<pid1>-*\d+)\s+(?P<pid2>-*\d+)")
        re_three_body_decay = re.compile(\
            "^\s+(?P<br>-*\d+\.\d+e(\+|-)\d+)\s+(?P<nda>\d+)\s+(?P<pid1>-*\d+)\s+(?P<pid2>-*\d+)\s+(?P<pid3>-*\d+)")

        # Define the decay pid, total width
        pid = 0
        total_width = 0

        # Go through lines in param_card
        for line in param_lines:
            if not line.strip() or line[0] == '#':
                continue
            line = line.lower()
            # Look for decay blocks
            decay_match = re_decay.match(line)
            if decay_match:
                pid = int(decay_match.group('pid'))
                total_width = float(decay_match.group('value'))
                self['decaywidth_list'][(pid, True)] = total_width
                continue
            # If no decay pid available, skip this line.
            if not pid:
                continue

            two_body_match = re_two_body_decay.match(line)
            three_body_match = re_three_body_decay.match(line)

            # Check three_body first!
            # Otherwise it will always to be two body.
            if three_body_match:
                # record the pids and br
                pid1 = int(three_body_match.group('pid1'))
                pid2 = int(three_body_match.group('pid2'))
                pid3 = int(three_body_match.group('pid3'))

                br = float(three_body_match.group('br'))
                final_ids = [pid1, pid2, pid3]
                amp = self.get_particle(pid).get_amplitude(final_ids)
                # If amplitude is found, record the ratio
                if amp:
                    amp['exa_decaywidth'] = br*total_width
                # If not found, show this info
                else:
                    logger.info('No amplitude for %d -> %d %d %d is found.' %\
                                       (pid, pid1, pid2, pid3))

                # Jump to next line. Do not match the two_body_decay
                continue

            # If not three-body, check two-body
            if two_body_match:
                # record the pids and br
                pid1 = int(two_body_match.group('pid1'))
                pid2 = int(two_body_match.group('pid2'))
                br = float(two_body_match.group('br'))
                final_ids = [pid1, pid2]
                amp = self.get_particle(pid).get_amplitude(final_ids)
                # If amplitude is found, record the ratio
                if amp:
                    amp['exa_decaywidth'] = br*total_width
                # If not found, show this info
                else:
                    logger.info('No amplitude for %d -> %d %d is found.' %\
                                       (pid, pid1, pid2))

                # Jump to next line. Do not match the three_body_decay
                continue


#===============================================================================
# Channel: A specialized Diagram object for decay
#===============================================================================

"""parameters for estimating the phase space area"""
c_psarea = 0.8

class Channel(base_objects.Diagram):
    """Channel: a diagram that describes a certain decay channel
                with on shell condition, apprximated matrix element, 
                phase space area, and decay width.
                There are several helper static methods.
                The check_idlegs will return the identical legs of the
                given vertex. The check_channels_equiv will check the
                equivalence of two channels.                
    """

    sorted_keys = ['vertices',
                   'orders',
                   'onshell', 'ini_pid', 'final_legs',
                   'has_idpart', 'id_part_list',
                   'tag', 'helastag', 'helas_number', 'std_diagram',
                   'apx_matrixelement_sq', 'apx_psarea', 'apx_decaywidth',
                   'apx_decaywidth_nextlevel', 'apx_width_calculated']#,
                   #'potential_gauge_dependence']

    def default_setup(self):
        """Default values for all properties"""
        self['vertices'] = base_objects.VertexList()
        self['orders'] = {}
        
        # New properties
        self['onshell'] = 0
        self['ini_pid'] = 0
        self['final_legs'] = base_objects.LegList()

        # This property denotes whether the channel has 
        # identical particles in it.
        self['has_idpart'] = False

        # (real) DiagramTag
        self['tag'] = []
        # IdentifyHelasTag
        self['helastag'] = []
        # the number of the corresponding helas calls
        self['helas_number'] = None
        # diagram written by IdentifyHelasTag
        self['std_diagram'] = None

        # old properties for check_channels_equiv, removable
        # The position of the identicle particles with pid as keys.
        self['id_part_list'] = {}
        self['final_mass_list'] = 0

        # Decay width related properties.
        self['apx_matrixelement_sq'] = 0.
        self['s_factor'] = 1
        self['apx_psarea'] = 0.
        self['apx_decaywidth'] = 0.
        # branch ratio is multiply by 100.
        self['apx_decaywidth_nextlevel'] = 0.
        self['apx_width_calculated'] = False

        # Potential gauge dependence
        #self['potential_gauge_dependence'] = False

        # Properties for abstract Channel
        # [pseudo_ab_inter_ids, pseudo_ab_intermediate_pids, pseudo_ab_finalids]
        # generated by add_ab_diagram in AbstractModel,
        # finalids is in the order of final_legs
        self['abstract_type'] = [[], [], []]
        self['fermionfactor'] = 1
        

    def filter(self, name, value):
        """Filter for valid diagram property values."""
        
        if name in ['apx_matrixelement_sq', 'apx_psarea', 
                    'apx_decaywidth', 'apx_br',
                    'apx_decaywidth_nextlevel']:
            if not isinstance(value, float):
                raise self.PhysicsObjectError, \
                    "Value %s is not a float" % str(value)
        
        if name == 'onshell' or name == 'has_idpart' or \
                name == 'apx_width_calculated':
            if not isinstance(value, bool):
                raise self.PhysicsObjectError, \
                        "%s is not a valid onshell condition." % str(value)

        return super(Channel, self).filter(name, value)
    
    def get(self, name, model=None):
        """ Check the onshell condition before the user get it. 
            And recalculate the apx_decaywidth_nextlevel if the 
            model is provided.
        """

        if name == 'onshell':
            logger.info("It is suggested to get onshell property from get_onshell function")

        if name == 'apx_decaywidth_nextlevel' and model:
            return self.get_apx_decaywidth_nextlevel(model)

        if name == 'tag' and not self['tag']:
            self['tag'] = diagram_generation.DiagramTag(self)

        if name == 'helastag' and not self['helastag'] and model:
            self['helastag'] = IdentifyHelasTag(self, model)
            

        return super(Channel, self).get(name)

    def get_sorted_keys(self):
        """Return particle property names as a nicely sorted list."""
        return self.sorted_keys

    def calculate_orders(self, model):
        """Calculate the actual coupling orders of this channel,
           negative vertex id is interepret as positive one 
           (the CPT counterpart)."""

        coupling_orders = {}
        for vertex in self['vertices']:
            if vertex.get('id') == 0: continue
            vid = vertex.get('id')
            couplings = model.get('interaction_dict')[abs(vertex.get('id'))].\
                        get('orders')
            for coupling in couplings.keys():
                try:
                    coupling_orders[coupling] += couplings[coupling]
                except:
                    coupling_orders[coupling] = couplings[coupling]

        self.set('orders', coupling_orders)

    def nice_string(self, decay_info=True):
        """ Add width/width_nextlevel to the nice_string"""
        mystr = super(Channel, self).nice_string()
        if self['vertices'] and decay_info:
            if self['onshell']:
                mystr +=" (width = %.3e)" % self['apx_decaywidth']
            else:
                mystr +=" (est. further width = %.3e)" % self['apx_decaywidth_nextlevel']              

            #if self['potential_gauge_dependence']:
            #    mystr +="*PGD"

        return mystr

    def initial_setups(self, model, force):
        """ Setup useful properties for a new channel.
            N.B. 'has_idpart' is not implemented because there is a faster way
            to get it from connect_channel_vertex. """

        self.get_initial_id(model)
        self.get_final_legs(force)
        self.get_onshell(model)
        self.calculate_orders(model)

    def get_anti_initial_id(self):
        """ Return the anti id of initial particle"""

        return self.get('vertices')[-1].get('legs')[-1].get('id')

    def get_initial_id(self, model=None):
        """ Return the id of initial particle"""

        if self['ini_pid'] != 0:
            return self['ini_pid']
        elif model:
            self['ini_pid'] = model.get_particle(self.get_anti_initial_id()).get_anti_pdg_code()
            return self['ini_pid']
        else:
            raise self.PhysicsObjectError, "No model is provided to get initial id."

    def check_gauge_dependence(self, model):
        """ compare the vertex ids to see if anyone of them may be illed."""

        parts = [model.get_particle(l['id']) for l in self.get_final_legs()]
        initial = model.get_particle(self.get_anti_initial_id())
        for i,part in enumerate(parts):
            if part.get('spin') % 2 == 0:
                continue
            if model['parameter_dict'][part.get('mass')] == 0:
                base = [l['id'] for l in self.get_final_legs()]
                base.pop(i)
                base.sort()
                
                for decay in initial.get_channels(len(base), True):
                    tmp = [l['id'] for l in decay.get_final_legs()]
                    tmp.sort()
                    if base == tmp:
                        return False
        return True
            


    def get_final_legs(self, force=False):
        """ Return a list of the final state legs."""

        # Reset leglist of force = True
        if force:
            self['final_legs'] = base_objects.LegList()

        if not self['final_legs']:
            for vert in self.get('vertices'):
                for leg in vert.get('legs')[:-1]:
                    if not leg.get('number') in [l.get('number') \
                                                 for l in self['final_legs']]\
                                                 and leg.get('number') > 1:
                        self['final_legs'].append(leg)

        return self['final_legs']
        
    def get_onshell(self, model):
        """ Evaluate the onshell condition with the aid of get_final_legs"""
        if not isinstance(self['onshell'], bool):
            # Check if model is valid
            if not isinstance(model, base_objects.Model):
                raise self.PhysicsObjectError, \
                    "The argument %s must be a model." % str(model)

            self['final_mass_list'] =sorted([abs(eval(model.get_particle(l.get('id')).get('mass'))) \
                                                for l in self.get_final_legs()])
            # Use the mass of antiparticle for mother mass
            ini_mass = abs(eval(model.get_particle(self.get_anti_initial_id()).get('mass')))
            # ini_mass = ini_mass.real
            self['onshell'] = ini_mass > sum(self['final_mass_list'])

        return self['onshell']
    
    def has_goldstone(self, model):

            # Check if model is valid
            assert isinstance(model, base_objects.Model), "The argument %s must be a model." % str(model)
                    
            return any(model.get_particle(l.get('id'))['type'] =='goldstone' \
                       for l in self.get_final_legs())


    def get_helas_properties(self, model):
        """ return helastag and std_diagram, construct them if necessary. """

        if not self['helastag'] and model:
            self['helastag'] = IdentifyHelasTag(self, model)

        if not self['std_diagram'] and model:
            self['std_diagram'] = self['helastag'].diagram_from_tag(model)

        return self['helastag'], self['std_diagram']

    def get_fermion_factor(self, model):
        """ Get the fermion_factor, same as get_fermion_factor 
            in HelasAmplitude."""

        # Record the fermion_order using numbers as keys
        order_dict = {}

        # Setup the order_dict of external legs
        for l in self.get_final_legs():
            if model.get_particle(l.get('id')).is_fermion:
                order_dict[l.get('id')] = [l.get('number'), []]
            else:
                order_dict[l.get('id')] = []

        # The initial leg is transformed as outgoing legs
        ini_part = model.get_particle(self.get_initial_id(model))
        if ini_part.is_fermion():
            order_dict[ini_part.get_pdg_code()] = [1, []]
        else:
            order_dict[ini_part.get_pdg_code()] = []

        # Find fermion_order in vertices except the identical one and 
        # initial one
        for i, vert in enumerate(self['vertices'][:-1]):

            # If the vert ended up with fermion, find its fermion_mother
            fermion_mother = None
            if ini_part.is_fermion():
                if vert['id'] > 0:
                    inter = model.get_interaction(vert['id'])
                else:
                    inter = model.get_interaction(model['conj_int_dict'][vert['id']])
                pdg_codes = [p.get_pdg_code() for p in inter['particles']]
                ini_index = pdg_codes.index(ini_part.get_anti_pdg_code())
                # part_index is the index of partner in vertex['legs']
                if ini_index % 2 ==0:
                    part_index = ini_index
                    mother_code = pdg_codes[ini_index+1]
                else:
                    part_index = ini_index - 1
                    mother_code = pdg_codes[ini_index-1]

                fermion_mother = vert['legs'][part_index]

                # Check if this vertex has the right order
                if fermion_mother['id'] != mother_code:
                    raise PhysicsObjectError, \
                        "The order of particle in interaction is not standard."

            # Refind the fermion_order except for initial vertex
            if i != len(self['vertices'])-2:
                order_dict[vert['legs'][-1]['id']] = \
                    self.get_fermion_order(vert, order_dict, fermion_mother)

        # Find fermion_factor in initial vertex
        vert = self['vertices'][-2]
        if vert['id'] > 0:
            inter = model.get_interaction(vert['id'])
        else:
            inter = model.get_interaction(model['conj_int_dict'][vert['id']])
        pdg_codes = [p.get_pdg_code() for p in inter['particles']]
        
        # Construct helper_vertex with 
        # legs = [ final+ initial legs in initial vertex, auxiliary_initial_leg]
        help_legs = base_objects.LegList(\
            [base_objects.Leg({'id': code}) for code in pdg_codes])
        help_legs.append(base_objects.Leg({'id': 0}))
        helper_vertex = base_objects.Vertex({'legs': help_legs})

        # Using this helper_vertex to get final fermion_order
        final_order = self.get_fermion_order(helper_vertex, order_dict, None)

        self['fermionfactor'] = self.sign_flips_to_order(final_order)

    @staticmethod    
    def get_fermion_order(vert, order_dict, fermion_mother):
        """ Get the fermion_order. similar to the get_fermion_order in
            HelasWavefunction. """

        new_order = []

        # Include the order of mother first
        if fermion_mother:
            new_order.extend(order_dict[fermion_mother['id']][1])

        # fermion_order of the fermionic legs
        fermions = [order_dict[l['id']] \
                        for l in vert['legs'][:-1] \
                        if len(order_dict[l['id']]) == 2 and \
                        l != fermion_mother]
        # fermion_order of the bosonic legs
        bosons = [order_dict[l['id']] \
                      for l in vert['legs'][:-1] \
                      if len(order_dict[l['id']]) == 1]

        # include fermionic legs first
        for i in range(0, len(fermions), 2):
            new_order.append(fermions[i][0])
            new_order.append(fermions[i+1][0])
            new_order.extend(fermions[i][1])
            new_order.extend(fermions[i+1][1])
        # Then for bosons
        for b in bosons:
            new_order.extend(b)

        # Add mother_number for fermion
        if fermion_mother:
            return [order_dict[fermion_mother['id']][0], new_order]
        else:
            return new_order

    @staticmethod
    def sign_flips_to_order(self, fermions):
        """Gives the sign corresponding to the number of flips needed
        to place the fermion numbers in order.
        N.B. copy from helas_object.py """

        # Perform bubble sort on the fermions, and keep track of
        # the number of flips that are needed

        nflips = 0

        for i in range(len(fermions) - 1):
            for j in range(i + 1, len(fermions)):
                if fermions[j] < fermions[i]:
                    tmp = fermions[i]
                    fermions[i] = fermions[j]
                    fermions[j] = tmp
                    nflips = nflips + 1

        return (-1) ** nflips

    @staticmethod
    def check_idlegs(vert):
        """ Helper function to check if the vertex has several identical legs.
            If id_legs exist, return a dict in the following format,
            {particle_id: [leg_index1, index2, ...]}
            Otherwise return False.
        """
        lindex_dict = {}
        id_part_list = {}
        # Record the occurence of each leg.
        for lindex, leg in enumerate(vert.get('legs')):
            try:
                lindex_dict[leg['id']].append(lindex)
            except KeyError:
                lindex_dict[leg['id']] = [lindex]

        for key, indexlist in lindex_dict.items():
            # If more than one index for a key, 
            # there are identical particles.
            if len(indexlist) > 1:
                # Record the index of vertex, vertex id (interaction id),
                # leg id, and the list of leg index.
                id_part_list[key] = indexlist

        return id_part_list

    @staticmethod
    def check_channels_equiv(channel_a, channel_b):
        """ Helper function to check if any channel is indeed identical to
            the given channel. (This may happens when identical particle in
            channel.) Use DiagramTag for full comparison.""" 

        return channel_a.get('tag') == channel_b.get('tag')

    # OBSELETE
    def get_idpartlist(self):
        """ Get the position of identical particles in this channel.
            The format of id_part_list is a dictionary with the vertex
            which has identical particles, value is the particle id and
            leg index of identicle particles. Eg.
            id_part_list = {(vertex_index_1, vertex id, pid_1): 
                           [index_1, index_2, ..],
                           (vertex_index_1, vertex id, pid_2): 
                           [index_1, index_2, ..],
                           (vertex_index_2,...): ...}
        """

        if not self['id_part_list']:
            # Check each vertex by check_idlegs
            for vindex, vert in enumerate(self.get('vertices')):
                # Use the id_part_list given by check_idlegs
                id_part_list = Channel.check_idlegs(vert)
                if id_part_list:
                    for key, idpartlist in id_part_list.items():
                        # Record the id_part_list if exists.
                        self['id_part_list'][(vindex, vert.get('id'), 
                                             key)] = id_part_list[key]
                        self['has_idpart'] = True

        return self['id_part_list']

#    # Obselete, replaced by DiagramTag
#    @staticmethod
#    def check_channels_equiv_rec(channel_a, vindex_a, channel_b, vindex_b):
#        """ The recursive function to check the equivalence of channels 
#            starting from the given vertex point.
#            Algorithm:
#            1. Check if the two vertices are the same (in id).
#            2. Compare each the non-identical legs. Either they are both 
#               final legs or their decay chain are the same. The comparision
#               of decay chain is via recursive call of check_channels_equiv_rec
#            3. Check all the identical particle legs, try to match each leg of b
#               for the legs of a of each kind of identical particles. 
#               If a leg of b is fit for one leg of a, do not match this leg of a
#               to other legs of b
#               If any one leg of b cannot be matched, return False.
#            4. If the two channels are the same for all the non-identical legs
#               and are the same for every kind of identical particle,
#               the two channels are the same from the given vertices."""
#        
#        # If vindex_a or vindex_b not in the normal range of index
#        # convert it. (e.g. vindex_b = -1)
#        vindex_a = vindex_a % len(channel_a.get('vertices'))        
#        vindex_b = vindex_b % len(channel_b.get('vertices'))
#
#        # First compare the id of the two vertices.
#        # Do not compare them directly because the number property of legs
#        # may be different.
#        # If vertex id is the same, then the legs id are all the same!
#        if channel_a.get('vertices')[vindex_a]['id'] != \
#                channel_b.get('vertices')[vindex_b]['id']:
#            return False
#
#        # Obsolete when idetical vertex is removed.
#        """# If the vertex is the initial vertex, start from the next vertex
#        if vindex_a == len(channel_a.get('vertices'))-1 and \
#                vindex_b == len(channel_b.get('vertices'))-1 :
#            return Channel.check_channels_equiv_rec(channel_a, -2,
#                                                    channel_b, -2)"""
#        
#        # Find the list of identical particles
#        id_part_list_a=Channel.check_idlegs(channel_a.get('vertices')[vindex_a])
#
#        result = True
#        # For each leg, find their decay chain and compare them.
#        for i, leg_a in \
#                enumerate(channel_a.get('vertices')[vindex_a].get('legs')[:-1]):
#            # The two channels are equivalent as long as the decay chain 
#            # of the two legs must be the same if they are not part of
#            # the identicle particles.
#            if not leg_a.get('id') in id_part_list_a.keys():
#                # The corresponding leg in channel_b
#                leg_b = channel_b.get('vertices')[vindex_b].get('legs')[i]
#
#                # If the 'is final' is inconsistent between a and b,
#                # return False. 
#                # If both are 'is final', end the comparision of these two legs.
#                leg_a_isfinal =  leg_a in channel_a.get_final_legs()
#                leg_b_isfinal =  leg_b in channel_b.get_final_legs()
#                if leg_a_isfinal or leg_b_isfinal:
#                    if leg_a_isfinal and leg_b_isfinal:
#                        continue
#                    else:
#                        # Return false if one is final leg 
#                        # while the other is not.
#                        return False
#                # The case with both legs are not final needs
#                # further analysis
#
#                # Find the next vertex index of the decay chain of 
#                # leg_a and leg_b.
#                for j in range(vindex_a-1, -1, -1):
#                    v = channel_a.get('vertices')[j]
#                    if leg_a in v.get('legs'):
#                        new_vid_a = j
#                        break
#                    
#                for j in range(vindex_b-1, -1, -1):
#                    v = channel_b.get('vertices')[j]
#                    if leg_b in v.get('legs'):
#                        new_vid_b = j
#                        break
#                
#                # Compare the decay chains of the two legs.
#                # If they are already different, return False
#                if not Channel.check_channels_equiv_rec(channel_a, new_vid_a,
#                                                        channel_b, new_vid_b):
#                    return False
#
#        # If the check can proceed out of the loop of legs,
#        # the decay chain is all the same for non-identicle particles.
#        # Return True if there is no identical particles.
#        if not id_part_list_a:
#            return True
#
#
#        # Check each kind of identicle particles
#        for pid, indices_a in id_part_list_a.items():
#            indices_b = copy.copy(indices_a)
#            # Match each leg of channel_b to every leg of channel_a
#            for index_b in indices_b:
#                # Suppose the fit fail
#                this_leg_fit = False
#                # setup leg_b                
#                leg_b = channel_b.get('vertices')[vindex_b].get('legs')[index_b]
#                # Search for match leg in legs from indices_a
#                for i, index_a in enumerate(indices_a):
#                    # setup leg_a
#                    leg_a = channel_a.get('vertices')[vindex_a].get('legs')[index_a]
#                    # Similar to non-identicle particles, but
#                    # we could not return False when one is final leg while
#                    # the other is not since this could due to the wrong
#                    # config used now.
#                    # If the leg is fit (both are final) stop the match of this
#                    # leg.
#                    leg_a_isfinal =  leg_a in channel_a.get_final_legs()
#                    leg_b_isfinal =  leg_b in channel_b.get_final_legs()
#                    if leg_a_isfinal or leg_b_isfinal:
#                        if leg_a_isfinal and leg_b_isfinal:
#                            this_leg_fit = True
#                            indices_a.pop(i)
#                            break
#                        else:
#                            continue
#
#                    # Get the vertex indices for the decay chain of the two
#                    # legs.
#                    for j in range(vindex_a-1, -1, -1):
#                        v = channel_a.get('vertices')[j]
#                        if leg_a in v.get('legs'):
#                            new_vid_a = j
#                            break
#                    for j in range(vindex_b-1, -1, -1):
#                        v = channel_b.get('vertices')[j]
#                        if leg_b in v.get('legs'):
#                            new_vid_b = j
#                            break
#
#                    # If any one of the pairs (leg_a, leg_b) is matched,
#                    # stop the match of this leg
#                    if Channel.check_channels_equiv_rec(channel_a,new_vid_a,
#                                                        channel_b,new_vid_b):
#                        this_leg_fit = True
#                        indices_a.pop(i)
#                        break
#
#                # If this_leg_fit is True, continue to match the next leg of
#                # channel_b. If this_leg_fit remains False, the match of this 
#                # leg cannot match to any leg of channel_a, return False
#                if not this_leg_fit:
#                    return False
#
#        # If no difference is found (i.e. return False),
#        # the two decay chain are the same eventually.
#        return True
#
#    # Helper function (obselete)    
#    @staticmethod
#    def generate_configs(id_part_list):
#        """ Generate all possible configuration for the identical particles in
#            the two channels. E.g. for legs of id=21, index= [1, 3, 5],
#            This function generate a dictionary
#            {leg id ( =21): [[1,3,5], [1,5,3], [3,1,5], [3,5,1], 
#                             [5,1,3], [5,3,1]]}
#            which gives all the possible between the id_legs 
#            in the two channels.
#        """
#        id_part_configs = {}
#        for leg_id, id_parts in id_part_list.items():
#            id_part_configs[leg_id] = [[]]
#
#            # For each index_a, pair it with an index_b.
#            for position, index_a in enumerate(id_parts):
#                # Initiate the new configs of next stage
#                id_part_configs_new = []
#
#                # For each configuration, try to find the index_b
#                # to pair with index_a in the new position.
#                for config in id_part_configs[leg_id]:
#                    # Try to pair index_a with index_b that has not been used
#                    # yet.
#                    for index_b in id_parts:
#                        if not index_b in config:
#                            config_new = copy.copy(config)
#                            config_new.append(index_b)
#                            id_part_configs_new.append(config_new)
#
#                # Finally, replace the configs by the new one.
#                id_part_configs[leg_id] = id_part_configs_new
#
#        return id_part_configs



    def simplify_lorentz(self, match):
        #
        all = match.group('var')
        var = all.split('(')[0]
        if var in ['P','PSlash']:
            nb = int(all.split(',')[1][:-1].strip())
            return '%%(q%s)s' % nb
        else:
            return '1'
        
    @classmethod
    def init_regular_expression(cls):
        dico = dict((`i`, '%s' % ','.join(["\s*-?'?[\w\s]*'?\s*"]*i)) for i in range(1,6))
        cls.lor_pattern = re.compile("""(?<![a-zA-Z])(?P<var>PSlash\(%(3)s\)|
                                        Gamma\(%(3)s\)|
                                        Sigma\(%(4)s\)|
                                        Gamma5\(%(2)s\)|
                                        C\(%(2)s\)|
                                        Epsilon\(%(4)s\)|
                                        Metric\(%(2)s\)|
                                        Identity\(%(2)s\)|
                                        ProjM\(%(2)s\)|
                                        ProjP\(%(2)s\)|
                                        P\(%(2)s\)
                                )
        """ % dico, re.VERBOSE)


    def get_apx_matrixelement_sq(self, model):
        """ Calculate the apx_matrixelement_sq, the estimation for each leg
            is in get_apx_fnrule.
            The color_multiplicity is first searching in the 
            color_multiplicity_def in model object. If no available result,
            use the get_color_multiplicity function.
            For off shell decay, this function will estimate the value
            as if it is on shell."""

        # To ensure the final_mass_list is setup, run get_onshell first
        self.get_onshell(model)
        # Setup the value of matrix element square and the average energy
        # q_dict is to record the flow of energy
        apx_m = 1
        ini_part = model.get_particle(self.get_initial_id())
        # Use the mass of antiparticle
        avg_q = (abs(eval(ini_part.get('mass'))) - sum(self['final_mass_list']))/len(self.get_final_legs())
        q_dict = {}

        # Estimate the width of normal onshell decay.
        if self.get_onshell(model):
            # Go through each vertex and assign factors to apx_m
            # Do not run the identical vertex
            for i, vert in enumerate(self['vertices']):
                
                # Total energy of this vertex
                q_total = 0
                # Color multiplcity
                final_color = []

                q_dict_lor = {}
                # Assign value to apx_m except the mother leg (last leg)
                for lnb,leg in enumerate(vert['legs'][:-1]):
                    # Case for final legs
                    if leg in self.get_final_legs():
                        mass  = abs(eval(model.get_particle(leg.get('id')).\
                                             get('mass')))
                        q_total += (mass + avg_q)
                        q_dict_lor['q%s'%(lnb+1)] = mass + avg_q
                        apx_m *= self.get_apx_fnrule(leg.get('id'),
                                                     avg_q+mass, True, model)
                    # If this is only internal leg, calculate the energy
                    # it accumulated. 
                    # (The value of this leg is assigned before.)
                    else:
                        q_total += q_dict[(leg.get('id'), leg.get('number'))]
                        q_dict_lor['q%s'%(lnb+1)]= q_dict[(leg.get('id'), leg.get('number'))]
                    # Record the color content
                    final_color.append(model.get_particle(leg.get('id')).\
                                           get('color'))
                else:
                    q_dict_lor['q%s'%(lnb+2)] = q_total

                # The energy for mother leg is sum of the energy of its product.
                # Set the q_dict
                q_dict[(vert.get('legs')[-1].get('id'), 
                        vert.get('legs')[-1].get('number'))] = q_total
                # Assign the value if the leg is not initial leg. (propagator),
                # for initial particle, q_total should be M.
                if i < len(self.get('vertices'))-1: 
                    apx_m *= self.get_apx_fnrule(vert.get('legs')[-1].get('id'),
                                                 q_total, False, model)

                else:
                    apx_m *=self.get_apx_fnrule(vert.get('legs')[-1].get('id'),
                                                abs(eval(ini_part.get('mass'))),
                                                True, model)

                # Evaluate the coupling strength
                vertex =  model.get('interaction_dict')[abs(vert.get('id'))]
                lorentz_factor = 0
                for key, v in vertex['couplings'].items():
                    if not hasattr(model, 'lorentz_dict'):
                        model.lorentz_dict = dict([(l.name, l) for l in model['lorentz']])
                        self.init_regular_expression()
                        
                    structure = model.lorentz_dict[vertex['lorentz'][key[1]]].structure 
                    new_structure = self.lor_pattern.sub(self.simplify_lorentz,
                                                         structure)
                    lor_value = eval(new_structure % q_dict_lor)
                    if lor_value == 0:
                        new_structure = new_structure.replace('-','+')
                        lor_value = eval(new_structure % q_dict_lor)
                    lorentz_factor += abs(eval(v))**2 * lor_value**2

                apx_m *= lorentz_factor

                # If final_color contain non-singlet,
                # get the color multiplicity.
                if any([i != 1 for i in final_color]):
                    ini_color = model.get_particle(vert.get('legs')[-1].get('id')).get('color')
                    # Try to find multiplicity directly from model
                    found = False
                    try:
                        color_configs = model.color_multiplicity_def(final_color)
                        for config in color_configs:
                            if config[0] == ini_color:
                                apx_m *= config[1]
                                found = True
                                break
                        # Call the get_color_multiplicity if no suitable
                        # configs in the color_dict.
                        if not found:
                            apx_m *= self.get_color_multiplicity(ini_color,
                                                                 final_color, 
                                                                 model, True)                         
                    # Call the get_color_multiplicity if the final_color
                    # cannot be found directly in the color_dict.
                    except KeyError:
                        apx_m *= self.get_color_multiplicity(ini_color,
                                                             final_color, 
                                                             model, True)
                        
        # A quick estimate of the next-level decay of a off-shell decay
        # Consider all legs are onshell.
        else:
            M = abs(eval(ini_part.get('mass')))
            # The avg_E is lower by one more particle in the next-level.
            avg_E = (M/(len(self.get_final_legs())+1.))

            # Go through each vertex and assign factors to apx_m
            # This will take all propagators into accounts.
            # Do not run the identical vertex
            for i, vert in enumerate(self['vertices']):
                
                # Assign the value if the leg is not initial leg.
                # q is assumed as 1M
                if i < len(self.get('vertices'))-1: 
                    apx_m *= self.get_apx_fnrule(vert.get('legs')[-1].get('id'),
                                                 1*M, False, model, True)

                # Assign the value to initial particle.
                else:
                    apx_m *= self.get_apx_fnrule(vert.get('legs')[-1].get('id'),
                                                 M, True, model)

                # Evaluate the coupling strength
                vertex =  model.get('interaction_dict')[abs(vert.get('id'))]
                lorentz_factor = 0
                q_dict_lor = {}
                for key, v in vertex['couplings'].items():
                    
                    if not hasattr(model, 'lorentz_dict'):
                        model.lorentz_dict = dict([(l.name, l) for l in model['lorentz']])
                        self.init_regular_expression()
                        
                    structure = model.lorentz_dict[vertex['lorentz'][key[1]]].structure 
                    new_structure = self.lor_pattern.sub(self.simplify_lorentz,
                                                         structure)
                    for i, part in enumerate(vertex['particles']):
                        mass  = abs(eval(part.get('mass')))
                        q_dict_lor['q%i' % (i+1)] = mass / 2
                    
                    lor_value = eval(new_structure % q_dict_lor)
                    # Avoid accidental zeros in lor_value
                    if lor_value == 0:
                        new_structure = new_structure.replace('-','+')
                        lor_value = eval(new_structure % q_dict_lor)

                    lorentz_factor += abs(eval(v))**2 * eval(new_structure % q_dict_lor)**2

                apx_m *= lorentz_factor

            # Calculate the contribution from final legs
            for leg in self.get_final_legs():
                apx_m *= self.get_apx_fnrule(leg.get('id'),
                                             avg_E, True, model)

        # For both on-shell and off-shell cases,
        # Correct the factor of spin/color sum of initial particle (average it)
        apx_m *= 1./(ini_part.get('spin'))
        self['apx_matrixelement_sq'] = apx_m
        return apx_m
            
    def get_apx_fnrule(self, pid, q, onshell, model, est = False):
        """ The library that provide the 'approximated Feynmann rule'
            q is the energy of the leg. The onshell label is to decide
            whether this particle is final or intermediate particle."""
        
        part = model.get('particle_dict')[pid]
        mass  = abs(eval(part.get('mass')))

        # Set the propagator value (square is for square of matrix element)
        # The width is included in the propagator.
        if onshell:
            value = 1.
        else:
            if not est:
                value = 1./((q ** 2 - mass ** 2) ** 2 + \
                                mass **2 * part.get('apx_decaywidth') **2)
            # Rough estimation on propagator. Avoid the large propagator when
            # q is close to mass
            else:
                m_large = max([q, mass])
                value = 1./(0.5* m_large**2)**2

        # Set the value according the particle type
        # vector boson case
        if part.get('spin') == 3:
            if onshell:
                # For massive vector boson
                if mass != 0. :
                    value *= (1+ (q/mass) **2)
                # For massless boson
                else:
                    value *= 1
            # The numerator of propagator.
            else:
                value *= (1 - 2* (q/mass) **2 + (q/mass) **4)
        # fermion case
        elif part.get('spin') == 2:
            if onshell:
                value *= 2.*q 
            else:
                value *= q **2
        #spin 3/2 case
        elif part.get('spin') == 4:
            if onshell:
                #spin1 part
                # For massive vector boson
                if mass != 0. :
                    value *= (1+ (q/mass) **2)
                # For massless boson
                else:
                    value *= 1
                #spin1/2 part
                value *= 2.*q 
            # The numerator of propagator.
            else:
                # For massive 
                if mass != 0. :
                    value *= 4/9 * (q)**2 * (1-q**2/mass**2)**2
                # For massless case
                else:
                    value *= q**2
        #spin2 case
        elif part.get('spin') == 5:
            if onshell:
                #spin1**2 
                # For massive vector boson
                if mass != 0. :
                    value *= (1+ (q/mass) **2)**2
                # For massless boson
                else:
                    value *= 1
            # The numerator of propagator.
            else:
                # For massive 
                if mass != 0. :
                    value *= (7/6-4/3*q**2/mass**2+2/3*q**4/mass**4)**2
                # For massless case
                else:
                    value *= 1      
                     

        # Do nothing for scalar
        return value

    def get_color_multiplicity(self, ini_color, final_color, model, base=False):
        """ Get the color multiplicity recursively of the given final_color.
            The multiplicity is obtained based on the color_multiplicity_def
            funtion in the model.
            If the color structure of final_color matches the 
            color_multiplicity_def, return the multiplicity.
            Otherwise, return 1 for the get_color_multiplicity with base = True
            or return 0 for the get_color_multiplicity with base = False."""
            
        # Combine the last two color factor to get the possible configs.
        color_configs = model.color_multiplicity_def([final_color.pop(),
                                                      final_color.pop()])
        c_factor = 1.
        # Try each config
        for config in color_configs:
            # The recursion ends when the length of the final_color now is 0.
            # (i.e. length = 2 before pop)
            if len(final_color) == 0:
                # If the final_color is consistent to ini_color, return the
                # nonzero multiplicity
                if config[0] == ini_color:
                    return config[1]

            else:
                # If next_final_color has more than one element,            
                # creaat a new final_color for recursion.
                next_final_color = copy.copy(final_color)
                next_final_color.append(config[0])

                # Call get_color_multiplicity with next_final_color as argument.
                c_factor = config[1]* self.get_color_multiplicity(ini_color,
                                                                  next_final_color,
                                                                  model)

                # If the c_factor is not zero, the color configs match successfully.
                # Return the c_factor.
                if c_factor != 0:
                    return c_factor

        # If no configs are satisfied...
        # Raise the warning message and return 1 for base get_color_multiplicity
        if base:
            logger.warning("Color structure %s in interaction is not included!" %str(final_color))
            return 1
        # return 0 for intermediate get_color_multiplicity.
        else:
            return 0
        

    def get_apx_psarea(self, model):
        """ Calculate the approximate phase space area. For off-shell case,
            it only approximate it as if it is on-shell.
            For on-shell case, it calls the recursive calculate_apx_psarea
            to get a more accurate estimation. """

        M = abs(eval(model.get_particle(self.get_initial_id()).get('mass')))

        # Off-shell channel only estimate the psarea if next level is onshell.
        if not self.get_onshell(model):
            # The power of extra integration for this level is
            # number of current final particle -2 
            # (3-body decay -> 1 integration)
            self['apx_psarea'] = 1/(8*math.pi)*\
                pow((c_psarea*(M/8./math.pi)**2), len(self.get_final_legs())-2)

        # For onshell case and psarea has not been calculated
        elif not self['apx_psarea']:
            # The initial particle mass
            mass_list = copy.copy(self['final_mass_list'])
            self['apx_psarea'] = self.calculate_apx_psarea(M, mass_list)

        return self['apx_psarea']

    def calculate_apx_psarea(self, M, mass_list):
        """Recursive function to calculate the apx_psarea.
           For the estimation of integration, it takes the final mass in the
           mass_list. The c_psarea is corrected in each integration.
           Symmetric factor of final state is corrected.
           """

        # For more than 2-body decay, estimate the integration from the
        # middle point with c_psarea factor for correction, then
        # calls the function itself with reduced mass_list.
        if len(mass_list) >2 :
            # Mass_n is the mass that use pop out to calculate the ps area.
            mass_n = mass_list.pop()
            # Mean value of the c.m. mass of the rest particles in mass_list
            M_eff_mean = ((M-mass_n) +sum(mass_list))/2
            # The range of the c.m. mass square
            delta_M_eff_sq = ((M-mass_n) ** 2-sum(mass_list) ** 2)
            # Recursive formula for ps area,
            # initial mass is replaced as the square root of M_eff_sq_mean
            return math.sqrt((M ** 2+mass_n ** 2-M_eff_mean**2) ** 2-\
                                 (2*M *mass_n) ** 2)* \
                                 self.calculate_apx_psarea(M_eff_mean, mass_list)*\
                                 delta_M_eff_sq*c_psarea* \
                                 1./(16*(math.pi ** 2)*(M ** 2))
            
        # for two particle decay the phase space area is known.
        else:
            # calculate the symmetric factor first
            self['s_factor'] =1
            id_list = sorted([l.get('id') for l in self.get_final_legs()])
            count =1
            for i, pid in enumerate(id_list):
                if i !=0 and id_list[i-1] == pid:
                    count += 1
                elif count != 1:
                    self['s_factor'] = self['s_factor'] * math.factorial(count)
                    count = 1

            # This complete the s_factor if the idparticle is in the last part
            # of list.
            if count != 1:
                self['s_factor'] = self['s_factor'] * math.factorial(count)
            return math.sqrt((M ** 2+mass_list[0] ** 2-mass_list[1] ** 2) ** 2-\
                                 (2* M *mass_list[0]) ** 2)* \
                                 1./(8*math.pi*(M ** 2)*self['s_factor'])


    def get_apx_decaywidth(self, model):
        """Calculate the apx_decaywidth
           formula: Gamma = ps_area* matrix element square * (1/2M)
           Note that it still simulate the value for off-shell decay."""

        # Return the value now if width has been calculated.
        if self['apx_width_calculated']:
            return self['apx_decaywidth']

        # Calculate width (use mass of antiparticle of initial particle)
        self['apx_decaywidth'] = self.get_apx_matrixelement_sq(model) * \
            self.get_apx_psarea(model)/ \
            (2*abs(eval(model.get_particle(self.get_initial_id()).get('mass'))))
        self['apx_width_calculated'] = True

        return self['apx_decaywidth']

    def get_apx_decaywidth_nextlevel(self, model):
        """ Estimate the sum of all the width of the next-level channels
            it developes."""

        M = abs(eval(model.get_particle(self.get_initial_id()).get('mass')))
        m_now = sum(self.get('final_mass_list'))
        avg_E = (M/(len(self.get_final_legs())+1.))
        # Ratio is the width of next-level channels over current channel.
        err = 0.
        for leg in self.get_final_legs():
            # Use only particle not anti-particle because anti-particle has no
            # width
            part = model.get_particle(abs(leg.get('id')))
            # For legs that are possible to decay.
            if (not part.get('is_stable')) and (M-m_now+part.get('2body_massdiff')) > 0.:
                # Suppose the further decay is two-body.
                # Formula: ratio = width_of_this_leg * (M/m_leg)**(-1) *
                #                  (2 * M * 8 * pi * (c_psarea* (M/8/pi)**2)) *
                #                  1/(leg_mleg(mleg)/leg_mleg(0.5M) *
                #                  Propagator of mleg(M)
                err += (part.get('apx_decaywidth')*\
                            (M/abs(eval(part.get('mass')))) **(-1) *\
                            (c_psarea*(M **3/4/math.pi)) / \
                            (self.get_apx_fnrule(leg.get('id'), avg_E,
                                                 True, model)*\
                                 self.get_apx_fnrule(leg.get('id'), 
                                                     abs(eval(part.get('mass'))), True, model))*\
                            self.get_apx_fnrule(leg.get('id'), M,
                                                False, model, True)
                        )

        self['apx_decaywidth_nextlevel'] = self.get_apx_decaywidth(model)*err

        return self['apx_decaywidth_nextlevel']


#===============================================================================
# ChannelList: List of all possible  channels for the decay
#===============================================================================
class ChannelList(base_objects.DiagramList):
    """List of decay Channel
    """

    def is_valid_element(self, obj):
        """ Test if the object is a valid Channel for the list. """
        return isinstance(obj, Channel)


#===============================================================================
# ChannelList: List of all possible  channels for the decay
#===============================================================================
class IdentifyHelasTag(diagram_generation.DiagramTag):
    """ DiagramTag daughter class to identify helas call. Model is necessary
        for initialization.
        Note: 1.) diagrams that are identical may have different type of
                  mother particles, e.g. t > b w+ and w+ > c s~ are the same
              2.) the output diagram still has the same leg number as its
                  mother diagram. This means the diagrams within the same 
                  amplitude still have consistent leg number.
                  But the leg numbers in a std_diagram may be different
                  from the numbers in model['helascalls'], though the structure
                  is identical.
              3.) in general, the last vertex of output diagram is not initial
                  vertex; however, the initial leg can be found by looking for
                  the end leg with number=1.
    """

    @staticmethod
    def link_from_leg(leg, model):
        """Returns the end link for a leg needed to identify Helas calls
        configs: ((spin, color), (leg id, leg number, leg state) 
        N.B.: only spin and color are used for comparison."""


        part = model.get_particle(leg.get('id'))


        return [((part.get('spin'), part.get('color')),
                 (leg.get('id'), leg.get('number'), leg.get('state')))]
        
    @staticmethod
    def vertex_id_from_vertex(vertex, last_vertex, model, ninitial):
        """Returns the info needed to identify Helas calls:
        ((interaction lorentz, color, keys of couplings), ineraction id)
        N.B.: interaction id is listed for reconstruction purpose,
        but not used for comparison."""

        inter = model.get_interaction(vertex.get('id'))

        return ((inter['lorentz'], inter['color'], 
                 sorted(inter['couplings'].keys())),
                inter['id'])

    @staticmethod
    def leg_from_link(link):
        """Return a leg from a link"""

        if link.end_link:
            # This is an external leg, info in links
            return base_objects.Leg({'number':link.links[0][1][1],
                                     'id':link.links[0][1][0],
                                     'state':link.links[0][1][2],
                                     'onshell':False})

        # This shouldn't happen
        assert False

    @staticmethod
    def id_from_vertex_id(vertex_id):
        """Return the numerical vertex id from a link.vertex_id"""
        return vertex_id[1]


#===============================================================================
# DecayAmplitude: An Amplitude like object contain Process and Channels
#===============================================================================
class DecayAmplitude(diagram_generation.Amplitude):
    """ DecayAmplitude is derived from Amplitude. It collects channels 
        with the same final states and create a Process object to describe it.
        This could be used to generate HELAS amplitude."""

    sorted_keys = ['process', 'diagrams', 'apx_decaywidth', 'apx_br',
                   'exa_decaywidth', 'part_sn_dict', 'inter_sn_dict',
                   'ab2real_dicts', 'has_mirror_process']#,
                   #'potential_gauge_dependence']

    def default_setup(self):
        """Default values for all properties. Property 'diagrams' is now
           as ChannelList object."""

        self['process'] = base_objects.Process()
        self['diagrams'] = ChannelList()
        self['apx_decaywidth'] = 0.
        self['apx_br'] = 0.
        self['exa_decaywidth'] = False
        #self['potential_gauge_dependence'] = False

        # Properties used in abstract amplitude
        self['part_sn_dict'] = {}
        self['inter_sn_dict'] = {}
        self['ab2real_dicts'] = Ab2RealDictList()
        # add has_mirror_process
        self['has_mirror_process'] = False

    def __init__(self, argument=None, model=None):
        """ Allow initialization with a Channel and DecayModel to create 
            the corresponding process."""

        if isinstance(argument, Channel) and isinstance(model, DecayModel):

            super(DecayAmplitude, self).__init__()

            # Set the corresponding process.
            self.set_process(argument, model)
            
            # Set diagram
            self.set('diagrams', ChannelList([argument]))
            #self['potential_gauge_dependence'] = argument['potential_gauge_dependence']

        else:
            super(DecayAmplitude, self).__init__(argument)

    def filter(self, name, value):
        """Filter for valid amplitude property values."""

        if name == 'process':
            if not isinstance(value, base_objects.Process):
                raise self.PhysicsObjectError, \
                        "%s is not a valid Process object." % str(value)
            # Reset the width and br
            self.reset_width_br()

        if name == 'diagrams':
            if not isinstance(value, ChannelList):
                raise self.PhysicsObjectError, \
                        "%s is not a valid ChannelList object." % str(value)
            # Reset the width and br
            self.reset_width_br()

        if name == 'apx_decaywidth' and name == 'apx_br':
            if not isinstance(value, float):
                raise self.PhysicsObjectError, \
                        "%s is not a float." % str(value)

        if name == 'exa_decaywidth':
            if not isinstance(value, float) and not isinstance(value,bool):
                raise self.PhysicsObjectError, \
                        "%s is not a float." % str(value)

        return True

    def get(self, name):
        """Get the value of the property name."""

        # When apx_decaywidth is requested, recalculate it if needed.
        # Calculate br in the same time.
        if name == 'apx_decaywidth' and not self[name]:
            model = self['process']['model']
            self['apx_decaywidth'] = sum([c.get('apx_decaywidth') \
                                  for c in self['diagrams']])

        # If apx_br is requested, recalculate from the apx_decaywidth if needed.
        if name == 'apx_br' and not self[name]:
            # diagrams should have already save their initial id.
            ini_part = self['process']['model'].\
                get_particle(self['diagrams'][0].get_initial_id())
            try:
                self['apx_br'] = self.get('apx_decaywidth')/ \
                    ini_part['apx_decaywidth']
            except ZeroDivisionError:
                pass
                #logger.warning("Try to get branch ratio from a zero width particle %s. No action proceed." % ini_part.get('name'))

        return super(DecayAmplitude, self).get(name)

    def get_sorted_keys(self):
        """Return DecayProcess property names as a nicely sorted list."""

        return self.sorted_keys

    def set_process(self, dia, model):
        """ Setup the process from the diagram,
            used in initial setup."""

        # Check the number of initial leg is 1
        if dia['vertices'][-1]['legs'][-1]['number'] != 1:
            raise self.PhysicsObjectError, \
                "The number of initial leg should be setted as 1."

        # Append the initial leg.
        leglist = base_objects.LegList([base_objects.Leg({'id': dia.get_initial_id(model),
                                                          'number':1,
                                                          'state': False})])
        # Extract legs from final legs of Channel.
        leglist.extend(base_objects.LegList(\
                copy.deepcopy(sorted([l for l in dia.get_final_legs()], 
                                     legcmp_bynumber))))
            
        # Set up process and model.
        self.set('process', base_objects.Process({'legs':leglist}))
        self['process'].set('model', model)

    def add_std_diagram(self, new_dia, model=None):
        """ Add new diagram into amplitude
            and check if the number identifiers 
            of outgoing legs are consistent with the process."""

        if not isinstance(new_dia, Channel):
            raise self.PhysicsObjectError,\
                "The argument should be Channel object."

        # Model must be provided in argument or in self['process']
        if not model:
            if not self['process']['model']:
                raise self.PhysicsObjectError,\
                    "No model is given either in argument or in process."
            else:
                model = self['process']['model']

        # If this amplitude has no process
        if not self['process']['legs']:
            self.set_process(new_dia, model)
            self['diagrams'].append(new_dia)
            return
        
        # Edit the potential_gauge_dependence property
#        if new_dia['potential_gauge_dependence']:
#            self['potential_gauge_dependence'] = True


        # non_std_number: number of new_dia
        # take final legs first
        non_std_numbers = [(l.get('id'),l.get('number')) \
                               for l in new_dia.get_final_legs()]


        # initial leg
        non_std_numbers.append((new_dia.get_initial_id(model), 1))
        non_std_numbers.sort(id_num_cmp)

        # std_number: numbers of legs in process
        std_numbers = [(l.get('id'),l.get('number')) \
                           for l in sorted(self['process']['legs'])]
        std_numbers.sort(id_num_cmp)

        # Return if the numbers in diagram is the same as process
        if non_std_numbers == std_numbers:
            self['diagrams'].append(new_dia)
            return

        # Conversion from non_std_number to std_number
        converted_dict = dict([(num[1], std_numbers[i][1])\
                                   for i, num in enumerate(non_std_numbers)])
        
        # 1st stage of converting all legs: change numbering without fixing
        # wrong number flows (e.g. number 3 2 > 3)
        all_numbers_goal = []
        for vert in new_dia.get('vertices'):
            for leg in vert.get('legs'):
                leg['number'] = converted_dict[leg['number']]
                all_numbers_goal.append((leg['id'], leg['number']))

        # 2nd stage of converting all legs: fixing illegal number flows.
        # (except for the first one)
        for pos,vert in enumerate(new_dia.get('vertices')[:-1]):
            lowest_num = vert.get('legs')[0]['number']
            for leg in vert.get('legs')[:-1]:
                if leg['number'] < lowest_num:
                    lowest_num = leg['number']

            mother_leg = vert.get('legs')[-1]
            old_id_number = (mother_leg['id'], mother_leg['number'])
            if old_id_number[1] != lowest_num:
                # Change the number of mother
                mother_leg['number'] = lowest_num                
                
                # Find the leg associated with the mother of this vertex,
                # and change its number to lowest number
                found = False
                for pre_vert in new_dia.get('vertices')[pos:]:
                    for child_leg in pre_vert['legs'][:-1]:
                        if (child_leg['id'], child_leg['number']) == old_id_number:
                            child_leg['number'] = lowest_num
                            found = True
                            break
                    if found:
                        break
            """
            # Put particles before antiparticles            
            vert['legs'][:-1] = sorted(vert['legs'][:-1],
                                       key=lambda leg: leg['id'],
                                       reverse = True)"""

        new_dia.initial_setups(model, True)

        # Add this standard diagram into diagrams
        self['diagrams'].append(new_dia)


    def reset_width_br(self):
        """ Reset the value of decay width and branch ratio.
            Automatically done in the set(filter) function.
            This is needed when content of diagrams or process are changed."""
        
        self['apx_decaywidth'] = 0.
        self['apx_br'] = 0.

    def nice_string(self, indent=0):
        """Returns a nicely formatted string of the amplitude content"""
#        if not self['potential_gauge_dependence']:
#            return self.get('process').nice_string(indent) + "\n" + \
#                self.get('diagrams').nice_string(indent)
#        else:
        return self.get('process').nice_string(indent) + ", *PGD" + "\n" + \
                self.get('diagrams').nice_string(indent)
            

    def decaytable_string(self, format='normal'):
        """ Write the string in the format for decay table.
            format = 'normal': show only branching ratio
                   = 'full'  : show branching ratio and all the channels."""

        output='   %.5e   %d' %(self.get('apx_br'),
                                len(self['process']['legs'])-1)
        output += ''.join(['%11d' %leg.get('id') \
                               for leg in self['process']['legs'][1:]])

        if format == 'cmp':
            if self.get('exa_decaywidth'):
                output += '\t%4.2f' % (self.get('apx_decaywidth')/self.get('exa_decaywidth'))
            else:
                output += '\tN/A'
        
        output += '   #Br(%s)' %self.get('process').input_string()

        # Gauge dependence tag
        #if self['potential_gauge_dependence']:
        #    output += ', *PGD\n'
        #else:
        output += '\n'
        
        # Output the channels if format is full
        if format=='full':
            # Set indent of the beginning for channels
            output += self.get('diagrams').nice_string(6)

        # Return process only, get rid off the final \n
        elif format == 'normal' or format == 'cmp':
            return output[:-1]

        # Raise error if format is wrong
        else:
            raise self.PhysicsObjectError,\
                "Format %s must be \'normal\' or \'full\'." % str(format)

        # Return output for full format case.
        return output


    def abstract_nice_string(self):
        """ Print the nice string of abstract amplitudes,
        which shows the real process of this type."""

        output = super(DecayAmplitude, self).nice_string()
        output += "\nReal process:  "
        output += ', '.join([ref['process'].input_string() for ref in self['ab2real_dicts']])

        return output


    def generate_variables_dicts(self, real_amp, ab2real_id = -1):
        """ Generate the dictionary of abstract to real variables,
        default is to use the latest Ab2RealDict."""

        ab_model = self['process']['model']
        real_model = real_amp['process']['model']
        ab2realdict = self['ab2real_dicts'][ab2real_id]

        # Set the initial particle dict
        ab2realdict['mass_dict']\
            [ab_model.get_particle(self['process'].get_initial_ids()[0])['mass']] = real_model.get_particle(real_amp['process'].get_initial_ids()[0])['mass']

        # Set the final particle dict
        for ab_pid, real_pid in ab2realdict['final_legs_dict'].items():
            if isinstance(real_pid, int):
                ab2realdict['mass_dict'][ab_model.get_particle(ab_pid)['mass']] = real_model.get_particle(real_pid)['mass']
            # For the case of identical particle
            else:
                for pid in real_pid:
                    ab2realdict['mass_dict'][ab_model.get_particle(ab_pid)['mass']] = real_model.get_particle(pid)['mass']

        # Set the interaction id dict and intermediate particles
        for  ab_dia_id, real_dia_id in ab2realdict['dia_sn_dict'].items():

            ab_dia = self['diagrams'][ab_dia_id]
            real_dia = real_amp['diagrams'][real_dia_id]
            
            # Set the intermediate particle and interaction,
            # except for the identical vertex 
            # (the initial particle will be resetted, but it should be fine.)
            for i, v in enumerate(ab_dia['vertices']):

                # Set intermediate particle
                ab2realdict['mass_dict']\
                    [ab_model.get_particle(v['legs'][-1]['id'])['mass']] \
                    = real_model.get_particle(real_dia['vertices'][i]['legs'][-1]['id'])['mass']
                
                # Set the interaction
                for ab_key, real_coup in \
                        ab_model['interaction_coupling_dict'][real_dia['vertices'][i]['id']].items():

                    ab_coup = ab_model.get_interaction(ab_dia['vertices'][i]['id'])['couplings'][ab_key]
                    ab2realdict['coup_dict'][ab_coup] = real_coup
        
#===============================================================================
# DecayAmplitudeList: An Amplitude like object contain Process and Channels
#===============================================================================
class DecayAmplitudeList(diagram_generation.AmplitudeList):
    """ List for DecayAmplitudeList objects.
    """

    def is_valid_element(self, obj):
        """ Test if the object is a valid DecayAmplitude for the list. """
        return isinstance(obj, DecayAmplitude)

    def get_amplitude(self, final_ids):
        """ Get the amplitudes with the given final particles if
        exist in this list. Otherwise, return None.
        Note: Use stored finallist in Channel!"""
        
        for amp in self:
            if sorted([abs(l.get('id')) for l in amp['diagrams'][0]\
                           .get_final_legs()]) == sorted(list_abs(final_ids)):
                return amp
            

        return None

    def nice_string(self):
        """ Nice string from Amplitude """
        mystr = '\n'.join([a.nice_string() for a in self])

        return mystr

    def abstract_nice_string(self):
        """ Nice string for abstract amplitude. """
        mystr = '\n'.join([a.abstract_nice_string() for a in self])

        return mystr


    def decaytable_string(self, format='normal'):
        """ Decaytable string for Amplitudes """

        # Get the level (n-body) of this amplitudelist and
        # the total width inside this DecayAmplitudeList
        level = len(self[0].get('diagrams')[0].get_final_legs())
        total_width = sum([amp.get('apx_decaywidth') for amp in self])

        # Print the header with the total width
        mystr = '## Contribution to total width of %d-body decay: %.5e \n' \
            %(level, total_width)
        mystr += '\n'.join([a.decaytable_string(format) for a in self])

        return mystr

#===============================================================================
# AbstractModel: A Model class in which objects are abstract
#===============================================================================
class AbstractModel(base_objects.Model):
    """ Model in which particles are different only in Lorentz structure
    and color. Interactions may be duplicated in order to generate amplitude
    with the same Lorentz structure vertices or diagrams
    """
    sorted_keys = ['name', 'particles', 'parameters', 'interactions', 
                   'couplings', 'lorentz',
                   'abstract_particles_dict', 'abstract_interactions_dict',
                   'particle_type_dict', 
                   'interaction_type_dict', 'interaction_coupling_dict',
                   'ab_matrix_elements'
                  ]

    def default_setup(self):
        """The particles is changed to ParticleList"""

        super(AbstractModel, self).default_setup()
        # Other properties
        self['particles'] = DecayParticleList()
        # Object type -> Abstract Object list of the given type
        self['abstract_particles_dict'] = {}
        self['abstract_interactions_dict'] = {}

        # Real interaction id -> (Abstract type (spin, color, self-anti),
        #                         pseudo_pdg_code)
        # keys are always positive!!!
        self['particle_type_dict'] = {}

        # Real interaction id -> Abstract type (pseudo abstract particlelist,
        #                                       lorentz,
        #                                       color)
        self['interaction_type_dict'] = {}

        self['interaction_coupling_dict'] = {}

        # Store all the matrix elements
        self['ab_matrix_elements'] = AbstractHelasMatrixElementList()

        self.spin_text_dict = {1:'S', 2: 'F', 3:'V', 5:'T'}
        

    def get_sorted_keys(self):
        return self.sorted_keys

    def get_particle_type(self, part, get_code=False):
        """ Return the tuple (spin, color, self_antipart) 
        of the given particle.
        NOTE: bosons are always treated as self-antipart. """

        if isinstance(part, base_objects.Particle):
            if get_code:
                return abs(self.get_particle_type_code(part))

            if part.is_boson():
                self_antipart = True
            else:
                self_antipart = part['self_antipart']

            return (part['spin'], part['color'], self_antipart)

        # Use particle_type_dict if input is pdg_code
        elif isinstance(part, int):
            if not get_code:
                return self['particle_type_dict'][abs(part)][0]
            else:
                return self['particle_type_dict'][abs(part)][1]
        else:
            raise self.PhysicsObjectError, \
                "Input should be Particle type or pdg_code(int)."


    def get_particle_type_code(self, part):
        """ Return the pseudo pdg_code :
        (9900000+1000*spin+100*color) for non-self-antipart;
        (9910000+1000*spin+100*color) for self-antipart (Majorana fermion)
        of the given particle."""

        # Special pdg_code for Majorana fermions
        if part.is_fermion():
            if part['self_antipart']:
                return 9910000+1000*part['spin']+100*part['color']
            else:
                return int(math.copysign(
                        9900000+1000*part['spin']+100*part['color'],
                        part.get_pdg_code()))
        else:
            return 9900000+1000*part['spin']+100*part['color']


    def add_ab_particle(self, pdg_code, force=False):
        """ Functions to add new particle according to the given particle
        spin and color information. This function will assign the correct
        serial number of new particle and change the abstract_particles."""

        # Check argument type
        if not force and not isinstance(pdg_code, int):
            raise self.PhysicsObjectError,\
                "Argument must be a pdg_code."

        # To ensure pdg_code is positive
        pdg_code = abs(pdg_code)

        # Setup new particle
        ab_part = DecayParticle()
        ab_part['spin'] = self.get_particle_type(pdg_code)[0]
        ab_part['color'] = self.get_particle_type(pdg_code)[1]
        ab_part['self_antipart'] = self.get_particle_type(pdg_code)[2]

        # Check the current abstract_particles
        if self.get_particle_type(pdg_code) \
                in self['abstract_particles_dict'].keys():
            # For existing type, record the serial number
            sn = len(self['abstract_particles_dict']\
                         [self.get_particle_type(pdg_code)])
        else:
            # Setup new type of abstract particle
            self['abstract_particles_dict']\
                [self.get_particle_type(pdg_code)]= DecayParticleList()
            sn = 0
            

        # Set other properties: name = text(spin)+color
        try:
            name_string = self.spin_text_dict[ab_part.get('spin')]
        except KeyError:
            name_string = 'P'

        ab_part['name']= \
            name_string +\
            str(ab_part.get('color')) +\
            '_%02d' %sn

        # mass = M'S or N''spin''color' 
        #                 % S for self-antipart, N for non-self-antipart
        mass_string = 'M'
        if ab_part['self_antipart']:
            mass_string = mass_string + 'S'
        else:
            mass_string = mass_string + 'N'
        mass_string = mass_string + name_string +\
            str(ab_part.get('color')) +\
            '_%02d' %sn
        ab_part['mass'] = mass_string

        # pdg_code = 99'1:S ,0: N'00'spin''color'
        #  e.g. w+: 9903100, gluino: 9912300
        ab_part['pdg_code'] = self.get_particle_type(pdg_code, get_code=True)+sn

        # Set anti-part properties if this is not a self-antiparticle
        if not ab_part['self_antipart']:
            ab_part['antiname'] = ab_part.get('name')+'~'

            # Create the anti-part
            anti_ab_part = DecayParticle(ab_part)
            anti_ab_part['is_part'] = False
        
        # Append ab_part into self['particles'] and abstract_particles_dict,
        self['particles'].append(ab_part)
        self['abstract_particles_dict'][self.get_particle_type(ab_part)].append(ab_part)
        logger.info("Add Abstract Particle %s according to #%s" \
                        %(ab_part.get('name'), pdg_code))


        # Reset the particle dictionary
        self.get('particle_dict')
        self['particle_dict'][ab_part['pdg_code']] = ab_part
        if not ab_part['self_antipart']:
            self['particle_dict'][-ab_part['pdg_code']] = anti_ab_part



    def setup_particle(self, part, force=False):
        """Add real particle into AbstractModel, 
        convert into abstract particle."""

        # Check argument type
        if not force and not isinstance(part, base_objects.Particle):
            raise self.PhysicsObjectError,\
                "Argument must be a Particle object."

        # Setup the particle_type_dict for all particle
        self['particle_type_dict'][part['pdg_code']] = \
            (self.get_particle_type(part), 
             self.get_particle_type(part, get_code=True))

        if self.get_particle_type(part) \
                not in self['abstract_particles_dict'].keys():
            # if not found in existing particle, create a new abstract particle
            self.add_ab_particle(part['pdg_code'], force)


    def setup_particles(self, part_list, force=False):
        """Add real particles into AbstractModel, 
        convert into abstract particles"""

        # Check argument type
        if not force and not isinstance(part_list, base_objects.ParticleList):
            raise self.PhysicsObjectError,\
                "Argument must be a ParticleList."

        # Call add_particle for each particle
        for part in part_list:
            self.setup_particle(part, force)

        # Reset dictionaries in the model
        self.reset_dictionaries()
        


    def get_particlelist_type(self, pdgcode_list, ignore_dup=True,
                              sn_dict={}, std_output=True):
        """ Return a list of the type of the given particlelist,
        and a dictionary records the number of each type of particle.
        The pdg_code will start from the number given by serial_number_dict.
        Note1: If ignore_dup is True,
               the same particles will not be assigned with the same id.
        Note2: std_output is used to get the same abstract_pidlist
               for a given real_pids, no matter what order is.
               Usage e.g.: pids of final_legs in one real amplitude.
               Not used for: comparing two real_pids.
        """

        pseudo_ab_particlelist = []
        # Used the given sn_dict or an empty one
        serial_number_dict = copy.copy(sn_dict)

        if not pdgcode_list:
            return pseudo_ab_particlelist, serial_number_dict
        
        # If standard output is required, input_pdgcode_list
        # is sorted.
        input_pdgcode_list = copy.copy(pdgcode_list)
        if std_output:
            # List of indices to record the original order
            original_order = range(len(pdgcode_list))

            # Sort according to (particle_type, real_pid) of input_pdgcode_list
            temp = zip(input_pdgcode_list, original_order)
            temp.sort(key=lambda item: \
                          (self.get_particle_type(item[0]),item[0]))
            try:
                input_pdgcode_list, original_order = zip(*temp)
            except:
                print original_order, pdgcode_list

        # Construct the abstract pid list
        for i, pdg_code in enumerate(input_pdgcode_list):

            # Use get_particle_type function
            part_type = self.get_particle_type(abs(pdg_code))
            # The sign of part_type_code depends on whether the
            # abstract particle is self-antipart or not
            part_type_code = self.get_particle_type(abs(pdg_code), True)
            if pdg_code < 0 and not part_type[2]:
                part_type_code = -part_type_code

            # Add the pseudo pdg code, if it is the first one in its type
            # appearing in the list.
            # There should be a abstract particle for this type,
            # if the setup_particles has been run correctly.
            if not part_type in serial_number_dict.keys():
                pseudo_ab_particlelist.append(\
                        part_type_code)

                serial_number_dict[part_type] = 1

            # If not ignore duplicate particles,
            # check if there are same particle in the list.
            else:
                set_new = True

                if not ignore_dup:
                    for j, previous_pdgcode in enumerate(input_pdgcode_list):
                        # find duplicate particle, use the abstract particle
                        # already exists.
                        # No need to update the serial_number_dict
                        if j == i:
                            break

                        if abs(previous_pdgcode) == \
                                abs(pdg_code):
                            pseudo_ab_particlelist.append(\
                                int(math.copysign(\
                                        pseudo_ab_particlelist[j],
                                        part_type_code)))
                                
                            set_new = False
                            break

                # Append new abstract particle if not appears before
                if set_new:
                    # If the abstract particle of the given serial number
                    # does not exist, add a new abstract particle and append
                    # it.
                    if serial_number_dict[part_type] >= \
                            len(self['abstract_particles_dict'][part_type]):
                        self.add_ab_particle(pdg_code, True)
                            
                    # Append the pdg_code into the list,
                    # starting from the s/n given by serial_number_dict
                    pseudo_ab_particlelist.append(\
                        int(math.copysign(\
                                abs(part_type_code)+\
                                    serial_number_dict[part_type],
                                part_type_code))
                        )

                    # Update the serial_number_dict
                    serial_number_dict[part_type] += 1

        # For std_output, reorder the pseudo_ab_particlelist
        # to be the same as pdgcode_list
        if std_output:
            temp = zip(original_order, pseudo_ab_particlelist)
            temp.sort()
            original_order, pseudo_ab_particlelist = zip(*temp)
            pseudo_ab_particlelist = list(pseudo_ab_particlelist)

        return pseudo_ab_particlelist, serial_number_dict


    def get_color_string(self, inter):
        """ Return the correct color string according to the 
        sorted order of particle. """
        
        pass


    def get_interaction_type(self, inter_id):
        """ Return the tuple (lorentz, particles_type) of the given interaction
        Raise error if type is not in quick reference dictionary."""

        # Check the quick reference dictionary.
        # If the lorentz type has been imported, use the type
        # provided by the dictionary. Otherwise, raise error.
        # Note: the keys, (lorentz, part_types), are well-sorted in
        # setup_interactions.
        try:
            return self['interaction_type_dict'][inter_id]
        except KeyError:
            raise self.PhysicsObjectError, \
                "Interaction #%d has not been imported into AbstractModel." \
                % inter_id


    def add_ab_interaction(self, inter_id, force=False, color = None):
        """ Functions to set new interaction according to the given interaction
        particles and lorentz sturcture. This function will assign the correct
        serial number of new interaction, change the abstract_interactions.
        Note: Add NEW type of interaction should use setup."""

        # Check argument type
        if not force and not isinstance(inter_id, int):
            raise self.PhysicsObjectError,\
                "Argument must be an Interaction id."
        
        # Setup new interaction
        ab_inter = base_objects.Interaction()
        inter_type = self.get_interaction_type(inter_id)

        # Check the current abstract_interactions
        if inter_type in \
                self['abstract_interactions_dict'].keys():
            # For existing type, record the serial number
            sn = len(self['abstract_interactions_dict'][inter_type])
        else:
            # Setup the new item, serial number = 0
            self['abstract_interactions_dict'][inter_type] = \
                base_objects.InteractionList()
            sn = 0


        # Type_sn is the serial number of the type
        if sn == 0:
            type_sn = len(self['abstract_interactions_dict'].keys())
        else:
            type_sn = \
                self['abstract_interactions_dict'][inter_type][0]['id']/1000

        # id = ___0__
        #        |  |_> the serial number
        #        | 
        #        |_> The serial number of the interaction type
        ab_inter['id'] = 1000*type_sn + sn

        # Get particle from the model.particle_dict
        # to ensure that the particle is DecayParticle, not Particle
        ab_inter['particles'] = DecayParticleList(\
            [self.get_particle(pid) for pid in inter_type[0]])
        ab_inter['lorentz'] = list(inter_type[1])

        # Retrieve the color information
        if color:
            ab_inter['color'] = color
        elif sn != 0:
            # Use the information before
            ab_inter['color'] = \
                self['abstract_interactions_dict'][inter_type][0]['color']
        else:
            raise self.PhysicsObjectError,\
                "Error to add interaction. No color information available."

        # couplings = G###cl##
        #              |  || |-> the serial number
        #              |  ||-> the lorentz identifier
        #              |  |-> the color identifier
        #              |-> The serial number of the interaction type
        ab_inter['couplings'] = {}
        for i, colr in enumerate(ab_inter['color']):
            for j, lorentz in enumerate(ab_inter['lorentz']):
                ab_inter['couplings'][(i, j)] = 'G%03d%1d%1d%02d' \
                    %(type_sn, i, j, sn)


            
        # Append ab_inter into self['interactions'] and 
        # abstract_interactions_dict
        self['interactions'].append(ab_inter)
        self['abstract_interactions_dict'][inter_type].append(ab_inter)

        # Reset the dictionary
        self['interaction_dict'][ab_inter['id']] = ab_inter
        

    def setup_interactions(self, inter_list, anti_dict, force=False):
        """Add real interactions into AbstractModel, 
        convert into abstract interactions.
        The lorentz and color 
        structures keep to be the union of all correpsonding
        real ones. 
        Construct the quick reference dictionary,
        and setup the coupling constants in the end."""

        # Check argument type
        if not force and not isinstance(inter_list, 
                                        base_objects.InteractionList):
            raise self.PhysicsObjectError,\
                "Argument must be an InteractionList."

        # Add all interactions except for
        # 1. radiation interaction, 2. interaction with all stable particles
        for inter in inter_list:
            # Check validity: remove stable particles
            final_list = [p.get('pdg_code')\
                              for p in inter['particles'] \
                              if not p['is_stable'] ]

            # Check validity: do not count duplicate particles,
            # i.e. those cannot decay through this interaction.
            __duplicate__ = True
            for pid in final_list:
                if final_list.count(pid) == 1:
                    __duplicate__ = False
                    break

            # Ignore all stable particle interaction (e.g. g g > g g) (len=0)
            # and interaction with all duplicate particles
            if len(final_list) == 0 or __duplicate__:
                continue


            # Use the get_particlelist_type to find the particles type of this
            # interaction.
            # NOTE: particles do not sort,
            #       ambiguity might happened!!
            abpart_types, sn = \
                self.get_particlelist_type([p.get_pdg_code() for p in\
                                                inter['particles']])

            # The key have to be sorted
            new_key = [tuple(abpart_types),
                       tuple(sorted(inter['lorentz'])),
                       tuple(sorted([str(colorstring) \
                                         for colorstring in inter['color']]))
                       ]
            color_list = inter['color']
            remove_list = []
            is_new_key = True

            # Check if this type is the subset of known lorentz type,
            # or contain the known lorentz type.
            # If so, union the two and continue search.
            # Do the same merge for color structure
            # Note: the abstract_interactions_dict always maintains to have
            # its keys that have no intersect lorentz with others.
            for key, ab_interlist in self['abstract_interactions_dict'].items():
                # Check lorentz if the particles are in the same types
                if key[0] == new_key[0]:
                    # Stop if this lorentz and color are already in
                    # included
                    if set(key[1]).issuperset(new_key[1]) and \
                            set(key[2]).issuperset(new_key[2]):
                        is_new_key = False
                        break
                    # Continue if no joit between the lorentz types
                    # and color strings either.
                    elif set(key[1]).isdisjoint(new_key[1]) and \
                            set(key[2]).isdisjoint(new_key[2]):
                        continue
                    # This key has a related lorentz structure or
                    # color string.
                    # Create the unions and prepare to remove this key
                    else:
                        union_lorentz = \
                            tuple(sorted(set(key[1]).union(inter['lorentz'])))
                        union_color = \
                            tuple(sorted(set(key[2]).union(new_key[2])))
                        new_key[1] = union_lorentz
                        new_key[2] = union_color

                        # Remember to remove this key later
                        remove_list.append(key)

                        # Collect the real color string
                        color_list.extend([c for c in ab_interlist[0]['color']\
                                               if not c in color_list])

            # Set temporate interaction_type_dict,
            # transform list into tuple in order to be eligible for key.
            self['interaction_type_dict'][inter['id']] = tuple(new_key)

            # If it is a new key, add interaction
            # in abstract_interactions_dict
            if is_new_key:
                
                # Use add_ab_interaction to get the correct format,
                # it will find type from interaction_type_dict
                self.add_ab_interaction(inter['id'], color = color_list)

                # Remove subset, if remove_list is not empty
                for remove_key in remove_list:
                    # Remove old abstract interactions in interaction list
                    for old_int in self['abstract_interactions_dict'][remove_key]:
                        self['interactions'].remove(old_int)
                    del self['abstract_interactions_dict'][remove_key]
                
        # Reset the id of all abstract interactions
        # (the deletion could cause some errors.)
        for i, ab_inter in enumerate(self['interactions']):
            
            type_sn = i+1
            # id = ___0__
            #        |  |_> the serial number
            #        | 
            #        |_> The serial number of the interaction type
            ab_inter['id'] = 1000*type_sn

            # couplings = G_______
            #              |  || |-> the serial number
            #              |  ||-> the lorentz identifier
            #              |  |-> the color identifier
            #              |-> The serial number of the interaction type
            ab_inter['couplings'] = {}
            for i, colr in enumerate(ab_inter['color']):
                for j, lorentz in enumerate(ab_inter['lorentz']):
                    ab_inter['couplings'][(i, j)] = 'G%03d%1d%1d00' \
                        %(type_sn, i, j)



        # Update the quick reference dict
        # and setup the interaction_coupling_dict
        for inter in inter_list:

            # Ignore interactions that cannot decay
            if inter['id'] not in self['interaction_type_dict'].keys():
                continue

            # Update the type dict
            old_type = self['interaction_type_dict'][inter['id']]
            for new_type in self['abstract_interactions_dict'].keys():
                if new_type[0] == old_type[0] and \
                        (set(old_type[1]) < set(new_type[1]) or \
                             set(old_type[2]) < set(new_type[2])):

                        self['interaction_type_dict'][inter['id']] = new_type


            # Construct the coupling dict
            self['interaction_coupling_dict'][inter['id']] = {}
            inter_type = self['interaction_type_dict'][inter['id']]            
            ab_inter = self['abstract_interactions_dict'][inter_type][0]
            for key, coup in inter['couplings'].items():
                color = inter['color'][key[0]]
                lorentz = inter['lorentz'][key[1]]
                ab_key = [0, 0]
                # Get new key for the coupling
                ab_key[0] = ab_inter['color'].index(color)
                ab_key[1] = ab_inter['lorentz'].index(lorentz)
                
                self['interaction_coupling_dict'][inter['id']][tuple(ab_key)]\
                    = coup


        # Update dict for anti-iteraction
        for inter in inter_list:
            # For possible anti-interaction, 
            # Update the interaction_type_dict, interaction_coupling_dict
            if inter['id'] in anti_dict.keys():
                    anti_inter_id = anti_dict[inter['id']]                    
            else:
                continue

            # property of -(inter_id) -> property of anti_inter_id
            self['interaction_type_dict'][-inter['id']] = \
                self['interaction_type_dict'][anti_inter_id]
            self['interaction_coupling_dict'][-inter['id']] = \
                self['interaction_coupling_dict'][anti_inter_id]


        # Reset dictionaries in the model
        self.reset_dictionaries()
        self.get('particle_dict')
        self.get('interaction_dict')


    def get_interactionlist_type(self, interid_list, ignore_dup=False,
                                 sn_dict={}):
        """ Return a list of the type of the given interactions,
        and a dictionary records the number of each type of interaction.
        The abstract interaction id will start from the number given 
        by serial_number_dict.
        Note: If ignore_dup is True,
        the same interactions will not assign the same id.
        """

        pseudo_ab_interlist = []
        # Used the given sn_dict or an empty one
        serial_number_dict = copy.copy(sn_dict)
        
        for i, inter_id in enumerate(interid_list):

            # Append identity vertices
            # Not being used by current standard of diagrams
            if inter_id == 0:
                pseudo_ab_interlist.append(0)
                continue

            # Use get_interaction_type function
            inter_type = self.get_interaction_type(inter_id)
            inter_type_code = \
                self['abstract_interactions_dict'][inter_type][0]['id']

            # Add the pseudo id, if it is the first one in its type
            # appearing in the list.
            # There should be a abstract interaction for this type,
            # if the setup_interactions has been run correctly.
            if not inter_type in serial_number_dict.keys():
                pseudo_ab_interlist.append(\
                        inter_type_code)

                serial_number_dict[inter_type] = 1

            # If not ignore duplicate interactions (default),
            # check if there are same interactions in the list.
            else:
                set_new = True

                if not ignore_dup:
                    for j, previous_id in enumerate(interid_list):
                        # find duplicate interaction, 
                        # use the abstract interaction already exists.
                        # No need to update the serial_number_dict
                        if j == i:
                            break

                        if previous_id == \
                                inter_id:
                            pseudo_ab_interlist.append(\
                                pseudo_ab_interlist[j])

                            set_new = False
                            break

                # Append new abstract interaction if not appears before
                if set_new:
                    # If the abstract interaction of the given serial number
                    # does not exist, add a new abstract interaction and append
                    # it.
                    if serial_number_dict[inter_type] >= \
                            len(self['abstract_interactions_dict'][inter_type]):
                        self.add_ab_interaction(inter_id, True)
                            
                    # Append the pdg_code into the list,
                    # starting from the s/n given by serial_number_dict
                    pseudo_ab_interlist.append(\
                        inter_type_code+\
                            serial_number_dict[inter_type]
                        )

                    # Update the serial_number_dict
                    serial_number_dict[inter_type] += 1
                        
                    
        return pseudo_ab_interlist, serial_number_dict


    # Helper function to construct the AbstractMatrixElement
    def compare_diagrams(self, ab_dia, real_dia, ab2realdict=None):
        """ Return True if the two diagrams are in the same abstract type.
        The ab_dia and real_dia must have the same topology in this algorithm.
        Algorithm: 
        a. Compare the pseudo-abstract interaction id list by ORDER
        b. Compare the pseudo-abstract pdg_code list by ORDER.
        """

        if not isinstance(ab_dia, Channel) or not isinstance(real_dia, Channel):
            raise self.PhysicsObjectError,\
                "The first two argument are not Channel objects."
        if ab2realdict != None and not isinstance(ab2realdict, Ab2RealDict):
            raise self.PhysicsObjectError,\
                "The final argument should be Ab2RealDict, otherwise should be omitted."
            

        # Interaction id list
        ab_inter_id_list = ab_dia['abstract_type'][0]
        real_inter_id_list = [v.get('id') for v in real_dia['vertices']]

        # Quick compare from the length.
        if len(ab_inter_id_list) != len(real_inter_id_list):
            return False

        # Full comparision of interaction type.
        # The duplicated intereactions are taken into consideration.
        if self.get_interactionlist_type(real_inter_id_list)[0] != \
                ab_inter_id_list:
            return False

        # Intermediate Particle id list
        real_pdgcode_list = [v.get('legs')[-1]['id'] \
                                 for v in real_dia['vertices'][:-1]]
        

        # Full comparision of intermediate particle type.
        # 1. Duplicated particles will be treated as different.
        # 2. Do not get the std_output. We want to compare two real pids,
        #    not to 
        if self.get_particlelist_type(real_pdgcode_list, std_output=False)[0]\
                != ab_dia['abstract_type'][1]:
            return False

        # Continue to check the final Particle id list
        ini_type = self.get_particle_type(real_dia.get_initial_id())
        if self.get_particle_type(real_dia.get_initial_id(), get_code=True) \
                != ab_dia.get_initial_id():
            return False

        # Final pids
        real_pdgcode_list = [l['id'] for l in real_dia.get_final_legs()]

        # Full comparision of particle type.
        # Duplicated particles will be treated as different.
        try:
            # Try to use final_legs_dict to save time
            temp = [ab2realdict['final_legs_dict'][l.get('id')] for l in ab_dia.get_final_legs()]
            assert_pidlist = []
            # Compare the full list, temp may contain sublist,
            # if one abstract pid corresponds to several real pids.
            for item in temp:
                if isinstance(item, list):
                    assert_pidlist.extend(item)
                else:
                    assert_pidlist.append(item)

            if assert_pidlist != real_pdgcode_list:
                return False

        except:
            if sorted(self.get_particlelist_type(real_pdgcode_list,
                                                 sn_dict = {ini_type: 1})[0])  != sorted(ab_dia['abstract_type'][2]):
                return False

            
        return True


    # Helper function
    def add_ab_diagram(self, ab_amp, real_dia):
        """ Add abstract diagram from real one into the abstract amplitude.
        The abstract_type of abstract diagram is constructed. """

        if not isinstance(ab_amp, DecayAmplitude):
            raise self.PhysicsObjectError,\
                "The first argument  is not DecayAmplitude object."
        if not isinstance(real_dia, Channel):
            raise self.PhysicsObjectError,\
                "The second argument is not Channel object."

        ab_dia = Channel({'onshell': True})
        ab_dia['vertices'] = copy.deepcopy(real_dia['vertices'])

        # Setup the interaction ids
        real_inter_id_list = [v.get('id') for v in real_dia['vertices']]

        # Setup the abstract interaction id type first
        new_inter_ids = self.get_interactionlist_type(real_inter_id_list)[0]
        ab_dia['abstract_type'][0] = new_inter_ids

        # Setup the abstract interaction ids from inter_sn_dict in ab_amp,
        # then recycle the sn record
        new_inter_ids, ab_amp['inter_sn_dict'] = \
            self.get_interactionlist_type(real_inter_id_list,
                                          sn_dict = ab_amp['inter_sn_dict'])
        for i, ab_inter_id in enumerate(new_inter_ids):
            ab_dia['vertices'][i]['id'] = ab_inter_id


        # Setup the initial legs, which pid should be the pid of its
        # antiparticle.
        # Also setup the 'ini_pid' property
        ini_type = self.get_particle_type(real_dia.get_initial_id())
        ini_code = self.get_particle_type(real_dia.get_initial_id(), 
                                          get_code=True)
        ab_dia['ini_pid'] = ini_code
        if ini_type[2]:
            ab_dia['vertices'][-1]['legs'][-1]['id'] = ini_code
        else:
            ab_dia['vertices'][-1]['legs'][-1]['id'] = -ini_code


        # Setup the final abstract particle id
        real_pdgcode_list = [l.get('id') for l in real_dia.get_final_legs()]
        ab_pid_list = self.get_particlelist_type(real_pdgcode_list, 
                                                 sn_dict={ini_type:1})[0]
        # The abstract diagram is copied from real diagram,
        # so the ab_pid_list, given according to the final_legs in
        # real diagram, can be used directly to abstract diagram.
        for i, leg in enumerate(ab_dia.get_final_legs()):
            leg['id'] = ab_pid_list[i]

        # Update the abstract_type
        ab_dia['abstract_type'][2] = ab_pid_list



        # Setup the intermediate Particle id list
        real_pdgcode_list = [v.get('legs')[-1]['id'] \
                                 for v in real_dia['vertices'][:-1]]

        # Setup the abstract particle id type first
        new_part_ids = self.get_particlelist_type(real_pdgcode_list)[0]
        ab_dia['abstract_type'][1] = new_part_ids

        # Setup the abstract interaction ids from the part_sn_dict of ab_amp,
        # then recycle the sn_dict.
        new_part_ids, ab_amp['part_sn_dict'] = \
            self.get_particlelist_type(real_pdgcode_list,
                                       sn_dict=ab_amp['part_sn_dict'])
        for i, ab_pid in enumerate(new_part_ids):
            ab_dia['vertices'][i]['legs'][-1]['id'] = ab_pid

            # Find the previous leg that connect with this one
            for v in ab_dia['vertices'][i+1:]:
                # Use the number to identify the same leg
                for l in v.get('legs')[:-1]:
                    if l['number'] == \
                            ab_dia['vertices'][i]['legs'][-1]['number']:
                        l['id'] = ab_pid
                        break                


        # Add this diagram into the amplitude.
        ab_amp.add_std_diagram(ab_dia, self)

            

    def generate_ab_amplitudes(self, amp_list):
        """ Generate the abstract Amplitudes from real amplitudes of the
        SAME initial particle, 
        then generating the AbstractMatrixElement. """

        # Skip empty list
        if not amp_list:
            return

        # Get the abstract initial id        
        ini_pdg = amp_list[0]['process'].get_initial_ids()[0]
        ab_ini_pdg = self.get_particle_type(ini_pdg, get_code=True)
        ab_ini = self.get_particle(ab_ini_pdg)

        for amp in amp_list:
            # Check if this abstract amplitude exists
            final_ids = [l.get('id') for l in amp['process'].get_final_legs()]
            # PDG code start from the initial one
            pseudo_pids, ini_sn_dict = \
                self.get_particlelist_type(final_ids,\
                                               sn_dict = {self.get_particle_type(ini_pdg): 1})
            ab_amp = ab_ini.get_amplitude(pseudo_pids)

            # Construct abstract Amplitude if not exist
            if not ab_amp:
                # Construct the process
                new_process = base_objects.Process({'model': self})
                new_process['legs'] = copy.deepcopy(amp['process']['legs'])
                i = 0
                for new_leg in new_process['legs']:
                    # Final legs
                    if new_leg.get('state'):
                        new_leg['id'] = pseudo_pids[i]
                        i += 1
                    else:
                        new_leg['id'] = ab_ini_pdg

                # Construct Amplitude
                ab_amp = DecayAmplitude({'part_sn_dict': ini_sn_dict})
                ab_amp['process'] = new_process
                try:
                    ab_ini.decay_amplitudes[len(final_ids)].append(ab_amp)
                except KeyError:
                    ab_ini.decay_amplitudes[len(final_ids)] = DecayAmplitudeList([ab_amp])

            # Create the Ab2RealDict for this real amplitude
            # Set the final_legs_dict, using the pseudo_pids.
            ab_amp['ab2real_dicts'].append(Ab2RealDict({'process':
                                                            amp['process'],
                                                        'ab_process':
                                                            ab_amp['process']}))
            ab_amp['ab2real_dicts'][-1].set_final_legs_dict()

            # Scanning the diagrams in real amplitude,
            # Create new diagrams if necessary
            for i, dia in enumerate(amp['diagrams']):
                not_exist = True
                # Check diagrams in abstract amplitude
                for j, ab_dia in enumerate(ab_amp['diagrams']):

                    if not j in ab_amp['ab2real_dicts'][-1]['dia_sn_dict'].keys() and self.compare_diagrams(ab_dia, dia, ab_amp['ab2real_dicts'][-1]):
                        
                        # Update the dia_sn_dict
                        ab_amp['ab2real_dicts'][-1]['dia_sn_dict'][j] = i
                        not_exist = False
                        break

                # Create new diagram if necessary
                if not_exist:
                    # the new abstract diagram corresponds to this real diagram
                    ab_amp['ab2real_dicts'][-1]['dia_sn_dict']\
                        [len(ab_amp['diagrams'])] = i
                    self.add_ab_diagram(ab_amp, dia)

                        
            # Construct the variable dicts
            ab_amp.generate_variables_dicts(amp)


    def generate_ab_matrixelements(self, ab_amplist):
        """ Generate abstract matrix elements for the given
        abstract amplitude list. """

        ab_matrix_elements = AbstractHelasMatrixElementList()

        for ab_amp in ab_amplist:
            ab_matrix_elements.append(AbstractHelasMatrixElement(ab_amp))

        return ab_matrix_elements


    def generate_ab_matrixelements_all(self):
        """ Generate all abstract matrix elements in this model and
        save inside ab_matrix_elements. """

        ab_matrix_elements = AbstractHelasMatrixElementList()

        # Generate all matrix elements in this model
        for part in self['particles']:            
            for clevel, amps in part.decay_amplitudes.items():
                ab_matrix_elements.extend(self.generate_ab_matrixelements(amps))

        self['ab_matrix_elements'] = ab_matrix_elements

        return ab_matrix_elements

#===============================================================================
# Ab2RealDict
#===============================================================================
class Ab2RealDict(base_objects.PhysicsObject):
    """A Reference dict to mapping the information of an abstract amplitude
    to a real amplitude."""


    sorted_keys = ['process', 'dia_sn_dict', 'mass_dict', 'coup_dict',
                   'final_legs_dict']

    def default_setup(self):
        """Default values for all properties."""

        # Dictionary of from serial numbers of diagrams in abstract amplitude
        # to real amplitude
        self['process'] = base_objects.Process()
        self['ab_process'] = base_objects.Process()
        self['dia_sn_dict'] = {}
        self['mass_dict'] = {}
        self['coup_dict'] = {}
        self['final_legs_dict'] = {}

    def set_final_legs_dict(self, ab_object=None, real_dia=None):
        """ Setup the final_legs_dict."""

        # Set real pids
        if real_dia:
            real_pids = [l.get('id') for l in real_dia.get_final_legs()]
            ini_pid = real_dia.get_initial_id()
        else:
            real_pids = [l.get('id') for l in self['process'].get_final_legs()]
            ini_pid = self['process'].get_initial_ids()[0]

        # Set ab pids
        if not ab_object:
            ab_model = self['ab_process']['model']
            ab_pids = ab_model.get_particlelist_type(real_pids, sn_dict = {ab_model.get_particle_type(ini_pid): 1})[0]
        elif isinstance(ab_object, AbstractModel):
            ab_pids = ab_object.get_particlelist_type(real_pids, sn_dict = {ab_object.get_particle_type(ini_pid): 1})[0]
        else:
            raise self.PhysicsObjectError, \
                "The first argument is not necessary, otherwise it should be AbstractModel."

        # Reset final_legs_dict
        self['final_legs_dict'] = {}

        # So far, each ab_pid only has one real_pid
        for k, ab_pid in enumerate(ab_pids):
            if not ab_pid in self['final_legs_dict'].keys():
                self['final_legs_dict'][ab_pid] = real_pids[k]
            elif isinstance(self['final_legs_dict'][ab_pid], list):
                logger.warning('multiple real id correspondence')
                self['final_legs_dict'][ab_pid].append(real_pids[k])
            else:
                logger.warning('multiple real id correspondence')
                self['final_legs_dict'][ab_pid] = \
                    [self['final_legs_dict'][ab_pid], real_pids[k]]

#===============================================================================
# Ab2RealDictList
#===============================================================================
class Ab2RealDictList(base_objects.PhysicsObjectList):
    """A InteractionList, with additional properties that stores the
       abstract interaction type."""

    def is_valid_element(self, obj):
        """ Test if the object is a valid Ab2RealDictList for the list. """
        return isinstance(obj, Ab2RealDict)


#===============================================================================
# AbstractHelasMatrixElement
#===============================================================================
class AbstractHelasMatrixElement(helas_objects.HelasMatrixElement):
    """A HelasMatrixElement that contains the Ab2RealDict."""


    def default_setup(self):
        """Default values for all properties."""

        super(AbstractHelasMatrixElement, self).default_setup()
        self['ab2real_dicts'] = Ab2RealDictList()


    # Customized constructor
    def __init__(self, amplitude=None, optimization=1,
                 decay_ids=[], gen_color=True):
        """Constructor for the AbstractHelasMatrixElement.
        If DecayAmplitude is provided, copy the Ab2RealDictList.
        """

        super(AbstractHelasMatrixElement, self).__init__()

        if isinstance(amplitude, DecayAmplitude):
            self['ab2real_dicts'] = amplitude['ab2real_dicts']


#===============================================================================
# AbstractHelasMatrixElementList
#===============================================================================
class AbstractHelasMatrixElementList(helas_objects.HelasMatrixElementList):
    """List of AbstractHelasMatrixElement objects
    """

    def is_valid_element(self, obj):
        """Test if object obj is a valid HelasMatrixElement for the list."""

        return isinstance(obj, AbstractHelasMatrixElement)

    
#===============================================================================
# Helper function
#===============================================================================
def list_abs(list):

    return [abs(item) for item in list]

def legcmp(x, y):
    """Define the leg comparison, useful when testEqual is execute"""
    mycmp = cmp(x['id'], y['id'])
    if mycmp == 0:
        mycmp = cmp(x['state'], y['state'])
    return mycmp

def legcmp_bynumber(x, y):
    """Define the leg comparison, useful in generation of process in DecayAmplitude."""
    mycmp = cmp(x['number'], y['number'])
    return mycmp

def id_num_cmp(x, y):
    """Define the leg (id, number) sort."""
    mycmp = cmp(x[0], y[0])
    if mycmp == 0:
        mycmp = cmp(x[1], y[1])
    return mycmp

def channelcmp_width(x, y):
    """ Sort the channels by their width."""
    if x['onshell']:
        mycmp = cmp(x['apx_decaywidth'], y['apx_decaywidth'])
    else:
        mycmp = cmp(x['apx_decaywidth_nextlevel'], y['apx_decaywidth_nextlevel'])
    return -mycmp

def channelcmp_final(x, y):
    """ Sort the channels by their final_mass_list. 
        This will be similar to sort by the final state particles."""

    mycmp = cmp(x['final_mass_list'], y['final_mass_list'])

    return -mycmp

def amplitudecmp_width(x, y):
    """ Sort the amplitudes by their width."""
    mycmp = cmp(x['apx_decaywidth'].real, y['apx_decaywidth'].real)

    return -mycmp

def part_type_cmp(x, y):
    """ Sort the abstract particle type."""
    mycmp = cmp(x[0], y[0])

    if mycmp == 0:
        mycmp = cmp(x[1], y[1])

    return mycmp

def part_cmp(x, y):
    """ Sort the particle according to signed pdg_code."""

    return cmp(x.get_pdg_code(), y.get_pdg_code())