~maddevelopers/mg5amcnlo/3.0.2_loopparams

################################################################################
#
# Copyright (c) 2009 The MadGraph5_aMC@NLO Development team and Contributors
#
# This file is a part of the MadGraph5_aMC@NLO 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 MadGraph5_aMC@NLO license which should accompany this 
# distribution.
#
# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
#
################################################################################
"""Definitions of all basic objects used in the core code: particle, 
interaction, model, leg, vertex, process, ..."""

import copy
import itertools
import logging
import math
import numbers
import os
import re
import StringIO
import madgraph.core.color_algebra as color
from madgraph import MadGraph5Error, MG5DIR, InvalidCmd
import madgraph.various.misc as misc 


logger = logging.getLogger('madgraph.base_objects')
pjoin = os.path.join

#===============================================================================
# PhysicsObject
#===============================================================================
class PhysicsObject(dict):
    """A parent class for all physics objects."""

    class PhysicsObjectError(Exception):
        """Exception raised if an error occurs in the definition
        or the execution of a physics object."""
        pass

    def __init__(self, init_dict={}):
        """Creates a new particle object. If a dictionary is given, tries to 
        use it to give values to properties."""

        dict.__init__(self)
        self.default_setup()

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


        for item in init_dict.keys():
            self.set(item, init_dict[item])
        

    def __getitem__(self, name):
        """ force the check that the property exist before returning the 
            value associated to value. This ensure that the correct error 
            is always raise
        """

        try:
            return dict.__getitem__(self, name)
        except KeyError:
            self.is_valid_prop(name) #raise the correct error


    def default_setup(self):
        """Function called to create and setup default values for all object
        properties"""
        pass

    def is_valid_prop(self, name):
        """Check if a given property name is valid"""

        assert isinstance(name, str), \
                                 "Property name %s is not a string" % repr(name)

        if name not in self.keys():
            raise self.PhysicsObjectError, \
                        """%s is not a valid property for this object: %s\n
    Valid property are %s""" % (name,self.__class__.__name__, self.keys())
        return True

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

        return self[name]

    def set(self, name, value, force=False):
        """Set the value of the property name. First check if value
        is a valid value for the considered property. Return True if the
        value has been correctly set, False otherwise."""
        if not __debug__ or force:
            self[name] = value
            return True

        if self.is_valid_prop(name):
            try:
                self.filter(name, value)
                self[name] = value
                return True
            except self.PhysicsObjectError, why:
                logger.warning("Property " + name + " cannot be changed:" + \
                                str(why))
                return False

    def filter(self, name, value):
        """Checks if the proposed value is valid for a given property
        name. Returns True if OK. Raises an error otherwise."""

        return True

    def get_sorted_keys(self):
        """Returns the object keys sorted in a certain way. By default,
        alphabetical."""

        return self.keys().sort()

    def __str__(self):
        """String representation of the object. Outputs valid Python 
        with improved format."""

        mystr = '{\n'
        for prop in self.get_sorted_keys():
            if isinstance(self[prop], str):
                mystr = mystr + '    \'' + prop + '\': \'' + \
                        self[prop] + '\',\n'
            elif isinstance(self[prop], float):
                mystr = mystr + '    \'' + prop + '\': %.2f,\n' % self[prop]
            else:
                mystr = mystr + '    \'' + prop + '\': ' + \
                        repr(self[prop]) + ',\n'
        mystr = mystr.rstrip(',\n')
        mystr = mystr + '\n}'

        return mystr

    __repr__ = __str__


#===============================================================================
# PhysicsObjectList
#===============================================================================
class PhysicsObjectList(list):
    """A class to store lists of physics object."""

    class PhysicsObjectListError(Exception):
        """Exception raised if an error occurs in the definition
        or execution of a physics object list."""
        pass

    def __init__(self, init_list=None):
        """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)
                
    def append(self, object):
        """Appends an element, but test if valid before."""
        
        assert self.is_valid_element(object), \
            "Object %s is not a valid object for the current list" % repr(object)

        list.append(self, object)
        

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

    def __str__(self):
        """String representation of the physics object list object. 
        Outputs valid Python with improved format."""

        mystr = '['

        for obj in self:
            mystr = mystr + str(obj) + ',\n'

        mystr = mystr.rstrip(',\n')

        return mystr + ']'

#===============================================================================
# Particle
#===============================================================================
class Particle(PhysicsObject):
    """The particle object containing the whole set of information required to
    univocally characterize a given type of physical particle: name, spin, 
    color, mass, width, charge,... The is_part flag tells if the considered
    particle object is a particle or an antiparticle. The self_antipart flag
    tells if the particle is its own antiparticle."""

    sorted_keys = ['name', 'antiname', 'spin', 'color',
                   'charge', 'mass', 'width', 'pdg_code',
                  'line', 'propagator',
                   'is_part', 'self_antipart', 'type', 'counterterm']

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

        self['name'] = 'none'
        self['antiname'] = 'none'
        self['spin'] = 1
        self['color'] = 1
        self['charge'] = 1.
        self['mass'] = 'ZERO'
        self['width'] = 'ZERO'
        self['pdg_code'] = 0
        #self['texname'] = 'none'
        #self['antitexname'] = 'none'
        self['line'] = 'dashed'
        #self['propagating'] = True -> removed in favor or 'line' = None
        self['propagator'] = ''
        self['is_part'] = True
        self['self_antipart'] = False
        # True if ghost, False otherwise
        #self['ghost'] = False
        self['type'] = '' # empty means normal can also be ghost or goldstone 
        # Counterterm defined as a dictionary with format:
        # ('ORDER_OF_COUNTERTERM',((Particle_list_PDG))):{laurent_order:CTCouplingName}
        self['counterterm'] = {}

    def get(self, name):
        
        if name == 'ghost':
            return self['type'] == 'ghost'
        elif name == 'goldstone':
            return self['type'] == 'goldstone'
        elif name == 'propagating':
            return self['line'] not in ['None',None] 
        else:
            return super(Particle, self).get(name)

    def set(self, name, value, force=False):
        
        if name in ['texname', 'antitexname']:
            return True
        elif name == 'propagating':
            if not value:
                return self.set('line', None, force=force)
            elif not self.get('line'):
                return self.set('line', 'dashed',force=force)
            return True
        elif name  in ['ghost', 'goldstone']:
            if self.get('type') == name:
                if value:
                    return True
                else:
                    return self.set('type', '', force=force)
            else:
                if value:
                    return self.set('type', name, force=force)
                else:
                    return True
        return super(Particle, self).set(name, value,force=force)
        

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

        if name in ['name', 'antiname']:
            # Forbid special character but +-~_
            p=re.compile('''^[\w\-\+~_]+$''')
            if not p.match(value):
                raise self.PhysicsObjectError, \
                        "%s is not a valid particle name" % value

        if name is 'ghost':
            if not isinstance(value,bool):
                raise self.PhysicsObjectError, \
                 "%s is not a valid bool for the 'ghost' attribute" % str(value)
                 
        if name is 'counterterm':
            if not isinstance(value,dict):
                raise self.PhysicsObjectError, \
                    "counterterm %s is not a valid dictionary" % repr(value)
            for key, val in value.items():
                if not isinstance(key,tuple):
                    raise self.PhysicsObjectError, \
                        "key %s is not a valid tuple for counterterm key" % repr(key)
                if not isinstance(key[0],str):
                    raise self.PhysicsObjectError, \
                        "%s is not a valid string" % repr(key[0])
                if not isinstance(key[1],tuple):
                    raise self.PhysicsObjectError, \
                        "%s is not a valid list" % repr(key[1])
                for elem in key[1]:
                    if not isinstance(elem,tuple):
                        raise self.PhysicsObjectError, \
                            "%s is not a valid list" % repr(elem)
                    for partPDG in elem:
                        if not isinstance(partPDG,int):
                            raise self.PhysicsObjectError, \
                                "%s is not a valid integer for PDG" % repr(partPDG)
                        if partPDG<=0:
                            raise self.PhysicsObjectError, \
                                "%s is not a valid positive PDG" % repr(partPDG)
                if not isinstance(val,dict):
                    raise self.PhysicsObjectError, \
                        "value %s is not a valid dictionary for counterterm value" % repr(val)
                for vkey, vvalue in val.items():
                    if vkey not in [0,-1,-2]:
                        raise self.PhysicsObjectError, \
                            "Key %s is not a valid laurent serie order" % repr(vkey)
                    if not isinstance(vvalue,str):
                        raise self.PhysicsObjectError, \
                            "Coupling %s is not a valid string" % repr(vvalue)
        if name is 'spin':
            if not isinstance(value, int):
                raise self.PhysicsObjectError, \
                    "Spin %s is not an integer" % repr(value)
            if (value < 1 or value > 5) and value != 99:
                raise self.PhysicsObjectError, \
                   "Spin %i not valid" % value

        if name is 'color':
            if not isinstance(value, int):
                raise self.PhysicsObjectError, \
                    "Color %s is not an integer" % repr(value)
            if value not in [1, 3, 6, 8]:
                raise self.PhysicsObjectError, \
                   "Color %i is not valid" % value

        if name in ['mass', 'width']:
            # Must start with a letter, followed by letters, digits or _
            p = re.compile('\A[a-zA-Z]+[\w\_]*\Z')
            if not p.match(value):
                raise self.PhysicsObjectError, \
                        "%s is not a valid name for mass/width variable" % \
                        value

        if name is 'pdg_code':
            if not isinstance(value, int):
                raise self.PhysicsObjectError, \
                    "PDG code %s is not an integer" % repr(value)

        if name is 'line':
            if not isinstance(value, str):
                raise self.PhysicsObjectError, \
                    "Line type %s is not a string" % repr(value)
            if value not in ['None','dashed', 'straight', 'wavy', 'curly', 'double','swavy','scurly','dotted']:
                raise self.PhysicsObjectError, \
                   "Line type %s is unknown" % value

        if name is 'charge':
            if not isinstance(value, float):
                raise self.PhysicsObjectError, \
                    "Charge %s is not a float" % repr(value)

        if name is 'propagating':
            if not isinstance(value, bool):
                raise self.PhysicsObjectError, \
                    "Propagating tag %s is not a boolean" % repr(value)

        if name in ['is_part', 'self_antipart']:
            if not isinstance(value, bool):
                raise self.PhysicsObjectError, \
                    "%s tag %s is not a boolean" % (name, repr(value))

        return True

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

        return self.sorted_keys

    # Helper functions

    def is_perturbating(self,order,model):
        """Returns wether this particle contributes in perturbation of the order passed
           in argument given the model specified. It is very fast for usual models"""
           
        for int in model['interactions'].get_type('base'):
            # We discard the interactions with more than one type of orders
            # contributing because it then doesn't necessarly mean that this
            # particle (self) is charged under the group corresponding to the
            # coupling order 'order'. The typical example is in SUSY which 
            # features a ' photon-gluon-squark-antisquark ' interaction which
            # has coupling orders QED=1, QCD=1 and would induce the photon
            # to be considered as a valid particle to circulate in a loop of
            # type "QCD".
            if len(int.get('orders'))>1:
                continue
            if order in int.get('orders').keys() and self.get('pdg_code') in \
              [part.get('pdg_code') for part in int.get('particles')]:
                return True
            
        return False
           
    def get_pdg_code(self):
        """Return the PDG code with a correct minus sign if the particle is its
        own antiparticle"""

        if not self['is_part'] and not self['self_antipart']:
            return - self['pdg_code']
        else:
            return self['pdg_code']

    def get_anti_pdg_code(self):
        """Return the PDG code of the antiparticle with a correct minus sign 
        if the particle is its own antiparticle"""

        if not self['self_antipart']:
            return - self.get_pdg_code()
        else:
            return self['pdg_code']

    def get_color(self):
        """Return the color code with a correct minus sign"""

        if not self['is_part'] and abs(self['color']) in [3, 6]:
            return - self['color']
        else:
            return self['color']

    def get_anti_color(self):
        """Return the color code of the antiparticle with a correct minus sign
        """

        if self['is_part'] and self['color'] not in [1, 8]:
            return - self['color']
        else:
            return self['color']
        
    def get_charge(self):
        """Return the charge code with a correct minus sign"""

        if not self['is_part']:
            return - self['charge']
        else:
            return self['charge']

    def get_anti_charge(self):
        """Return the charge code of the antiparticle with a correct minus sign
        """

        if self['is_part']:
            return - self['charge']
        else:
            return self['charge']
        
    def get_name(self):
        """Return the name if particle, antiname if antiparticle"""

        if not self['is_part'] and not self['self_antipart']:
            return self['antiname']
        else:
            return self['name']

    def get_helicity_states(self, allow_reverse=True):
        """Return a list of the helicity states for the onshell particle"""

        spin = self.get('spin')
        if spin ==1:
            # Scalar
            res = [ 0 ]
        elif spin == 2:
            # Spinor
            res = [ -1, 1 ]                
        elif spin == 3 and self.get('mass').lower() == 'zero':
            # Massless vector
            res = [ -1, 1 ]
        elif spin == 3:
            # Massive vector
            res = [ -1, 0, 1 ]
        elif spin == 4 and self.get('mass').lower() == 'zero':
            # Massless tensor
            res = [-3, 3]
        elif spin == 4:
            # Massive tensor
            res = [-3, -1, 1, 3]
        elif spin == 5 and self.get('mass').lower() == 'zero':
            # Massless tensor
            res = [-2, -1, 1, 2]
        elif spin in [5, 99]:
            # Massive tensor
            res = [-2, -1, 0, 1, 2]
        else:
            raise self.PhysicsObjectError, \
              "No helicity state assignment for spin %d particles" % spin
                  
        if allow_reverse and not self.get('is_part'):
            res.reverse()


        return res

    def is_fermion(self):
        """Returns True if this is a fermion, False if boson"""

        return self['spin'] % 2 == 0

    def is_boson(self):
        """Returns True if this is a boson, False if fermion"""

        return self['spin'] % 2 == 1

#===============================================================================
# ParticleList
#===============================================================================
class ParticleList(PhysicsObjectList):
    """A class to store lists of particles."""

    def is_valid_element(self, obj):
        """Test if object obj is a valid Particle for the list."""
        return isinstance(obj, Particle)
                    
    def get_copy(self, name):
        """Try to find a particle with the given name. Check both name
        and antiname. If a match is found, return the a copy of the 
        corresponding particle (first one in the list), with the 
        is_part flag set accordingly. None otherwise."""
        
        assert isinstance(name, str)
        
        part = self.find_name(name)
        if not part:
            # Then try to look for a particle with that PDG
            try:
                pdg = int(name)
            except ValueError:
                return None

            for p in self:
                if p.get_pdg_code()==pdg:
                    part = copy.copy(p)
                    part.set('is_part', True)
                    return part
                elif p.get_anti_pdg_code()==pdg:
                    part = copy.copy(p)
                    part.set('is_part', False)
                    return part

            return None
        part = copy.copy(part)
          
        if part.get('name') == name:
            part.set('is_part', True)
            return part
        elif part.get('antiname') == name:
            part.set('is_part', False)
            return part
        return None

    def find_name(self, name):
        """Try to find a particle with the given name. Check both name
        and antiname. If a match is found, return the a copy of the 
        corresponding particle (first one in the list), with the 
        is_part flag set accordingly. None otherwise."""

        assert isinstance(name, str), "%s is not a valid string" % str(name) 

        for part in self:
            if part.get('name') == name:
                return part
            elif part.get('antiname') == name:
                return part

        return None

    def generate_ref_dict(self):
        """Generate a dictionary of part/antipart pairs (as keys) and
        0 (as value)"""

        ref_dict_to0 = {}

        for part in self:
            ref_dict_to0[(part.get_pdg_code(), part.get_anti_pdg_code())] = [0]
            ref_dict_to0[(part.get_anti_pdg_code(), part.get_pdg_code())] = [0]

        return ref_dict_to0

    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.deepcopy(particle)
                antipart.set('is_part', False)
                particle_dict[antipart.get_pdg_code()] = antipart

        return particle_dict


#===============================================================================
# Interaction
#===============================================================================
class Interaction(PhysicsObject):
    """The interaction object containing the whole set of information 
    required to univocally characterize a given type of physical interaction: 
    
    particles: a list of particle ids
    color: a list of string describing all the color structures involved
    lorentz: a list of variable names describing all the Lorentz structure
             involved
    couplings: dictionary listing coupling variable names. The key is a
               2-tuple of integers referring to color and Lorentz structures
    orders: dictionary listing order names (as keys) with their value
    """

    sorted_keys = ['id', 'particles', 'color', 'lorentz', 'couplings',
                   'orders','loop_particles','type','perturbation_type']

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

        self['id'] = 0
        self['particles'] = []
        self['color'] = []
        self['lorentz'] = []
        self['couplings'] = { (0, 0):'none'}
        self['orders'] = {}
        # The type of interactions can be 'base', 'UV' or 'R2'.
        # For 'UV' or 'R2', one can always specify the loop it corresponds
        # to by a tag in the second element of the list. If the tag is an
        # empty list, then the R2/UV interaction will be recognized only
        # based on the nature of the identity of the particles branching
        # off the loop and the loop orders. 
        # Otherwise, the tag can be specified and it will be used when 
        # identifying the R2/UV interaction corresponding to a given loop
        # generated.
        # The format is [(lp1ID,int1ID),(lp1ID,int1ID),(lp1ID,int1ID),etc...]
        # Example of a tag for the following loop
        #
        #             ___34_____   The ';' line is a gluon with ID 21
        #          45/   ;         The '|' line is a d-quark with ID 1
        #     ------<    ;         The numbers are the interactions ID
        #            \___;______   The tag for this loop would be:
        #                12          ((21,34),(1,45),(1,12))
        #                         
        # This tag is equivalent to all its cyclic permutations. This is why
        # it must be specified in the canonical order which is defined with 
        # by putting in front of the tag the lowest 2-tuple it contains.
        # (the order relation is defined by comparing the particle ID first
        # and the interaction ID after in case the particle ID are the same).
        # In case there are two identical lowest 2-tuple in the tag, the
        # tag chosen is such that it has the lowest second 2-tuple. The procedure
        # is repeated again with the subsequent 2-tuple until there is only
        # one cyclic permutation remaining and the ambiguity is resolved.
        # This insures to have one unique unambiguous canonical tag chosen.
        # In the example above, it would be:
        #       ((1,12),(21,34),(1,45))
        # PS: Notice that in the UFO model, the tag-information is limited to 
        # the minimally relevant one which are the loop particles specified in
        # in the attribute below. In this case, 'loop_particles' is the list of 
        # all the loops giving this same counterterm contribution. 
        # Each loop being represented by a set of the PDG of the particles 
        # (not repeated) constituting it. In the example above, it would simply
        # be (1,21). In the UFO, if the loop particles are not specified then
        # MG5 will account for this counterterm only once per concerned vertex.
        # Taking the example of the three gluon vertex counterterm, one can
        # possibly have in the ufo:
        #                VertexB = blabla, loop_particles = (b)
        #                VertexT = blabla, loop_particles = (t)
        # or 
        #                VertexALL = blabla, loop_particles = ()
        # In the first case UFO specifies the specific counterterm to the three-
        # gluon loop with the bottom running in (VertexB) and with the top running
        # in (VertexT). So MG5 will associate these counterterm vertices once to
        # each of the two loop.
        # In the case where UFO defined VertexALL, then whenever MG5 encounters
        # a triangle three-gluon loop (say the bottom one), it will associate to
        # it the vertex VertexALL but will not do so again when encountering the 
        # same loop with the top quark running in. This, because it assumes that
        # the UFO vertexALL comprises all contributions already.
        
        self['loop_particles']=[[]]
        self['type'] = 'base'
        self['perturbation_type'] = None

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

        if name == 'id':
            #Should be an integer
            if not isinstance(value, int):
                raise self.PhysicsObjectError, \
                        "%s is not a valid integer" % str(value)

        if name == 'particles':
            #Should be a list of valid particle names
            if not isinstance(value, ParticleList):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list of particles" % str(value)

        if name == 'perturbation_type':
            if value!=None and not isinstance(value, str):
                raise self.PhysicsObjectError, \
                        "%s is not a valid string" % str(value)            

        if name == 'type':
            #Should be a string
            if not isinstance(value, str):
                raise self.PhysicsObjectError, \
                        "%s is not a valid string" % str(value)
        if name == 'loop_particles':
            if isinstance(value,list):
                for l in value:
                    if isinstance(l,list):
                        for part in l:
                            if not isinstance(part,int):
                                raise self.PhysicsObjectError, \
                                    "%s is not a valid integer" % str(part)
                            if part<0:
                                raise self.PhysicsObjectError, \
                                    "%s is not a valid positive integer" % str(part)

        if name == 'orders':
            #Should be a dict with valid order names ask keys and int as values
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                        "%s is not a valid dict for coupling orders" % \
                                                                    str(value)
            for order in value.keys():
                if not isinstance(order, str):
                    raise self.PhysicsObjectError, \
                        "%s is not a valid string" % str(order)
                if not isinstance(value[order], int):
                    raise self.PhysicsObjectError, \
                        "%s is not a valid integer" % str(value[order])

        if name in ['color']:
            #Should be a list of list strings
            if not isinstance(value, list):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list of Color Strings" % str(value)
            for mycolstring in value:
                if not isinstance(mycolstring, color.ColorString):
                    raise self.PhysicsObjectError, \
                            "%s is not a valid list of Color Strings" % str(value)

        if name in ['lorentz']:
            #Should be a list of list strings
            if not isinstance(value, list):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list of strings" % str(value)
            for mystr in value:
                if not isinstance(mystr, str):
                    raise self.PhysicsObjectError, \
                        "%s is not a valid string" % str(mystr)

        if name == 'couplings':
            #Should be a dictionary of strings with (i,j) keys
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                        "%s is not a valid dictionary for couplings" % \
                                                                str(value)

            for key in value.keys():
                if not isinstance(key, tuple):
                    raise self.PhysicsObjectError, \
                        "%s is not a valid tuple" % str(key)
                if len(key) != 2:
                    raise self.PhysicsObjectError, \
                        "%s is not a valid tuple with 2 elements" % str(key)
                if not isinstance(key[0], int) or not isinstance(key[1], int):
                    raise self.PhysicsObjectError, \
                        "%s is not a valid tuple of integer" % str(key)
                if not isinstance(value[key], str):
                    raise self.PhysicsObjectError, \
                        "%s is not a valid string" % value[key]

        return True

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

        return self.sorted_keys 
                
    def is_perturbating(self, orders_considered):
        """ Returns if this interaction comes from the perturbation of one of
        the order listed in the argument """
        
        if self['perturbation_type']==None:
            return True
        else:
            return (self['perturbation_type'] in orders_considered)
                
    def is_R2(self):
        """ Returns if the interaction is of R2 type."""

        # Precaution only useful because some tests have a predefined model
        # bypassing the default_setup and for which type was not defined.
        if 'type' in self.keys():
            return (len(self['type'])>=2 and self['type'][:2]=='R2')
        else:
            return False

    def is_UV(self):
        """ Returns if the interaction is of UV type."""

        # Precaution only useful because some tests have a predefined model
        # bypassing the default_setup and for which type was not defined.
        if 'type' in self.keys():
            return (len(self['type'])>=2 and self['type'][:2]=='UV')
        else:
            return False
        
    def is_UVmass(self):
        """ Returns if the interaction is of UVmass type."""

        # Precaution only useful because some tests have a predefined model
        # bypassing the default_setup and for which type was not defined.
        if 'type' in self.keys():
            return (len(self['type'])>=6 and self['type'][:6]=='UVmass')
        else:
            return False
        
    def is_UVloop(self):
        """ Returns if the interaction is of UVmass type."""

        # Precaution only useful because some tests have a predefined model
        # bypassing the default_setup and for which type was not defined.
        if 'type' in self.keys():
            return (len(self['type'])>=6 and self['type'][:6]=='UVloop')
        else:
            return False
        
    def is_UVtree(self):
        """ Returns if the interaction is of UVmass type."""

        # Precaution only useful because some tests have a predefined model
        # bypassing the default_setup and for which type was not defined.
        if 'type' in self.keys():
            return (len(self['type'])>=6 and self['type'][:6]=='UVtree')
        else:
            return False
        
    def is_UVCT(self):
        """ Returns if the interaction is of the UVCT type which means that 
        it has been selected as a possible UV counterterm interaction for this
        process. Such interactions are marked by having the 'UVCT_SPECIAL' order
        key in their orders."""

        # Precaution only useful because some tests have a predefined model
        # bypassing the default_setup and for which type was not defined.
        if 'UVCT_SPECIAL' in self['orders'].keys():
            return True
        else:
            return False
        
    def get_epsilon_order(self):
        """ Returns 0 if this interaction contributes to the finite part of the
        amplitude and 1 (2) is it contributes to its single (double) pole """
        
        if 'type' in self.keys():
            if '1eps' in self['type']:
                return 1
            elif '2eps' in self['type']:
                return 2
            else:
                return 0
        else:
            return 0

    def generate_dict_entries(self, ref_dict_to0, ref_dict_to1):
        """Add entries corresponding to the current interactions to 
        the reference dictionaries (for n>0 and n-1>1)"""

        # Create n>0 entries. Format is (p1,p2,p3,...):interaction_id.
        # We are interested in the unordered list, so use sorted()

        pdg_tuple = tuple(sorted([p.get_pdg_code() for p in self['particles']]))
        if pdg_tuple not in ref_dict_to0.keys():
            ref_dict_to0[pdg_tuple] = [self['id']]
        else:
            ref_dict_to0[pdg_tuple].append(self['id'])

        # Create n-1>1 entries. Note that, in the n-1 > 1 dictionary,
        # the n-1 entries should have opposite sign as compared to
        # interaction, since the interaction has outgoing particles,
        # while in the dictionary we treat the n-1 particles as
        # incoming

        for part in self['particles']:

            # We are interested in the unordered list, so use sorted()
            pdg_tuple = tuple(sorted([p.get_pdg_code() for (i, p) in \
                                      enumerate(self['particles']) if \
                                      i != self['particles'].index(part)]))
            pdg_part = part.get_anti_pdg_code()
            if pdg_tuple in ref_dict_to1.keys():
                if (pdg_part, self['id']) not in  ref_dict_to1[pdg_tuple]:
                    ref_dict_to1[pdg_tuple].append((pdg_part, self['id']))
            else:
                ref_dict_to1[pdg_tuple] = [(pdg_part, self['id'])]

    def get_WEIGHTED_order(self, model):
        """Get the WEIGHTED order for this interaction, for equivalent
        3-particle vertex. Note that it can be fractional."""

        return float(sum([model.get('order_hierarchy')[key]*self.get('orders')[key]\
                          for key in self.get('orders')]))/ \
               max((len(self.get('particles'))-2), 1)

    def __str__(self):
        """String representation of an interaction. Outputs valid Python 
        with improved format. Overrides the PhysicsObject __str__ to only
        display PDG code of involved particles."""

        mystr = '{\n'

        for prop in self.get_sorted_keys():
            if isinstance(self[prop], str):
                mystr = mystr + '    \'' + prop + '\': \'' + \
                        self[prop] + '\',\n'
            elif isinstance(self[prop], float):
                mystr = mystr + '    \'' + prop + '\': %.2f,\n' % self[prop]
            elif isinstance(self[prop], ParticleList):
                mystr = mystr + '    \'' + prop + '\': [%s],\n' % \
                   ','.join([str(part.get_pdg_code()) for part in self[prop]])
            else:
                mystr = mystr + '    \'' + prop + '\': ' + \
                        repr(self[prop]) + ',\n'
        mystr = mystr.rstrip(',\n')
        mystr = mystr + '\n}'

        return mystr

    def canonical_repr(self):
        """ Returns a string representation that allows to order CT vertices in 
        a diagram so that identical HelasMatrix element (i.e. typically different
        flavors with identical masses and couplings) can be matched even though
        the order of the couplings specified in the UFO is different. """

        return '%s|%s'%(self['type'],
           '&'.join( sorted('%s_%s_%s'%(self['color'][k[0]],self['lorentz'][k[1]],v)
                                          for k,v in self['couplings'].items() ) ) )

    def get_canonical_couplings_keys_order(self):
        """ Returns a list of the keys to the 'couplings' dictionary, canonically 
        ordered so so that identical HelasMatrix element (i.e. typically different
        flavors with identical masses and couplings) can be matched even though
        the order of the couplings specified in the UFO is different. """

        return sorted(self['couplings'].keys(), key=lambda k:
        '%s_%s_%s'%(self['color'][k[0]],self['lorentz'][k[1]],self['couplings'][k]))


#===============================================================================
# InteractionList
#===============================================================================
class InteractionList(PhysicsObjectList):
    """A class to store lists of interactionss."""

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

        return isinstance(obj, Interaction)

    def generate_ref_dict(self,useR2UV=False, useUVCT=False):
        """Generate the reference dictionaries from interaction list.
        Return a list where the first element is the n>0 dictionary and
        the second one is n-1>1."""

        ref_dict_to0 = {}
        ref_dict_to1 = {}
        buffer = {}

        for inter in self:
            if useR2UV or (not inter.is_UV() and not inter.is_R2() and \
                           not inter.is_UVCT()):
                inter.generate_dict_entries(ref_dict_to0, ref_dict_to1)
            if useUVCT and inter.is_UVCT():
                inter.generate_dict_entries(ref_dict_to0, ref_dict_to1)
                
        return [ref_dict_to0, ref_dict_to1]

    def generate_dict(self):
        """Generate a dictionary from interaction id to interaction.
        """

        interaction_dict = {}

        for inter in self:
            interaction_dict[inter.get('id')] = inter

        return interaction_dict

    def synchronize_interactions_with_particles(self, particle_dict):
        """Make sure that the particles in the interactions are those
        in the particle_dict, and that there are no interactions
        refering to particles that don't exist. To be called when the
        particle_dict is updated in a model.
        """

        iint = 0
        while iint < len(self):
            inter = self[iint]
            particles = inter.get('particles')
            try:
                for ipart, part in enumerate(particles):
                    particles[ipart] = particle_dict[part.get_pdg_code()]
                iint += 1
            except KeyError:
                # This interaction has particles that no longer exist
                self.pop(iint)

    def get_type(self, type):
        """ return all interactions in the list of type 'type' """
        return InteractionList([int for int in self if int.get('type')==type])

    def get_R2(self):
        """ return all interactions in the list of type R2 """
        return InteractionList([int for int in self if int.is_R2()])

    def get_UV(self):
        """ return all interactions in the list of type UV """
        return InteractionList([int for int in self if int.is_UV()])

    def get_UVmass(self):
        """ return all interactions in the list of type UVmass """
        return InteractionList([int for int in self if int.is_UVmass()])

    def get_UVtree(self):
        """ return all interactions in the list of type UVtree """
        return InteractionList([int for int in self if int.is_UVtree()])
    
    def get_UVloop(self):
        """ return all interactions in the list of type UVloop """
        return InteractionList([int for int in self if int.is_UVloop()])

#===============================================================================
# Model
#===============================================================================
class Model(PhysicsObject):
    """A class to store all the model information."""
    
    mg5_name = False #store if particle name follow mg5 convention
    
    def default_setup(self):

        self['name'] = ""
        self['particles'] = ParticleList()
        self['interactions'] = InteractionList()
        self['parameters'] = None
        self['functions'] = None
        self['couplings'] = None
        self['lorentz'] = None
        self['particle_dict'] = {}
        self['interaction_dict'] = {}
        self['ref_dict_to0'] = {}
        self['ref_dict_to1'] = {}
        self['got_majoranas'] = None
        self['order_hierarchy'] = {}
        self['conserved_charge'] = set()
        self['coupling_orders'] = None
        self['expansion_order'] = None
        self['version_tag'] = None # position of the directory (for security)
        self['gauge'] = [0, 1]
        self['case_sensitive'] = True
        # attribute which might be define if needed
        #self['name2pdg'] = {'name': pdg}
        
        

    def filter(self, name, value):
        """Filter for model property values"""

        if name in ['name']:
            if not isinstance(value, str):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a string" %type(value)

        elif name == 'particles':
            if not isinstance(value, ParticleList):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a ParticleList object" % \
                                                            type(value)
        elif name == 'interactions':
            if not isinstance(value, InteractionList):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a InteractionList object" % \
                                                            type(value)
        elif name == 'particle_dict':
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a dictionary" % \
                                                        type(value)
        elif name == 'interaction_dict':
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a dictionary" % type(value)

        elif name == 'ref_dict_to0':
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a dictionary" % type(value)
                    
        elif name == 'ref_dict_to1':
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a dictionary" % type(value)

        elif name == 'got_majoranas':
            if not (isinstance(value, bool) or value == None):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a boolean" % type(value)

        elif name == 'conserved_charge':
            if not (isinstance(value, set)):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a set" % type(value)

        elif name == 'version_tag':
            if not (isinstance(value, str)):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a string" % type(value)

        elif name == 'order_hierarchy':
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a dictionary" % \
                                                            type(value)
            for key in value.keys():
                if not isinstance(value[key],int):
                    raise self.PhysicsObjectError, \
                        "Object of type %s is not an integer" % \
                                                            type(value[key])
        elif name == 'gauge':
            if not (isinstance(value, list)):
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a list" % type(value)

        elif name == 'case_sensitive':
            if not value in [True ,False]:
                raise self.PhysicsObjectError, \
                    "Object of type %s is not a boolean" % type(value)
            

        return True

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

        if (name == 'ref_dict_to0' or name == 'ref_dict_to1') and \
                                                                not self[name]:
            if self['interactions']:
                [self['ref_dict_to0'], self['ref_dict_to1']] = \
                            self['interactions'].generate_ref_dict()
                self['ref_dict_to0'].update(
                                self['particles'].generate_ref_dict())

        if (name == 'particle_dict') and not self[name]:
            if self['particles']:
                self['particle_dict'] = self['particles'].generate_dict()
            if self['interactions']:
                self['interactions'].synchronize_interactions_with_particles(\
                                                          self['particle_dict'])
        if name == 'modelpath':
            modeldir = self.get('version_tag').rsplit('##',1)[0]
            if os.path.exists(modeldir):
                modeldir = os.path.expanduser(modeldir)
                return modeldir
            else:
                raise Exception, "path %s not valid anymore." % modeldir
            #modeldir = os.path.join(os.path.dirname(modeldir),
            #                        os.path.basename(modeldir).rsplit("-",1)[0])
            #if os.path.exists(modeldir):
            #    return modeldir 
            #raise Exception, 'Invalid Path information: %s' % self.get('version_tag')          
        elif name == 'modelpath+restriction':
            modeldir = self.get('version_tag').rsplit('##',1)[0]
            modelname = self['name']            
            if not  os.path.exists(modeldir):
                raise Exception, "path %s not valid anymore" % modeldir
            modeldir = os.path.dirname(modeldir)
            modeldir = pjoin(modeldir, modelname)
            modeldir = os.path.expanduser(modeldir)
            return modeldir
        elif name == 'restrict_name':
            modeldir = self.get('version_tag').rsplit('##',1)[0]
            modelname = self['name']            
            basename = os.path.basename(modeldir)
            restriction = modelname[len(basename)+1:]
            return restriction

        if (name == 'interaction_dict') and not self[name]:
            if self['interactions']:
                self['interaction_dict'] = self['interactions'].generate_dict()

        if (name == 'got_majoranas') and self[name] == None:
            if self['particles']:
                self['got_majoranas'] = self.check_majoranas()

        if (name == 'coupling_orders') and self[name] == None:
            if self['interactions']:
                self['coupling_orders'] = self.get_coupling_orders()

        if (name == 'order_hierarchy') and not self[name]:
            if self['interactions']:
                self['order_hierarchy'] = self.get_order_hierarchy()    

        if (name == 'expansion_order') and self[name] == None:
            if self['interactions']:
                self['expansion_order'] = \
                   dict([(order, -1) for order in self.get('coupling_orders')])
                   
        if (name == 'name2pdg') and 'name2pdg' not in self:
            self['name2pdg'] = {}
            for p in self.get('particles'):
                self['name2pdg'][p.get('antiname')] = -1*p.get('pdg_code')
                self['name2pdg'][p.get('name')] =  p.get('pdg_code')
                
        return Model.__bases__[0].get(self, name) # call the mother routine

    def set(self, name, value, force = False):
        """Special set for particles and interactions - need to
        regenerate dictionaries."""

        if name == 'particles':
            # Ensure no doublets in particle list
            make_unique(value)
            # Reset dictionaries
            self['particle_dict'] = {}
            self['ref_dict_to0'] = {}
            self['got_majoranas'] = None

        if name == 'interactions':
            # Ensure no doublets in interaction list
            make_unique(value)
            # Reset dictionaries
            self['interaction_dict'] = {}
            self['ref_dict_to1'] = {}
            self['ref_dict_to0'] = {}
            self['got_majoranas'] = None
            self['coupling_orders'] = None
            self['order_hierarchy'] = {}
            self['expansion_order'] = None

        if name == 'name2pdg':
            self['name2pgg'] = value
            return

        result = Model.__bases__[0].set(self, name, value, force) # call the mother routine

        if name == 'particles':
            # Recreate particle_dict
            self.get('particle_dict')

        return result

    def actualize_dictionaries(self):
        """This function actualizes the dictionaries"""

        [self['ref_dict_to0'], self['ref_dict_to1']] = \
                self['interactions'].generate_ref_dict()
        self['ref_dict_to0'].update(
                                self['particles'].generate_ref_dict())

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

        return ['name', 'particles', 'parameters', 'interactions',
                'couplings','lorentz', 'gauge']

    def get_particle(self, id):
        """Return the particle corresponding to the id / name"""
        
        try:
            return self["particle_dict"][id]
        except Exception:
            if isinstance(id, int):
                try:
                    return self.get("particle_dict")[id]
                except Exception, error:
                    return None
            else:
                if not hasattr(self, 'name2part'):
                    self.create_name2part()
                try: 
                    return self.name2part[id]
                except:
                    return None

    def create_name2part(self):
        """create a dictionary name 2 part"""
        
        self.name2part = {}
        for part in self.get("particle_dict").values():
            self.name2part[part.get('name')] = part
            self.name2part[part.get('antiname')] = part
            
    def get_lorentz(self, name):
        """return the lorentz object from the associate name"""
        if hasattr(self, 'lorentz_name2obj'):
            return self.lorentz_name2obj[name]  
        else:
            self.create_lorentz_dict()
            return self.lorentz_name2obj[name]

    def create_lorentz_dict(self):
        """create the dictionary linked to the lorentz structure"""
        self.lorentz_name2obj = {}
        self.lorentz_expr2name = {}
        if not self.get('lorentz'):
            return
        for lor in self.get('lorentz'):
            self.lorentz_name2obj[lor.name] = lor
            self.lorentz_expr2name[lor.structure] = lor.name

    def get_interaction(self, id):
        """Return the interaction corresponding to the id"""

        try:
            return self.get("interaction_dict")[id]
        except Exception:
            return None

    def get_parameter(self, name):
        """Return the parameter associated to the name NAME"""
        
        # If information is saved
        if hasattr(self, 'parameters_dict') and self.parameters_dict:
            try:
                return self.parameters_dict[name]
            except Exception:
                # try to reload it before crashing 
                pass
            
        # Else first build the dictionary
        self.parameters_dict = {}
        for data in self['parameters'].values():
            [self.parameters_dict.__setitem__(p.name,p) for p in data]
        
        return self.parameters_dict[name]

    def get_coupling_orders(self):
        """Determine the coupling orders of the model"""
        return set(sum([i.get('orders').keys() for i in \
                        self.get('interactions')], []))

    def get_order_hierarchy(self):
        """Set a default order hierarchy for the model if not set by the UFO."""
        # Set coupling hierachy
        hierarchy = dict([(order, 1) for order in self.get('coupling_orders')])
        # Special case for only QCD and QED couplings, unless already set
        if self.get('coupling_orders') == set(['QCD', 'QED']):
            hierarchy['QED'] = 2
        return hierarchy


    def get_nflav(self):
        """returns the number of light quark flavours in the model."""
        return len([p for p in self.get('particles') \
                if p['spin'] == 2 and p['is_part'] and \
                p ['color'] != 1 and p['mass'].lower() == 'zero'])


    def get_quark_pdgs(self):
        """returns the PDG codes of the light quarks and antiquarks"""
        pdg_list = [p['pdg_code'] for p in self.get('particles') \
                       if p['spin'] == 2 and \
                       p['color'] == 3 and \
                       p['charge'] != 0. and p['mass'].lower() == 'zero']

        for p in pdg_list[:]:
            if not self.get('particle_dict')[p]['self_antipart']:
                pdg_list.append(self.get('particle_dict')[p].get_anti_pdg_code())
                
        return sorted(pdg_list)


    def get_nleps(self):
        """returns the number of light lepton flavours in the model."""
        return len([p for p in self.get('particles') \
                if p['spin'] == 2 and p['is_part'] and \
                p['color'] == 1 and \
                p['charge'] != 0. and p['mass'].lower() == 'zero'])


    def get_lepton_pdgs(self):
        """returns the PDG codes of the light leptons and antileptons"""
        pdg_list = [p['pdg_code'] for p in self.get('particles') \
                       if p['spin'] == 2 and \
                       p['color'] == 1 and \
                       p['charge'] != 0. and p['mass'].lower() == 'zero']

        for p in pdg_list[:]:
            if not self.get('particle_dict')[p]['self_antipart']:
                pdg_list.append(self.get('particle_dict')[p].get_anti_pdg_code())
                
        return sorted(pdg_list)

    
    def get_particles_hierarchy(self):
        """Returns the order hierarchies of the model and the
        particles which have interactions in at least this hierarchy
        (used in find_optimal_process_orders in MultiProcess diagram
        generation):

        Check the coupling hierarchy of the model. Assign all
        particles to the different coupling hierarchies so that a
        particle is considered to be in the highest hierarchy (i.e.,
        with lowest value) where it has an interaction.
        """
        
        # Find coupling orders in model
        coupling_orders = self.get('coupling_orders')
        # Loop through the different coupling hierarchy values, so we
        # start with the most dominant and proceed to the least dominant
        hierarchy = sorted(list(set([self.get('order_hierarchy')[k] for \
                                     k in coupling_orders])))

        # orders is a rising list of the lists of orders with a given hierarchy
        orders = []
        for value in hierarchy:
            orders.append([ k for (k, v) in \
                            self.get('order_hierarchy').items() if \
                            v == value ])

        # Extract the interaction that correspond to the different
        # coupling hierarchies, and the corresponding particles
        interactions = []
        particles = []
        for iorder, order in enumerate(orders):
            sum_orders = sum(orders[:iorder+1], [])
            sum_interactions = sum(interactions[:iorder], [])
            sum_particles = sum([list(p) for p in particles[:iorder]], [])
            # Append all interactions that have only orders with at least
            # this hierarchy
            interactions.append([i for i in self.get('interactions') if \
                                 not i in sum_interactions and \
                                 not any([k not in sum_orders for k in \
                                          i.get('orders').keys()])])
            # Append the corresponding particles, excluding the
            # particles that have already been added
            particles.append(set(sum([[p.get_pdg_code() for p in \
                                      inter.get('particles') if \
                                       p.get_pdg_code() not in sum_particles] \
                                      for inter in interactions[-1]], [])))

        return particles, hierarchy

    def get_max_WEIGHTED(self):
        """Return the maximum WEIGHTED order for any interaction in the model,
        for equivalent 3-particle vertices. Note that it can be fractional."""

        return max([inter.get_WEIGHTED_order(self) for inter in \
                        self.get('interactions')])
            

    def check_majoranas(self):
        """Return True if there is fermion flow violation, False otherwise"""

        if any([part.is_fermion() and part.get('self_antipart') \
                for part in self.get('particles')]):
            return True

        # No Majorana particles, but may still be fermion flow
        # violating interactions
        for inter in self.get('interactions'):
            # Do not look at UV Wfct renormalization counterterms
            if len(inter.get('particles'))==1:
                continue
            fermions = [p for p in inter.get('particles') if p.is_fermion()]
            for i in range(0, len(fermions), 2):
                if fermions[i].get('is_part') == \
                   fermions[i+1].get('is_part'):
                    # This is a fermion flow violating interaction
                    return True
        # No fermion flow violations
        return False

    def reset_dictionaries(self):
        """Reset all dictionaries and got_majoranas. This is necessary
        whenever the particle or interaction content has changed. If
        particles or interactions are set using the set routine, this
        is done automatically."""

        self['particle_dict'] = {}
        self['ref_dict_to0'] = {}
        self['got_majoranas'] = None
        self['interaction_dict'] = {}
        self['ref_dict_to1'] = {}
        self['ref_dict_to0'] = {}
        
    def pass_particles_name_in_mg_default(self):
        """Change the name of the particles such that all SM and MSSM particles
        follows the MG convention"""
        
        self.mg5_name = True
        
        # Check that default name/antiname is not already use 
        def check_name_free(self, name):
            """ check if name is not use for a particle in the model if it is 
            raise an MadGraph5error"""
            part = self['particles'].find_name(name)
            if part: 
                error_text = \
                '%s particles with pdg code %s is in conflict with MG ' + \
                'convention name for particle %s.\n Use -modelname in order ' + \
                'to use the particles name defined in the model and not the ' + \
                'MadGraph5_aMC@NLO convention'
                
                raise MadGraph5Error, error_text % \
                                     (part.get_name(), part.get_pdg_code(), pdg)                

        default = self.load_default_name()

        for pdg in default.keys():
            part = self.get_particle(pdg)
            if not part:
                continue
            antipart = self.get_particle(-pdg)
            name = part.get_name()
            if name != default[pdg]:
                check_name_free(self, default[pdg])
                if part.get('is_part'):
                    part.set('name', default[pdg])
                    if antipart:
                        antipart.set('name', default[pdg])
                    else:
                        part.set('antiname', default[pdg])                        
                else:
                    part.set('antiname', default[pdg])
                    if antipart:
                        antipart.set('antiname', default[pdg])
        
        #additional check for the Higgs in the mssm
        if self.get('name') == 'mssm' or self.get('name').startswith('mssm-'):
            part = self.get_particle(25)
            part.set('name', 'h1')
            part.set('antiname', 'h1')


            
    def change_parameter_name_with_prefix(self, prefix='mdl_'):
        """ Change all model parameter by a given prefix.
        Modify the parameter if some of them are identical up to the case"""
        
        lower_dict={}
        duplicate = set()
        keys = self.get('parameters').keys()
        for key in keys:
            for param in self['parameters'][key]:
                lower_name = param.name.lower()
                if not lower_name:
                    continue
                try:
                    lower_dict[lower_name].append(param)
                except KeyError:
                    lower_dict[lower_name] = [param]
                else:
                    duplicate.add(lower_name)
                    logger.debug('%s is defined both as lower case and upper case.' 
                                                                   % lower_name)
        
        if prefix == '' and  not duplicate:
            return
                
        re_expr = r'''\b(%s)\b'''
        to_change = []
        change={}
        # recast all parameter in prefix_XX
        for key in keys:
            for param in self['parameters'][key]:
                value = param.name.lower()
                if value in ['as','mu_r', 'zero','aewm1','g']:
                    continue
                elif value.startswith(prefix):
                    continue
                elif value in duplicate:
                    continue # handle later
                elif value:
                    change[param.name] = '%s%s' % (prefix,param.name)
                    to_change.append(param.name)
                    param.name = change[param.name]
            
        for value in duplicate:
            for i, var in enumerate(lower_dict[value]):
                to_change.append(var.name)
                new_name = '%s%s%s' % (prefix, var.name.lower(), 
                                                  ('__%d'%(i+1) if i>0 else ''))
                change[var.name] = new_name
                var.name = new_name
                to_change.append(var.name)
        assert 'zero' not in to_change
        replace = lambda match_pattern: change[match_pattern.groups()[0]]
        
        if not to_change:
            return
        
        if 'parameter_dict' in self:
            new_dict = dict( (change[name] if (name in change) else name, value) for
                             name, value in self['parameter_dict'].items())
            self['parameter_dict'] = new_dict
    
        if hasattr(self,'map_CTcoup_CTparam'):
            # If the map for the dependence of couplings to CTParameters has
            # been defined, we must apply the renaming there as well. 
            self.map_CTcoup_CTparam = dict( (coup_name, 
            [change[name] if (name in change) else name for name in params]) 
                  for coup_name, params in self.map_CTcoup_CTparam.items() )

        i=0
        while i*1000 <= len(to_change): 
            one_change = to_change[i*1000: min((i+1)*1000,len(to_change))]
            i+=1
            rep_pattern = re.compile('\\b%s\\b'% (re_expr % ('\\b|\\b'.join(one_change))))
            
            # change parameters
            for key in keys:
                if key == ('external',):
                    continue
                for param in self['parameters'][key]:
                    param.expr = rep_pattern.sub(replace, param.expr)
            # change couplings
            for key in self['couplings'].keys():
                for coup in self['couplings'][key]:
                    coup.expr = rep_pattern.sub(replace, coup.expr)

            # change form-factor
            ff = [l.formfactors for l in self['lorentz'] if hasattr(l, 'formfactors')]
            ff = set(sum(ff,[])) # here we have the list of ff used in the model
            for f in ff:
                f.value = rep_pattern.sub(replace, f.value)

            # change mass/width
            for part in self['particles']:
                if str(part.get('mass')) in one_change:
                    part.set('mass', rep_pattern.sub(replace, str(part.get('mass'))))
                if str(part.get('width')) in one_change:
                    part.set('width', rep_pattern.sub(replace, str(part.get('width'))))  
                if  hasattr(part, 'partial_widths'):
                    for key, value in part.partial_widths.items():    
                        part.partial_widths[key] = rep_pattern.sub(replace, value)
                
        #ensure that the particle_dict is up-to-date
        self['particle_dict'] =''
        self.get('particle_dict') 

        

    def get_first_non_pdg(self):
        """Return the first positive number that is not a valid PDG code"""
        return [c for c in range(1, len(self.get('particles')) + 1) if \
                c not in self.get('particle_dict').keys()][0]
                

    def write_param_card(self, filepath=None):
        """Write out the param_card, and return as string."""
        
        import models.write_param_card as writer
        if not filepath:
            out = StringIO.StringIO() # it's suppose to be written in a file
        else:
            out = filepath
        param = writer.ParamCardWriter(self, filepath=out)
        if not filepath:
            return out.getvalue()
        else:
            return param
        
    @ staticmethod
    def load_default_name():
        """ load the default for name convention """

        logger.info('Change particles name to pass to MG5 convention')    
        default = {}
        for line in open(os.path.join(MG5DIR, 'input', \
                                                 'particles_name_default.txt')):
            line = line.lstrip()
            if line.startswith('#'):
                continue
            
            args = line.split()
            if len(args) != 2:
                logger.warning('Invalid syntax in interface/default_name:\n %s' % line)
                continue
            default[int(args[0])] = args[1].lower()
        
        return default

    def change_electroweak_mode(self, mode):
        """Change the electroweak mode. The only valid mode now is external.
        Where in top of the default MW and sw2 are external parameters."""

        assert mode in ["external",set(['mz','mw','alpha'])]
        
        try:
            W = self.get('particle_dict')[24]
        except KeyError:
            raise InvalidCmd('No W particle in the model impossible to '+
                                                        'change the EW scheme!')
        
        if mode=='external':
            MW = self.get_parameter(W.get('mass'))
            if not isinstance(MW, ParamCardVariable):
                newMW = ParamCardVariable(MW.name, MW.value, 'MASS', [24])
                if not newMW.value:
                    newMW.value = 80.385
                #remove the old definition
                self.get('parameters')[MW.depend].remove(MW)
                # add the new one
                self.add_param(newMW, ['external'])
                
            # Now check for sw2. if not define bypass this
            try:
                sw2 = self.get_parameter('sw2')
            except KeyError:
                try:
                    sw2 = self.get_parameter('mdl_sw2')
                except KeyError:
                    sw2=None
    
            if sw2:
                newsw2 = ParamCardVariable(sw2.name,sw2.value, 'SMINPUTS', [4])
                if not newsw2.value:
                    newsw2.value = 0.222246485786
                #remove the old definition
                self.get('parameters')[sw2.depend].remove(sw2)
                # add the new one
                self.add_param(newsw2, ['external'])
            # Force a refresh of the parameter dictionary
            self.parameters_dict = None
            return True

        elif mode==set(['mz','mw','alpha']):
            # For now, all we support is to go from mz, Gf, alpha to mz, mw, alpha
            W = self.get('particle_dict')[24]
            mass = self.get_parameter(W.get('mass'))
            mass_expr = 'cmath.sqrt(%(prefix)sMZ__exp__2/2. + cmath.sqrt('+\
              '%(prefix)sMZ__exp__4/4. - (%(prefix)saEW*cmath.pi*%(prefix)s'+\
              'MZ__exp__2)/(%(prefix)sGf*%(prefix)ssqrt__2)))'
            if 'external' in mass.depend:
                # Nothing to be done
                return True
            match = False
            if mass.expr == mass_expr%{'prefix':''}:
                prefix = ''
                match = True
            elif mass.expr == mass_expr%{'prefix':'mdl_'}:
                prefix = 'mdl_'
                match = True
            if match:
                MW = ParamCardVariable(mass.name, mass.value, 'MASS', [24])
                if not MW.value:
                    MW.value = 80.385                
                self.get('parameters')[('external',)].append(MW)
                self.get('parameters')[mass.depend].remove(mass)
                # Make Gf an internal parameter
                new_param = ModelVariable('Gf',
                '-%(prefix)saEW*%(prefix)sMZ**2*cmath.pi/(cmath.sqrt(2)*%(MW)s**2*(%(MW)s**2 - %(prefix)sMZ**2))' %\
                {'MW': mass.name,'prefix':prefix}, 'complex', mass.depend)
                Gf = self.get_parameter('%sGf'%prefix)
                self.get('parameters')[('external',)].remove(Gf)
                self.add_param(new_param, ['%saEW'%prefix])
                # Force a refresh of the parameter dictionary
                self.parameters_dict = None
                return True
            else:
                return False

    def change_mass_to_complex_scheme(self, toCMS=True):
        """modify the expression changing the mass to complex mass scheme"""
        
        # 1) Change the 'CMSParam' of loop_qcd_qed model to 1.0 so as to remove
        #    the 'real' prefix fromall UVCT wf renormalization expressions.
        #    If toCMS is False, then it makes sure CMSParam is 0.0 and returns
        #    immediatly.
        # 2) Find All input parameter mass and width associated
        #   Add a internal parameter and replace mass with that param
        # 3) Find All mass fixed by the model and width associated
        #   -> Both need to be fixed with a real() /Imag()
        # 4) Find All width set by the model
        #   -> Need to be set with a real()
        # 5) Fix the Yukawa mass to the value of the complex mass/ real mass
        # 6) Loop through all expression and modify those accordingly
        #    Including all parameter expression as complex
        
        try:
            CMSParam = self.get_parameter('CMSParam')
        except KeyError:
            try:
                CMSParam = self.get_parameter('mdl_CMSParam')
            except KeyError:
                CMSParam = None
                
        # Handle the case where we want to make sure the CMS is turned off
        if not toCMS:
            if CMSParam:
                CMSParam.expr = '0.0'
            return            
        
        # Now handle the case where we want to turn to CMS.    
        if CMSParam:
            CMSParam.expr = '1.0'
        
        to_change = {}
        mass_widths = [] # parameter which should stay real
        for particle in self.get('particles'):
            m = particle.get('width')
            if m in mass_widths:
                continue
            mass_widths.append(particle.get('width'))
            mass_widths.append(particle.get('mass'))
            width = self.get_parameter(particle.get('width'))
            if (isinstance(width.value, (complex,float)) and abs(width.value)==0.0) or \
                                                    width.name.lower() =='zero':
                #everything is fine since the width is zero
                continue
            if not isinstance(width, ParamCardVariable):
                width.expr = 're(%s)' % width.expr
            mass = self.get_parameter(particle.get('mass'))
            if (isinstance(width.value, (complex,float)) and abs(width.value)!=0.0) or \
                                                    mass.name.lower() != 'zero':
                # special SM treatment to change the gauge scheme automatically.
                if particle.get('pdg_code') == 24 and isinstance(mass, 
                                                                 ModelVariable):
                    status = self.change_electroweak_mode(
                                                   set(['mz','mw','alpha']))
                    # Use the newly defined parameter for the W mass
                    mass = self.get_parameter(particle.get('mass'))
                    if not status:
                        logger.warning('The W mass is not an external '+
                        'parameter in this model and the automatic change of'+
                        ' electroweak scheme changed. This is not advised for '+
                                            'applying the complex mass scheme.')

                # Add A new parameter CMASS
                #first compute the dependencies (as,...)
                depend = list(set(mass.depend + width.depend))
                if len(depend)>1 and 'external' in depend:
                    depend.remove('external')
                depend = tuple(depend)
                if depend == ('external',):
                    depend = ()
                
                # Create the new parameter
                if isinstance(mass, ParamCardVariable):
                    New_param = ModelVariable('CMASS_'+mass.name,
                        'cmath.sqrt(%(mass)s**2 - complex(0,1) * %(mass)s * %(width)s)' \
                              % {'mass': mass.name, 'width': width.name}, 
                        'complex', depend)              
                else:
                    New_param = ModelVariable('CMASS_'+mass.name,
                        mass.expr, 'complex', depend)
                    # Modify the treatment of the width in this case
                    if not isinstance(width, ParamCardVariable):
                        width.expr = '- im(%s**2) / cmath.sqrt(re(%s**2))' % (mass.expr, mass.expr)
                    else:
                        # Remove external parameter from the param_card
                        New_width = ModelVariable(width.name,
                        '-1 * im(CMASS_%s**2) / %s' % (mass.name, mass.name), 'real', mass.depend)
                        self.get('parameters')[('external',)].remove(width)
                        self.add_param(New_param, (mass,))
                        self.add_param(New_width, (New_param,))
                        mass.expr = 'cmath.sqrt(re(%s**2))' % mass.expr                
                        to_change[mass.name] = New_param.name
                        continue                        
                        
                    mass.expr = 're(%s)' % mass.expr                
                self.add_param(New_param, (mass, width))
                to_change[mass.name] = New_param.name
        
        # Remove the Yukawa and fix those accordingly to the mass/complex mass
        yukawas = [p for p in self.get('parameters')[('external',)] 
                                              if p.lhablock.lower() == 'yukawa']
        for yukawa in yukawas:
            # clean the pevious parameter
            self.get('parameters')[('external',)].remove(yukawa)
            
            particle = self.get_particle(yukawa.lhacode[0])
            mass = self.get_parameter(particle.get('mass'))
            
            # add the new parameter in the correct category
            if mass.depend == ('external',):
                depend = ()
            else:
                depend = mass.depend
                
            New_param = ModelVariable(yukawa.name, mass.name, 'real', depend)
            
            # Add it in the model at the correct place (for the dependences)
            if mass.name in to_change:
                expr = 'CMASS_%s' % mass.name
            else:
                expr = mass.name
            param_depend = self.get_parameter(expr)
            self.add_param(New_param, [param_depend])
            
        if not to_change:
            return
            
            
        # So at this stage we still need to modify all parameters depending of
        # particle's mass. In addition all parameter (but mass/width/external 
        # parameter) should be pass in complex mode.
        pat = '|'.join(to_change.keys())
        pat = r'(%s)\b' % pat
        pat = re.compile(pat)
        def replace(match):
            return to_change[match.group()]
        
        # Modify the parameters
        for dep, list_param in self['parameters'].items():
            for param in list_param:
                if param.name.startswith('CMASS_') or param.name in mass_widths or\
                              isinstance(param, ParamCardVariable):
                    continue
                param.type = 'complex'
#                print param.expr,  to_change
                
                param.expr = pat.sub(replace, param.expr)
        
        # Modify the couplings        
        for dep, list_coup in self['couplings'].items():
            for coup in list_coup:                
                coup.expr = pat.sub(replace, coup.expr)
                
    def add_param(self, new_param, depend_param):
        """add the parameter in the list of parameter in a correct position"""
            
        pos = 0
        for i,param in enumerate(self.get('parameters')[new_param.depend]):
            if param.name in depend_param:
                pos = i + 1
        self.get('parameters')[new_param.depend].insert(pos, new_param)


    #def __repr__(self):
    #    """ """
    #    raise Exception
    #    return "Model(%s)" % self.get_name()
    #__str__ = __repr__
################################################################################
# Class for Parameter / Coupling
################################################################################
class ModelVariable(object):
    """A Class for storing the information about coupling/ parameter"""
    
    def __init__(self, name, expression, type, depend=()):
        """Initialize a new parameter/coupling"""
        
        self.name = name
        self.expr = expression # python expression
        self.type = type # real/complex
        self.depend = depend # depend on some other parameter -tuple-
        self.value = None
    
    def __eq__(self, other):
        """Object with same name are identical, If the object is a string we check
        if the attribute name is equal to this string"""
        
        try:
            return other.name == self.name
        except Exception:
            return other == self.name

class ParamCardVariable(ModelVariable):
    """ A class for storing the information linked to all the parameter 
    which should be define in the param_card.dat"""
    
    depend = ('external',)
    type = 'real'
    
    def __init__(self, name, value, lhablock, lhacode):
        """Initialize a new ParamCardVariable
        name: name of the variable
        value: default numerical value
        lhablock: name of the block in the param_card.dat
        lhacode: code associate to the variable
        """
        self.name = name
        self.value = value 
        self.lhablock = lhablock
        self.lhacode = lhacode


#===============================================================================
# Classes used in diagram generation and process definition:
#    Leg, Vertex, Diagram, Process
#===============================================================================

#===============================================================================
# Leg
#===============================================================================
class Leg(PhysicsObject):
    """Leg object: id (Particle), number, I/F state, flag from_group
    """

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

        self['id'] = 0
        self['number'] = 0
        # state: True = final, False = initial (boolean to save memory)
        self['state'] = True
        #self['loop_line'] = False
        self['loop_line'] = False
        # from_group: Used in diagram generation
        self['from_group'] = True
        # onshell: decaying leg (True), forbidden s-channel (False), none (None)
        self['onshell'] = None

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

        if name in ['id', 'number']:
            if not isinstance(value, int):
                raise self.PhysicsObjectError, \
                        "%s is not a valid integer for leg id" % str(value)

        if name == 'state':
            if not isinstance(value, bool):
                raise self.PhysicsObjectError, \
                        "%s is not a valid leg state (True|False)" % \
                                                                    str(value)

        if name == 'from_group':
            if not isinstance(value, bool) and value != None:
                raise self.PhysicsObjectError, \
                        "%s is not a valid boolean for leg flag from_group" % \
                                                                    str(value)

        if name == 'loop_line':
            if not isinstance(value, bool) and value != None:
                raise self.PhysicsObjectError, \
                    "%s is not a valid boolean for leg flag loop_line" % \
                                                                    str(value)

        if name == 'onshell':
            if not isinstance(value, bool) and value != None:
                raise self.PhysicsObjectError, \
                        "%s is not a valid boolean for leg flag onshell" % \
                                                                    str(value)
        return True

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

        return ['id', 'number', 'state', 'from_group', 'loop_line', 'onshell']

    def is_fermion(self, model):
        """Returns True if the particle corresponding to the leg is a
        fermion"""

        assert isinstance(model, Model), "%s is not a model" % str(model)

        return model.get('particle_dict')[self['id']].is_fermion()

    def is_incoming_fermion(self, model):
        """Returns True if leg is an incoming fermion, i.e., initial
        particle or final antiparticle"""

        assert isinstance(model, Model), "%s is not a model" % str(model)

        part = model.get('particle_dict')[self['id']]
        return part.is_fermion() and \
               (self.get('state') == False and part.get('is_part') or \
                self.get('state') == True and not part.get('is_part'))

    def is_outgoing_fermion(self, model):
        """Returns True if leg is an outgoing fermion, i.e., initial
        antiparticle or final particle"""

        assert isinstance(model, Model), "%s is not a model" % str(model)        
        
        part = model.get('particle_dict')[self['id']]
        return part.is_fermion() and \
               (self.get('state') == True and part.get('is_part') or \
                self.get('state') == False and not part.get('is_part'))

    # Helper function. We don't overload the == operator because it might be useful
    # to define it differently than that later.

    def same(self, leg):
        """ Returns true if the leg in argument has the same ID and the same numer """

        # In case we want to check this leg with an integer in the tagging procedure, 
        # then it only has to match the leg number.
        if isinstance(leg,int):
            if self['number']==leg:
                return True
            else:
                return False

        # If using a Leg object instead, we also want to compare the other relevant
        # properties.
        elif isinstance(leg, Leg):
            if self['id']==leg.get('id') and \
               self['number']==leg.get('number') and \
               self['loop_line']==leg.get('loop_line') :
                return True
            else:
                return False

        else :
            return False

    # Make sure sort() sorts lists of legs according to 'number'
    def __lt__(self, other):
        return self['number'] < other['number']

#===============================================================================
# LegList
#===============================================================================
class LegList(PhysicsObjectList):
    """List of Leg objects
    """

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

        return isinstance(obj, Leg)

    # Helper methods for diagram generation

    def from_group_elements(self):
        """Return all elements which have 'from_group' True"""

        return filter(lambda leg: leg.get('from_group'), self)

    def minimum_one_from_group(self):
        """Return True if at least one element has 'from_group' True"""

        return len(self.from_group_elements()) > 0

    def minimum_two_from_group(self):
        """Return True if at least two elements have 'from_group' True"""

        return len(self.from_group_elements()) > 1

    def can_combine_to_1(self, ref_dict_to1):
        """If has at least one 'from_group' True and in ref_dict_to1,
           return the return list from ref_dict_to1, otherwise return False"""
        if self.minimum_one_from_group():
            return ref_dict_to1.has_key(tuple(sorted([leg.get('id') for leg in self])))
        else:
            return False

    def can_combine_to_0(self, ref_dict_to0, is_decay_chain=False):
        """If has at least two 'from_group' True and in ref_dict_to0,
        
        return the vertex (with id from ref_dict_to0), otherwise return None

        If is_decay_chain = True, we only allow clustering of the
        initial leg, since we want this to be the last wavefunction to
        be evaluated.
        """
        if is_decay_chain:
            # Special treatment - here we only allow combination to 0
            # if the initial leg (marked by from_group = None) is
            # unclustered, since we want this to stay until the very
            # end.
            return any(leg.get('from_group') == None for leg in self) and \
                   ref_dict_to0.has_key(tuple(sorted([leg.get('id') \
                                                      for leg in self])))

        if self.minimum_two_from_group():
            return ref_dict_to0.has_key(tuple(sorted([leg.get('id') for leg in self])))
        else:
            return False

    def get_outgoing_id_list(self, model):
        """Returns the list of ids corresponding to the leglist with
        all particles outgoing"""

        res = []

        assert isinstance(model, Model), "Error! model not model"


        for leg in self:
            if leg.get('state') == False:
                res.append(model.get('particle_dict')[leg.get('id')].get_anti_pdg_code())
            else:
                res.append(leg.get('id'))

        return res
    
    def sort(self,*args, **opts):
        """Match with FKSLegList"""
        Opts=copy.copy(opts)
        if 'pert' in Opts.keys():
            del Opts['pert']
        return super(LegList,self).sort(*args, **Opts)


#===============================================================================
# MultiLeg
#===============================================================================
class MultiLeg(PhysicsObject):
    """MultiLeg object: ids (Particle or particles), I/F state
    """

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

        self['ids'] = []
        self['state'] = True

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

        if name == 'ids':
            if not isinstance(value, list):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list" % str(value)
            for i in value:
                if not isinstance(i, int):
                    raise self.PhysicsObjectError, \
                          "%s is not a valid list of integers" % str(value)

        if name == 'state':
            if not isinstance(value, bool):
                raise self.PhysicsObjectError, \
                        "%s is not a valid leg state (initial|final)" % \
                                                                    str(value)

        return True

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

        return ['ids', 'state']

#===============================================================================
# LegList
#===============================================================================
class MultiLegList(PhysicsObjectList):
    """List of MultiLeg objects
    """

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

        return isinstance(obj, MultiLeg)

#===============================================================================
# Vertex
#===============================================================================
class Vertex(PhysicsObject):
    """Vertex: list of legs (ordered), id (Interaction)
    """
    
    sorted_keys = ['id', 'legs']
    
    # This sets what are the ID's of the vertices that must be ignored for the
    # purpose of the multi-channeling. 0 and -1 are ID's of various technical
    # vertices which have no relevance from the perspective of the diagram 
    # topology, while -2 is the ID of a vertex that results from a shrunk loop
    # (for loop-induced integration with MadEvent) and one may or may not want
    # to consider these higher point loops for the purpose of the multi-channeling.
    # So, adding -2 to the list below makes sur that all loops are considered
    # for multichanneling.
    ID_to_veto_for_multichanneling = [0,-1,-2]
    
    # For loop-induced integration, considering channels from up to box loops 
    # typically leads to better efficiencies. Beyond that, it is detrimental 
    # because the phase-space generation is not suited to map contact interactions
    # This parameter controls up to how many legs should loop-induced diagrams
    # be considered for multichanneling.
    # Notice that, in the grouped subprocess case mode, if -2 is not added to 
    # the list ID_to_veto_for_multichanneling then all loop are considered by 
    # default and the constraint below is not applied.
    max_n_loop_for_multichanneling = 4
    
    def default_setup(self):
        """Default values for all properties"""

        # The 'id' of the vertex corresponds to the interaction ID it is made of.
        # Notice that this 'id' can take the special values :
        #  -1 : A two-point vertex which either 'sews' the two L-cut particles
        #       together or simply merges two wavefunctions to create an amplitude
        #       (in the case of tree-level diagrams).
        #  -2 : The id given to the ContractedVertices (i.e. a shrunk loop) so 
        #       that it can be easily identified when constructing the DiagramChainLinks.
        self['id'] = 0
        self['legs'] = LegList()

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

        if name == 'id':
            if not isinstance(value, int):
                raise self.PhysicsObjectError, \
                        "%s is not a valid integer for vertex id" % str(value)

        if name == 'legs':
            if not isinstance(value, LegList):
                raise self.PhysicsObjectError, \
                        "%s is not a valid LegList object" % str(value)

        return True

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

        return self.sorted_keys  #['id', 'legs']

    def nice_string(self):
        """return a nice string"""
        
        mystr = []
        for leg in self['legs']:
            mystr.append( str(leg['number']) + '(%s)' % str(leg['id']))
        mystr = '(%s,id=%s ,obj_id:%s)' % (', '.join(mystr), self['id'], id(self))
        
        return(mystr)


    def get_s_channel_id(self, model, ninitial):
        """Returns the id for the last leg as an outgoing
        s-channel. Returns 0 if leg is t-channel, or if identity
        vertex. Used to check for required and forbidden s-channel
        particles."""

        leg = self.get('legs')[-1]

        if ninitial == 1:
            # For one initial particle, all legs are s-channel
            # Only need to flip particle id if state is False            
            if leg.get('state') == True:
                return leg.get('id')
            else:
                return model.get('particle_dict')[leg.get('id')].\
                       get_anti_pdg_code()

        # Number of initial particles is at least 2
        if self.get('id') == 0 or \
           leg.get('state') == False:
            # identity vertex or t-channel particle
            return 0

        if leg.get('loop_line'):
            # Loop lines never count as s-channel
            return 0

        # Check if the particle number is <= ninitial
        # In that case it comes from initial and we should switch direction
        if leg.get('number') > ninitial:
            return leg.get('id')
        else:
            return model.get('particle_dict')[leg.get('id')].\
                       get_anti_pdg_code()

        ## Check if the other legs are initial or final.
        ## If the latter, return leg id, if the former, return -leg id
        #if self.get('legs')[0].get('state') == True:
        #    return leg.get('id')
        #else:
        #    return model.get('particle_dict')[leg.get('id')].\
        #               get_anti_pdg_code()

#===============================================================================
# VertexList
#===============================================================================
class VertexList(PhysicsObjectList):
    """List of Vertex objects
    """

    orders = {}

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

        return isinstance(obj, Vertex)

    def __init__(self, init_list=None, orders=None):
        """Creates a new list object, with an optional dictionary of
        coupling orders."""

        list.__init__(self)

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

        if isinstance(orders, dict):
            self.orders = orders

#===============================================================================
# ContractedVertex
#===============================================================================
class ContractedVertex(Vertex):
    """ContractedVertex: When contracting a loop to a given vertex, the created
    vertex object is then a ContractedVertex object which has additional 
    information with respect to a regular vertex object. For example, it contains
    the PDG of the particles attached to it. (necessary because the contracted
    vertex doesn't have an interaction ID which would allow to retrieve such
    information).
    """

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

        self['PDGs'] = []
        self['loop_tag'] = tuple()
        self['loop_orders'] = {}
        super(ContractedVertex, self).default_setup()

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

        if name == 'PDGs':
            if isinstance(value, list):
                for elem in value:
                    if not isinstance(elem,int):
                        raise self.PhysicsObjectError, \
                            "%s is not a valid integer for leg PDG" % str(elem)
            else:
                raise self.PhysicsObjectError, \
                  "%s is not a valid list for contracted vertex PDGs"%str(value)                
        if name == 'loop_tag':
            if isinstance(value, tuple):
                for elem in value:
                    if not (isinstance(elem,int) or isinstance(elem,tuple)):
                        raise self.PhysicsObjectError, \
                          "%s is not a valid int or tuple for loop tag element"%str(elem)
            else:
                raise self.PhysicsObjectError, \
                  "%s is not a valid tuple for a contracted vertex loop_tag."%str(value)
        if name == 'loop_orders':
            Interaction.filter(Interaction(), 'orders', value)
        else:
            return super(ContractedVertex, self).filter(name, value)

        return True

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

        return super(ContractedVertex, self).get_sorted_keys()+['PDGs']

#===============================================================================
# Diagram
#===============================================================================
class Diagram(PhysicsObject):
    """Diagram: list of vertices (ordered)
    """

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

        self['vertices'] = VertexList()
        self['orders'] = {}

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

        if name == 'vertices':
            if not isinstance(value, VertexList):
                raise self.PhysicsObjectError, \
                        "%s is not a valid VertexList object" % str(value)

        if name == 'orders':
            Interaction.filter(Interaction(), 'orders', value)

        return True

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

        return ['vertices', 'orders']
    
    def nice_string(self):
        """Returns a nicely formatted string of the diagram content."""

        pass_sanity = True
        if self['vertices']:
            mystr = '('
            for vert in self['vertices']:
                used_leg = [] 
                mystr = mystr + '('
                for leg in vert['legs'][:-1]:
                    mystr = mystr + str(leg['number']) + '(%s)' % str(leg['id']) + ','
                    used_leg.append(leg['number'])
                if __debug__ and len(used_leg) != len(set(used_leg)):
                    pass_sanity = False
                    responsible = id(vert)
                    
                if self['vertices'].index(vert) < len(self['vertices']) - 1:
                    # Do not want ">" in the last vertex
                    mystr = mystr[:-1] + '>'
                mystr = mystr + str(vert['legs'][-1]['number']) + '(%s)' % str(vert['legs'][-1]['id']) + ','
                mystr = mystr + 'id:' + str(vert['id']) + '),'
                                
            mystr = mystr[:-1] + ')'
            mystr += " (%s)" % (",".join(["%s=%d" % (key, self['orders'][key]) \
                                     for key in sorted(self['orders'].keys())]))
            
            if not pass_sanity:
                raise Exception, "invalid diagram: %s. vert_id: %s" % (mystr, responsible) 
                
            return mystr
        else:
            return '()'
    
    def calculate_orders(self, model):
        """Calculate the actual coupling orders of this diagram. Note
        that the special order WEIGTHED corresponds to the sum of
        hierarchys for the couplings."""

        coupling_orders = dict([(c, 0) for c in model.get('coupling_orders')])
        weight = 0
        for vertex in self['vertices']:
            if vertex.get('id') in [0,-1]: continue
            if vertex.get('id') == -2:
                couplings = vertex.get('loop_orders')
            else:
                couplings = model.get('interaction_dict')[vertex.get('id')].\
                                                                   get('orders')
            for coupling in couplings:
                coupling_orders[coupling] += couplings[coupling]
            weight += sum([model.get('order_hierarchy')[c]*n for \
                              (c,n) in couplings.items()])
        coupling_orders['WEIGHTED'] = weight
        self.set('orders', coupling_orders)

    def pass_squared_order_constraints(self, diag_multiplier, squared_orders,
                                                               sq_orders_types):
        """ Returns wether the contributiong consisting in the current diagram 
        multiplied by diag_multiplier passes the *positive* squared_orders 
        specified ( a dictionary ) of types sq_order_types (a dictionary whose 
        values are the relational operator used to define the constraint of the
        order in key)."""
        
        for order, value in squared_orders.items():
            if value<0:
                continue
            combined_order = self.get_order(order) + \
                                                diag_multiplier.get_order(order)
            if ( sq_orders_types[order]=='==' and combined_order != value ) or \
               ( sq_orders_types[order] in ['=', '<='] and combined_order > value) or \
               ( sq_orders_types[order]=='>' and combined_order <= value) :
                return False
        return True

    def get_order(self, order):
        """Return the order of this diagram. It returns 0 if it is not present."""

        try:
            return self['orders'][order]
        except Exception:
            return 0

    def get_contracted_loop_diagram(self, struct_rep=None):
        """ Returns a Diagram which correspond to the loop diagram with the 
        loop shrunk to a point. Of course for a instance of base_objects.Diagram
        one must simply return self."""
        
        return self
        
    def get_external_legs(self):
        """ Return the list of external legs of this diagram """
        
        external_legs = LegList([])
        for leg in sum([vert.get('legs') for vert in self.get('vertices')],[]):
            if not leg.get('number') in [l.get('number') for l in external_legs]:
               external_legs.append(leg) 
               
        return external_legs
        
    def renumber_legs(self, perm_map, leg_list):
        """Renumber legs in all vertices according to perm_map"""

        vertices = VertexList()
        min_dict = copy.copy(perm_map)
        # Dictionary from leg number to state
        state_dict = dict([(l.get('number'), l.get('state')) for l in leg_list])
        # First renumber all legs in the n-1->1 vertices
        for vertex in self.get('vertices')[:-1]:
            vertex = copy.copy(vertex)
            leg_list = LegList([copy.copy(l) for l in vertex.get('legs')])
            for leg in leg_list[:-1]:
                leg.set('number', min_dict[leg.get('number')])
                leg.set('state', state_dict[leg.get('number')])
            min_number = min([leg.get('number') for leg in leg_list[:-1]])
            leg = leg_list[-1]
            min_dict[leg.get('number')] = min_number
            # resulting leg is initial state if there is exactly one
            # initial state leg among the incoming legs
            state_dict[min_number] = len([l for l in leg_list[:-1] if \
                                          not l.get('state')]) != 1
            leg.set('number', min_number)
            leg.set('state', state_dict[min_number])
            vertex.set('legs', leg_list)
            vertices.append(vertex)
        # Now renumber the legs in final vertex
        vertex = copy.copy(self.get('vertices')[-1])
        leg_list = LegList([copy.copy(l) for l in vertex.get('legs')])
        for leg in leg_list:
            leg.set('number', min_dict[leg.get('number')])
            leg.set('state', state_dict[leg.get('number')])
        vertex.set('legs', leg_list)
        vertices.append(vertex)
        # Finally create new diagram
        new_diag = copy.copy(self)
        new_diag.set('vertices', vertices)
        state_dict = {True:'T',False:'F'}
        return new_diag

    def get_vertex_leg_numbers(self, 
                        veto_inter_id=Vertex.ID_to_veto_for_multichanneling,
                        max_n_loop=0):
        """Return a list of the number of legs in the vertices for
        this diagram. 
        This function is only used for establishing the multi-channeling, so that
        we exclude from it all the fake vertices and the vertices resulting from
        shrunk loops (id=-2)"""


        if max_n_loop == 0:
            max_n_loop = Vertex.max_n_loop_for_multichanneling
        
        res = [len(v.get('legs')) for v in self.get('vertices') if (v.get('id') \
                                  not in veto_inter_id) or (v.get('id')==-2 and 
                                                 len(v.get('legs'))>max_n_loop)]

        return res
    
    def get_num_configs(self, model, ninitial):
        """Return the maximum number of configs from this diagram,
        given by 2^(number of non-zero width s-channel propagators)"""

        s_channels = [v.get_s_channel_id(model,ninitial) for v in \
                              self.get('vertices')[:-1]]
        num_props = len([i for i in s_channels if i != 0 and \
                         model.get_particle(i).get('width').lower() != 'zero'])
        
        if num_props < 1:
            return 1
        else:
            return 2**num_props
        
    def get_flow_charge_diff(self, model):
        """return the difference of total diff of charge occuring on the 
        lofw of the initial parton. return [None,None] if the two initial parton
        are connected and the (partial) value if None if the initial parton is 
        not a fermiom"""
        
        import madgraph.core.drawing as drawing
        drawdiag = drawing.FeynmanDiagram(self, model)
        drawdiag.load_diagram()
        out = []
        
        for v in drawdiag.initial_vertex:
            init_part = v.lines[0]
            if not init_part.is_fermion():
                out.append(None)
                continue
            
            init_charge = model.get_particle(init_part.id).get('charge')
            
            l_last = init_part
            v_last = v
            vcurrent = l_last.end
            if vcurrent == v:
                vcurrent = l_last.begin
            security =0
            while not vcurrent.is_external():
                if security > 1000:
                    raise Exception, 'wrong diagram'
                next_l = [l for l in vcurrent.lines if l is not l_last and l.is_fermion()][0]
                next_v = next_l.end
                if next_v == vcurrent:
                    next_v = next_l.begin
                l_last, vcurrent = next_l, next_v
            if vcurrent in drawdiag.initial_vertex:
                return [None, None]
            
            out.append(model.get_particle(l_last.id).get('charge') - init_charge)    
        return out
                

#===============================================================================
# DiagramList
#===============================================================================
class DiagramList(PhysicsObjectList):
    """List of Diagram objects
    """

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

        return isinstance(obj, Diagram)

    def nice_string(self, indent=0):
        """Returns a nicely formatted string"""
        mystr = " " * indent + str(len(self)) + ' diagrams:\n'
        for i, diag in enumerate(self):
            mystr = mystr + " " * indent + str(i+1) + "  " + \
                    diag.nice_string() + '\n'
        return mystr[:-1]

    # Helper function

    def get_max_order(self,order):
        """ Return the order of the diagram in the list with the maximum coupling
        order for the coupling specified """
        max_order=-1

        for diag in self:
            if order in diag['orders'].keys():
                if max_order==-1 or diag['orders'][order] > max_order:
                    max_order = diag['orders'][order]

        return max_order

    def apply_negative_sq_order(self, ref_diag_list, order, value, order_type):
        """ This function returns a fitlered version of the diagram list self
        which satisfy the negative squared_order constraint 'order' with negative
        value 'value' and of type 'order_type', assuming that the diagram_list
        it must be squared against is 'reg_diag_list'. It also returns the
        new postive target squared order which correspond to this negative order
        constraint. Example: u u~ > d d~ QED^2<=-2 means that one wants to
        pick terms only up to the the next-to-leading order contributiong in QED,
        which is QED=2 in this case, so that target_order=4 is returned."""
        
        # First we must compute all contributions to that order
        target_order = min(ref_diag_list.get_order_values(order))+\
                                  min(self.get_order_values(order))+2*(-value-1)
        
        new_list = self.apply_positive_sq_orders(ref_diag_list, 
                                       {order:target_order}, {order:order_type})
        
        return new_list, target_order
        
    def apply_positive_sq_orders(self, ref_diag_list, sq_orders, sq_order_types):
        """ This function returns a filtered version of self which contain
        only the diagram which satisfy the positive squared order constraints
        sq_orders of type sq_order_types and assuming that the diagrams are
        multiplied with those of the reference diagram list ref_diag_list."""
                
        new_diag_list = DiagramList()
        for tested_diag in self:
            for ref_diag in ref_diag_list:
                if tested_diag.pass_squared_order_constraints(ref_diag,
                                                      sq_orders,sq_order_types):
                    new_diag_list.append(tested_diag)
                    break
        return new_diag_list

    def filter_constrained_orders(self, order, value, operator):
        """ This function modifies the current object and remove the diagram
        which do not obey the condition """ 
        
        new = []
        for tested_diag in self:
            if operator == '==':
                if tested_diag['orders'][order] == value:
                    new.append(tested_diag)
            elif operator == '>':
                if tested_diag['orders'][order] > value:
                    new.append(tested_diag)                
        self[:] = new
        return self


    def get_min_order(self,order):
        """ Return the order of the diagram in the list with the mimimum coupling
        order for the coupling specified """
        min_order=-1
        for diag in self:
            if order in diag['orders'].keys():
                if min_order==-1 or diag['orders'][order] < min_order:
                    min_order = diag['orders'][order]
            else:
                return 0

        return min_order

    def get_order_values(self, order):
        """ Return the list of possible values appearing in the diagrams of this
        list for the order given in argument """

        values=set([])
        for diag in self:
            if order in diag['orders'].keys():
                values.add(diag['orders'][order])
            else:
                values.add(0)  

        return list(values)

#===============================================================================
# Process
#===============================================================================
class Process(PhysicsObject):
    """Process: list of legs (ordered)
                dictionary of orders
                model
                process id
    """

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

        self['legs'] = LegList()
        # These define the orders restrict the born and loop amplitudes.
        self['orders'] = {}
        self['model'] = Model()
        # Optional number to identify the process
        self['id'] = 0
        self['uid'] = 0 # should be a uniq id number
        # Required s-channels are given as a list of id lists. Only
        # diagrams with all s-channels in any of the lists are
        # allowed. This enables generating e.g. Z/gamma as s-channel
        # propagators.
        self['required_s_channels'] = []
        self['forbidden_onsh_s_channels'] = []
        self['forbidden_s_channels'] = []
        self['forbidden_particles'] = []
        self['is_decay_chain'] = False
        self['overall_orders'] = {}
        # Decay chain processes associated with this process
        self['decay_chains'] = ProcessList()
        # Legs with decay chains substituted in
        self['legs_with_decays'] = LegList()
        # Loop particles if the process is to be computed at NLO
        self['perturbation_couplings']=[]        
        # These orders restrict the order of the squared amplitude.
        # This dictionary possibly contains a key "WEIGHTED" which
        # gives the upper bound for the total weighted order of the
        # squared amplitude.
        self['squared_orders'] = {}
        # The squared order (sqorders) constraints above can either be upper 
        # bound (<=) or exact match (==) depending on how they were specified
        # in the user input. This choice is stored in the dictionary below.
        # Notice that the upper bound is the default
        self['sqorders_types'] = {}
        # other type of constraint at amplitude level
        self['constrained_orders'] = {} # {QED: (4,'>')}
        self['has_born'] = True
        # The NLO_mode is always None for a tree-level process and can be
        # 'all', 'real', 'virt' for a loop process.
        self['NLO_mode'] = 'tree'
        # in the context of QED or QED+QCD perturbation, it is useful to
        # keep track of the orders that have been explicitly asked by the 
        # user, because other borns will appear used for the subtraction
        # of singularities
        self['born_orders'] = {}
        # The user might want to have the individual matrix element evaluations
        # for specific values of the coupling orders. The list below specifies
        # what are the coupling names which need be individually treated.
        # For example, for the process p p > j j [] QED=2 (QED=2 is 
        # then a squared order constraint), then QED will appear in the 
        # 'split_orders' list so that the subroutine in matrix.f return the
        # evaluation of the matrix element individually for the pure QCD 
        # contribution 'QCD=4 QED=0', the pure interference 'QCD=2 QED=2' and
        # the pure QED contribution of order 'QCD=0 QED=4'.
        self['split_orders'] = []

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

        if name in ['legs', 'legs_with_decays'] :
            if not isinstance(value, LegList):
                raise self.PhysicsObjectError, \
                        "%s is not a valid LegList object" % str(value)

        if name in ['orders', 'overall_orders','squared_orders', 'born_orders']:
            Interaction.filter(Interaction(), 'orders', value)

        if name == 'constrained_orders':
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                        "%s is not a valid dictionary" % str(value)            

        if name == 'sqorders_types':
            if not isinstance(value, dict):
                raise self.PhysicsObjectError, \
                        "%s is not a valid dictionary" % str(value)
            for order in value.keys()+value.values():
                if not isinstance(order, str):
                    raise self.PhysicsObjectError, \
                          "%s is not a valid string" % str(value)

        if name == 'split_orders':
            if not isinstance(value, list):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list" % str(value)
            for order in value:
                if not isinstance(order, str):
                    raise self.PhysicsObjectError, \
                          "%s is not a valid string" % str(value)

        if name == 'model':
            if not isinstance(value, Model):
                raise self.PhysicsObjectError, \
                        "%s is not a valid Model object" % str(value)
        if name in ['id', 'uid']:
            if not isinstance(value, int):
                raise self.PhysicsObjectError, \
                    "Process %s %s is not an integer" % (name, repr(value))

        if name == 'required_s_channels':
            if not isinstance(value, list):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list" % str(value)
            for l in value:
                if not isinstance(l, list):
                    raise self.PhysicsObjectError, \
                          "%s is not a valid list of lists" % str(value)
                for i in l:
                    if not isinstance(i, int):
                        raise self.PhysicsObjectError, \
                              "%s is not a valid list of integers" % str(l)
                    if i == 0:
                        raise self.PhysicsObjectError, \
                          "Not valid PDG code %d for s-channel particle" % i

        if name in ['forbidden_onsh_s_channels', 'forbidden_s_channels']:
            if not isinstance(value, list):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list" % str(value)
            for i in value:
                if not isinstance(i, int):
                    raise self.PhysicsObjectError, \
                          "%s is not a valid list of integers" % str(value)
                if i == 0:
                    raise self.PhysicsObjectError, \
                      "Not valid PDG code %d for s-channel particle" % str(value)

        if name == 'forbidden_particles':
            if not isinstance(value, list):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list" % str(value)
            for i in value:
                if not isinstance(i, int):
                    raise self.PhysicsObjectError, \
                          "%s is not a valid list of integers" % str(value)
                if i <= 0:
                    raise self.PhysicsObjectError, \
                      "Forbidden particles should have a positive PDG code" % str(value)

        if name == 'perturbation_couplings':
            if not isinstance(value, list):
                raise self.PhysicsObjectError, \
                        "%s is not a valid list" % str(value)
            for order in value:
                if not isinstance(order, str):
                    raise self.PhysicsObjectError, \
                          "%s is not a valid string" % str(value)

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

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

        if name == 'decay_chains':
            if not isinstance(value, ProcessList):
                raise self.PhysicsObjectError, \
                        "%s is not a valid ProcessList" % str(value)

        if name == 'NLO_mode':
            import madgraph.interface.madgraph_interface as mg
            if value not in mg.MadGraphCmd._valid_nlo_modes:
                raise self.PhysicsObjectError, \
                        "%s is not a valid NLO_mode" % str(value)
        return True

    def has_multiparticle_label(self):
        """ A process, not being a ProcessDefinition never carries multiple
        particles labels"""
        
        return False

    def set(self, name, value):
        """Special set for forbidden particles - set to abs value."""

        if name == 'forbidden_particles':
            try:
                value = [abs(i) for i in value]
            except Exception:
                pass

        if name == 'required_s_channels':
            # Required s-channels need to be a list of lists of ids
            if value and isinstance(value, list) and \
               not isinstance(value[0], list):
                value = [value]

        return super(Process, self).set(name, value) # call the mother routine

    def get_squared_order_type(self, order):
        """ Return what kind of squared order constraint was specified for the
        order 'order'."""

        if order in self['sqorders_types'].keys():
            return self['sqorders_types'][order]
        else:
            # Default behavior '=' is interpreted as upper bound '<='
            return '='

    def get(self, name):
        """Special get for legs_with_decays"""
        
        if name == 'legs_with_decays':
            self.get_legs_with_decays()

        if name == 'sqorders_types':
            # We must make sure that there is a type for each sqorder defined
            for order in self['squared_orders'].keys():
                if order not in self['sqorders_types']:
                    # Then assign its type to the default '='
                    self['sqorders_types'][order]='='
                    
        return super(Process, self).get(name) # call the mother routine

    

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

        return ['legs', 'orders', 'overall_orders', 'squared_orders',
                'constrained_orders',
                'model', 'id', 'required_s_channels', 
                'forbidden_onsh_s_channels', 'forbidden_s_channels',
                'forbidden_particles', 'is_decay_chain', 'decay_chains',
                'legs_with_decays', 'perturbation_couplings', 'has_born', 
                'NLO_mode', 'split_orders', 'born_orders']

    def nice_string(self, indent=0, print_weighted = True, prefix=True):
        """Returns a nicely formated string about current process
        content. Since the WEIGHTED order is automatically set and added to 
        the user-defined list of orders, it can be ommitted for some info
        displays."""

        if isinstance(prefix, bool) and prefix:
            mystr = " " * indent + "Process: "
        elif isinstance(prefix, str):
            mystr = prefix
        else:
            mystr = ""
        prevleg = None
        for leg in self['legs']:
            mypart = self['model'].get('particle_dict')[leg['id']]
            if prevleg and prevleg['state'] == False \
                   and leg['state'] == True:
                # Separate initial and final legs by >
                mystr = mystr + '> '
                # Add required s-channels
                if self['required_s_channels'] and \
                       self['required_s_channels'][0]:
                    mystr += "|".join([" ".join([self['model'].\
                                       get('particle_dict')[req_id].get_name() \
                                                for req_id in id_list]) \
                                    for id_list in self['required_s_channels']])
                    mystr = mystr + ' > '

            mystr = mystr + mypart.get_name() + ' '
            #mystr = mystr + '(%i) ' % leg['number']
            prevleg = leg

        # Add orders
        if self['orders']:
            to_add = []
            for key in sorted(self['orders'].keys()):
                if not print_weighted and key == 'WEIGHTED':
                    continue
                value = int(self['orders'][key])
                if key in self['squared_orders']:
                    if self.get_squared_order_type(key) in ['<=', '==', '='] and \
                        self['squared_orders'][key] == value:
                        continue 
                    if self.get_squared_order_type(key) in ['>'] and value == 99:
                        continue
                if key in self['constrained_orders']:
                    if value == self['constrained_orders'][key][0] and\
                       self['constrained_orders'][key][1] in ['=', '<=', '==']:
                        continue
                if value == 0:
                    to_add.append('%s=0' % key)
                else:
                    to_add.append('%s<=%s' % (key,value))
                 
            if to_add:
                mystr = mystr + " ".join(to_add) + ' '

        if self['constrained_orders']:
            mystr = mystr + " ".join('%s%s%d' % (key, 
              self['constrained_orders'][key][1], self['constrained_orders'][key][0]) 
                    for key in sorted(self['constrained_orders'].keys()))  + ' '

        # Add perturbation_couplings
        if self['perturbation_couplings']:
            mystr = mystr + '[ '
            if self['NLO_mode']!='tree':
                if self['NLO_mode']=='virt' and not self['has_born']:
                    mystr = mystr + 'sqrvirt = '
                else:
                    mystr = mystr + self['NLO_mode'] + ' = '
            for order in self['perturbation_couplings']:
                mystr = mystr + order + ' '
            mystr = mystr + '] '

        # Add squared orders
        if self['squared_orders']:
            to_add = []
            for key in sorted(self['squared_orders'].keys()):
                if not print_weighted and key == 'WEIGHTED':
                    continue
                if key in self['constrained_orders']:
                    if self['constrained_orders'][key][0] == self['squared_orders'][key]/2 and \
                       self['constrained_orders'][key][1] == self.get_squared_order_type(key):
                        continue
                to_add.append(key + '^2%s%d'%\
                (self.get_squared_order_type(key),self['squared_orders'][key]))
            
            if to_add:
                mystr = mystr + " ".join(to_add) + ' '
            

        # Add forbidden s-channels
        if self['forbidden_onsh_s_channels']:
            mystr = mystr + '$ '
            for forb_id in self['forbidden_onsh_s_channels']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        # Add double forbidden s-channels
        if self['forbidden_s_channels']:
            mystr = mystr + '$$ '
            for forb_id in self['forbidden_s_channels']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        # Add forbidden particles
        if self['forbidden_particles']:
            mystr = mystr + '/ '
            for forb_id in self['forbidden_particles']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        # Remove last space
        mystr = mystr[:-1]

        if self.get('id') or self.get('overall_orders'):
            mystr += " @%d" % self.get('id')
            if self.get('overall_orders'):
                mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \
                       for key in sorted(self['orders'])]) + ' '
        
        if not self.get('decay_chains'):
            return mystr

        for decay in self['decay_chains']:
            mystr = mystr + '\n' + \
                    decay.nice_string(indent + 2).replace('Process', 'Decay')

        return mystr

    def input_string(self):
        """Returns a process string corresponding to the input string
        in the command line interface."""

        mystr = ""
        prevleg = None

        for leg in self['legs']:
            mypart = self['model'].get('particle_dict')[leg['id']]
            if prevleg and prevleg['state'] == False \
                   and leg['state'] == True:
                # Separate initial and final legs by ">"
                mystr = mystr + '> '
                # Add required s-channels
                if self['required_s_channels'] and \
                       self['required_s_channels'][0]:
                    mystr += "|".join([" ".join([self['model'].\
                                       get('particle_dict')[req_id].get_name() \
                                                for req_id in id_list]) \
                                    for id_list in self['required_s_channels']])
                    mystr = mystr + '> '

            mystr = mystr + mypart.get_name() + ' '
            #mystr = mystr + '(%i) ' % leg['number']
            prevleg = leg

        if self['orders']:
            mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
                       for key in self['orders']]) + ' '

        # Add squared orders
        if self['squared_orders']:
            mystr = mystr + " ".join([key + '^2=' + repr(self['squared_orders'][key]) \
                       for key in self['squared_orders']]) + ' '

        # Add perturbation orders
        if self['perturbation_couplings']:
            mystr = mystr + '[ '
            if self['NLO_mode']:
                mystr = mystr + self['NLO_mode']
                if not self['has_born']:
                    mystr = mystr + '^2'
                mystr = mystr + '= '
                
            for order in self['perturbation_couplings']:
                mystr = mystr + order + ' '
            mystr = mystr + '] '


        # Add forbidden s-channels
        if self['forbidden_onsh_s_channels']:
            mystr = mystr + '$ '
            for forb_id in self['forbidden_onsh_s_channels']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        # Add double forbidden s-channels
        if self['forbidden_s_channels']:
            mystr = mystr + '$$ '
            for forb_id in self['forbidden_s_channels']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        # Add forbidden particles
        if self['forbidden_particles']:
            mystr = mystr + '/ '
            for forb_id in self['forbidden_particles']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        # Remove last space
        mystr = mystr[:-1]

        if self.get('overall_orders'):
            mystr += " @%d" % self.get('id')
            if self.get('overall_orders'):
                mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \
                       for key in sorted(self['orders'])]) + ' '
        
        if not self.get('decay_chains'):
            return mystr

        for decay in self['decay_chains']:
            paren1 = ''
            paren2 = ''
            if decay.get('decay_chains'):
                paren1 = '('
                paren2 = ')'
            mystr += ', ' + paren1 + decay.input_string() + paren2

        return mystr

    def base_string(self):
        """Returns a string containing only the basic process (w/o decays)."""

        mystr = ""
        prevleg = None
        for leg in self.get_legs_with_decays():
            mypart = self['model'].get('particle_dict')[leg['id']]
            if prevleg and prevleg['state'] == False \
                   and leg['state'] == True:
                # Separate initial and final legs by ">"
                mystr = mystr + '> '
            mystr = mystr + mypart.get_name() + ' '
            prevleg = leg

        # Remove last space
        return mystr[:-1]

    def shell_string(self, schannel=True, forbid=True, main=True, pdg_order=False,
                                                                print_id = True):
        """Returns process as string with '~' -> 'x', '>' -> '_',
        '+' -> 'p' and '-' -> 'm', including process number,
        intermediate s-channels and forbidden particles,
        pdg_order allow to order to leg order by pid."""

        mystr = ""
        if not self.get('is_decay_chain') and print_id:
            mystr += "%d_" % self['id']
        
        prevleg = None
        if pdg_order:
            legs = [l for l in self['legs'][1:]]
            def order_leg(l1,l2):
                id1 = l1.get('id')
                id2 = l2.get('id')
                return id2-id1
            legs.sort(cmp=order_leg)
            legs.insert(0, self['legs'][0])
        else:
            legs = self['legs']
        
        
        for leg in legs:
            mypart = self['model'].get('particle_dict')[leg['id']]
            if prevleg and prevleg['state'] == False \
                   and leg['state'] == True:
                # Separate initial and final legs by ">"
                mystr = mystr + '_'
                # Add required s-channels
                if self['required_s_channels'] and \
                       self['required_s_channels'][0] and schannel:
                    mystr += "_or_".join(["".join([self['model'].\
                                       get('particle_dict')[req_id].get_name() \
                                                for req_id in id_list]) \
                                    for id_list in self['required_s_channels']])
                    mystr = mystr + '_'
            if mypart['is_part']:
                mystr = mystr + mypart['name']
            else:
                mystr = mystr + mypart['antiname']
            prevleg = leg

        # Check for forbidden particles
        if self['forbidden_particles'] and forbid:
            mystr = mystr + '_no_'
            for forb_id in self['forbidden_particles']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name()

        # Replace '~' with 'x'
        mystr = mystr.replace('~', 'x')
        # Replace '+' with 'p'
        mystr = mystr.replace('+', 'p')
        # Replace '-' with 'm'
        mystr = mystr.replace('-', 'm')
        # Just to be safe, remove all spaces
        mystr = mystr.replace(' ', '')

        for decay in self.get('decay_chains'):
            mystr = mystr + "_" + decay.shell_string(schannel,forbid, main=False,
                                                     pdg_order=pdg_order)

        # Too long name are problematic so restrict them to a maximal of 70 char
        if len(mystr) > 64 and main:
            if schannel and forbid:
                out = self.shell_string(True, False, True, pdg_order)
            elif schannel:
                out = self.shell_string(False, False, True, pdg_order)
            else:
                out = mystr[:64]
            if not out.endswith('_%s' % self['uid']):    
                out += '_%s' % self['uid']
            return out

        return mystr

    def shell_string_v4(self):
        """Returns process as v4-compliant string with '~' -> 'x' and
        '>' -> '_'"""

        mystr = "%d_" % self['id']
        prevleg = None
        for leg in self.get_legs_with_decays():
            mypart = self['model'].get('particle_dict')[leg['id']]
            if prevleg and prevleg['state'] == False \
                   and leg['state'] == True:
                # Separate initial and final legs by ">"
                mystr = mystr + '_'
            if mypart['is_part']:
                mystr = mystr + mypart['name']
            else:
                mystr = mystr + mypart['antiname']
            prevleg = leg

        # Replace '~' with 'x'
        mystr = mystr.replace('~', 'x')
        # Just to be safe, remove all spaces
        mystr = mystr.replace(' ', '')

        return mystr

    # Helper functions

    def are_negative_orders_present(self):
        """ Check iteratively that no coupling order constraint include negative
        values."""

        if any(val<0 for val in self.get('orders').values()+\
                                           self.get('squared_orders').values()):
            return True
        
        for procdef in self['decay_chains']:
            if procdef.are_negative_orders_present():
                return True

        return False

    def are_decays_perturbed(self):
        """ Check iteratively that the decayed processes are not perturbed """
        
        for procdef in self['decay_chains']:
            if procdef['perturbation_couplings'] or procdef.are_decays_perturbed():
                return True
        return False
    
    def decays_have_squared_orders(self):
        """ Check iteratively that the decayed processes are not perturbed """
        
        for procdef in self['decay_chains']:
            if procdef['squared_orders']!={} or procdef.decays_have_squared_orders():
                return True
        return False
    
    def get_ninitial(self):
        """Gives number of initial state particles"""

        return len(filter(lambda leg: leg.get('state') == False,
                           self.get('legs')))

    def get_initial_ids(self):
        """Gives the pdg codes for initial state particles"""

        return [leg.get('id') for leg in \
                filter(lambda leg: leg.get('state') == False,
                       self.get('legs'))]

    def get_initial_pdg(self, number):
        """Return the pdg codes for initial state particles for beam number"""

        legs = filter(lambda leg: leg.get('state') == False and\
                       leg.get('number') == number,
                       self.get('legs'))
        if not legs:
            return None
        else:
            return legs[0].get('id')
        
    def get_initial_final_ids(self):
        """return a tuple of two tuple containing the id of the initial/final
           state particles. Each list is ordered"""
           
        initial = []
        final = [l.get('id') for l in self.get('legs')\
              if l.get('state') or initial.append(l.get('id'))]
        initial.sort()
        final.sort()
        return (tuple(initial), tuple(final))
    
    def get_final_ids_after_decay(self):
        """Give the pdg code of the process including decay"""
        
        finals = self.get_final_ids()
        for proc in self.get('decay_chains'):
            init = proc.get_initial_ids()[0]
            #while 1:
            try:
                pos = finals.index(init)
            except:
                break
            finals[pos] = proc.get_final_ids_after_decay()
        output = []
        for d in finals:
            if isinstance(d, list):
                output += d
            else:
                output.append(d)
        
        return output
    

    def get_final_legs(self):
        """Gives the final state legs"""

        return filter(lambda leg: leg.get('state') == True,
                       self.get('legs'))
    
    def get_final_ids(self):
        """Gives the pdg codes for final state particles"""

        return [l.get('id') for l in self.get_final_legs()]
    
                
    def get_legs_with_decays(self):
        """Return process with all decay chains substituted in."""

        if self['legs_with_decays']:
            return self['legs_with_decays']

        legs = copy.deepcopy(self.get('legs'))
        org_decay_chains = copy.copy(self.get('decay_chains'))
        sorted_decay_chains = []
        # Sort decay chains according to leg order
        for leg in legs:
            if not leg.get('state'): continue
            org_ids = [l.get('legs')[0].get('id') for l in \
                           org_decay_chains]
            if leg.get('id') in org_ids:
                sorted_decay_chains.append(org_decay_chains.pop(\
                                        org_ids.index(leg.get('id'))))
        assert not org_decay_chains
        ileg = 0
        for decay in sorted_decay_chains:
            while legs[ileg].get('state') == False or \
                      legs[ileg].get('id') != decay.get('legs')[0].get('id'):
                ileg = ileg + 1
            decay_legs = decay.get_legs_with_decays()
            legs = legs[:ileg] + decay_legs[1:] + legs[ileg+1:]
            ileg = ileg + len(decay_legs) - 1

        # Replace legs with copies
        legs = [copy.copy(l) for l in legs]

        for ileg, leg in enumerate(legs):
            leg.set('number', ileg + 1)
            
        self['legs_with_decays'] = LegList(legs)

        return self['legs_with_decays']
    
    def get_tag(self):
        """return the tag for standalone call"""
        
        initial = []    #filled in the next line
        final = [l.get('id') for l in self.get('legs')\
              if l.get('state') or initial.append(l.get('id'))]
        decay_finals = self.get_final_ids_after_decay()
        decay_finals.sort()
        tag = (tuple(initial), tuple(decay_finals))
        return tag
        

    def list_for_sort(self):
        """Output a list that can be compared to other processes as:
        [id, sorted(initial leg ids), sorted(final leg ids),
        sorted(decay list_for_sorts)]"""

        sorted_list =  [self.get('id'),
                        sorted(self.get_initial_ids()),
                        sorted(self.get_final_ids())]
        
        if self.get('decay_chains'):
            sorted_list.extend(sorted([d.list_for_sort() for d in \
                                       self.get('decay_chains')]))

        return sorted_list

    def compare_for_sort(self, other):
        """Sorting routine which allows to sort processes for
        comparison. Compare only process id and legs."""

        if self.list_for_sort() > other.list_for_sort():
            return 1
        if self.list_for_sort() < other.list_for_sort():
            return -1
        return 0
        
    def identical_particle_factor(self):
        """Calculate the denominator factor for identical final state particles
        """

        final_legs = filter(lambda leg: leg.get('state') == True, \
                              self.get_legs_with_decays())

        identical_indices = {}
        for leg in final_legs:
            if leg.get('id') in identical_indices:
                identical_indices[leg.get('id')] = \
                                    identical_indices[leg.get('id')] + 1
            else:
                identical_indices[leg.get('id')] = 1
        return reduce(lambda x, y: x * y, [ math.factorial(val) for val in \
                        identical_indices.values() ], 1)

    def check_expansion_orders(self):
        """Ensure that maximum expansion orders from the model are
        properly taken into account in the process"""

        # Ensure that expansion orders are taken into account
        expansion_orders = self.get('model').get('expansion_order')
        orders = self.get('orders')
        sq_orders = self.get('squared_orders')
        
        tmp = [(k,v) for (k,v) in expansion_orders.items() if 0 < v < 99]
        for (k,v) in tmp:  
            if k in orders:
                if v < orders[k]:
                    if k in sq_orders.keys() and \
                                             (sq_orders[k]>v or sq_orders[k]<0):
                        logger.warning(
'''The process with the squared coupling order (%s^2%s%s) specified can potentially 
recieve contributions with powers of the coupling %s larger than the maximal 
value allowed by the model builder (%s). Hence, MG5_aMC sets the amplitude order
for that coupling to be this maximal one. '''%(k,self.get('sqorders_types')[k],
                                             self.get('squared_orders')[k],k,v))
                    else:
                        logger.warning(
'''The coupling order (%s=%s) specified is larger than the one allowed 
             by the model builder. The maximal value allowed is %s. 
             We set the %s order to this value''' % (k,orders[k],v,k))
                    orders[k] = v
            else:
                orders[k] = v

    def __eq__(self, other):
        """Overloading the equality operator, so that only comparison
        of process id and legs is being done, using compare_for_sort."""

        if not isinstance(other, Process):
            return False

        return self.compare_for_sort(other) == 0

    def __ne__(self, other):
        return not self.__eq__(other)

#===============================================================================
# ProcessList
#===============================================================================
class ProcessList(PhysicsObjectList):
    """List of Process objects
    """

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

        return isinstance(obj, Process)

    def nice_string(self, indent = 0):
        """Returns a nicely formatted string of the matrix element processes."""

        mystr = "\n".join([p.nice_string(indent) for p in self])

        return mystr

#===============================================================================
# ProcessDefinition
#===============================================================================
class ProcessDefinition(Process):
    """ProcessDefinition: list of multilegs (ordered)
                          dictionary of orders
                          model
                          process id
    """

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

        super(ProcessDefinition, self).default_setup()

        self['legs'] = MultiLegList()
        # Decay chain processes associated with this process
        self['decay_chains'] = ProcessDefinitionList()
        if 'legs_with_decays' in self: del self['legs_with_decays']

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

        if name == 'legs':
            if not isinstance(value, MultiLegList):
                raise self.PhysicsObjectError, \
                        "%s is not a valid MultiLegList object" % str(value)
        elif name == 'decay_chains':
            if not isinstance(value, ProcessDefinitionList):
                raise self.PhysicsObjectError, \
                        "%s is not a valid ProcessDefinitionList" % str(value)

        else:
            return super(ProcessDefinition, self).filter(name, value)

        return True

    def has_multiparticle_label(self):
        """ Check that this process definition will yield a single process, as
        each multileg only has one leg"""
        
        for process in self['decay_chains']:
            if process.has_multiparticle_label():
                return True
            
        for mleg in self['legs']:
            if len(mleg['ids'])>1:
                return True
        
        return False

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

        keys = super(ProcessDefinition, self).get_sorted_keys()
        keys.remove('legs_with_decays')                                  

        return keys

    def get_minimum_WEIGHTED(self):
        """Retrieve the minimum starting guess for WEIGHTED order, to
        use in find_optimal_process_orders in MultiProcess diagram
        generation (as well as particles and hierarchy). The algorithm:

        1) Pick out the legs in the multiprocess according to the
        highest hierarchy represented (so don't mix particles from
        different hierarchy classes in the same multiparticles!)

        2) Find the starting maximum WEIGHTED order as the sum of the
        highest n-2 weighted orders

        3) Pick out required s-channel particle hierarchies, and use
        the highest of the maximum WEIGHTED order from the legs and
        the minimum WEIGHTED order extracted from 2*s-channel
        hierarchys plus the n-2-2*(number of s-channels) lowest
        leg weighted orders.
        """

        model = self.get('model')
        
        # Extract hierarchy and particles corresponding to the
        # different hierarchy levels from the model
        particles, hierarchy = model.get_particles_hierarchy()
        # Find legs corresponding to the different orders
        # making sure we look at lowest hierarchy first for each leg
        max_order_now = []
        new_legs =  copy.copy(self.get('legs'))
        import madgraph.core.base_objects as base_objects
        for parts, value in zip(particles, hierarchy):
            ileg = 0
            while ileg < len(new_legs):
                if any([id in parts for id in new_legs[ileg].get('ids')]):
                    max_order_now.append(value)
                    new_legs.pop(ileg)
                else:
                    ileg += 1

        # Now remove the two lowest orders to get maximum (since the
        # number of interactions is n-2)
        max_order_now = sorted(max_order_now)[2:]

        # Find s-channel propagators corresponding to the different orders
        max_order_prop = []
        for idlist in self.get('required_s_channels'):
            max_order_prop.append([0,0])
            for id in idlist:
                for parts, value in zip(particles, hierarchy):
                    if id in parts:
                        max_order_prop[-1][0] += 2*value
                        max_order_prop[-1][1] += 1
                        break

        if max_order_prop:
            if len(max_order_prop) >1:
                max_order_prop = min(*max_order_prop, key=lambda x:x[0])
            else:
                max_order_prop = max_order_prop[0]

            # Use either the max_order from the external legs or
            # the maximum order from the s-channel propagators, plus
            # the appropriate lowest orders from max_order_now
            max_order_now = max(sum(max_order_now),
                                max_order_prop[0] + \
                                sum(max_order_now[:-2 * max_order_prop[1]]))
        else:
            max_order_now = sum(max_order_now)            

        return max_order_now, particles, hierarchy

    def __iter__(self):
        """basic way to loop over all the process definition. 
           not used by MG which used some smarter version (use by ML)"""
        
        isids = [leg['ids'] for leg in self['legs'] \
                 if leg['state'] == False]
        fsids = [leg['ids'] for leg in self['legs'] \
                 if leg['state'] == True]

        red_isidlist = []
        # Generate all combinations for the initial state
        for prod in itertools.product(*isids):
            islegs = [Leg({'id':id, 'state': False}) for id in prod]  
            if tuple(sorted(prod)) in red_isidlist:
                    continue        
            red_isidlist.append(tuple(sorted(prod)))      
            red_fsidlist = []
            for prod in itertools.product(*fsids):
                # Remove double counting between final states
                if tuple(sorted(prod)) in red_fsidlist:
                    continue        
                red_fsidlist.append(tuple(sorted(prod)))
                leg_list = [copy.copy(leg) for leg in islegs]
                leg_list.extend([Leg({'id':id, 'state': True}) for id in prod])
                legs = LegList(leg_list)
                process = self.get_process_with_legs(legs)
                yield process

    def nice_string(self, indent=0, print_weighted=False, prefix=True):
        """Returns a nicely formated string about current process
        content"""

        if prefix:
            mystr = " " * indent + "Process: "
        else:
            mystr=""
        prevleg = None
        for leg in self['legs']:
            myparts = \
                   "/".join([self['model'].get('particle_dict')[id].get_name() \
                             for id in leg.get('ids')])
            if prevleg and prevleg['state'] == False \
                   and leg['state'] == True:
                # Separate initial and final legs by ">"
                mystr = mystr + '> '
                # Add required s-channels
                if self['required_s_channels'] and \
                       self['required_s_channels'][0]:
                    mystr += "|".join([" ".join([self['model'].\
                                       get('particle_dict')[req_id].get_name() \
                                                for req_id in id_list]) \
                                    for id_list in self['required_s_channels']])
                    mystr = mystr + '> '

            mystr = mystr + myparts + ' '
            #mystr = mystr + '(%i) ' % leg['number']
            prevleg = leg

        # Add forbidden s-channels
        if self['forbidden_onsh_s_channels']:
            mystr = mystr + '$ '
            for forb_id in self['forbidden_onsh_s_channels']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        # Add double forbidden s-channels
        if self['forbidden_s_channels']:
            mystr = mystr + '$$ '
            for forb_id in self['forbidden_s_channels']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        # Add forbidden particles
        if self['forbidden_particles']:
            mystr = mystr + '/ '
            for forb_id in self['forbidden_particles']:
                forbpart = self['model'].get('particle_dict')[forb_id]
                mystr = mystr + forbpart.get_name() + ' '

        if self['orders']:
            mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
                       for key in sorted(self['orders'])]) + ' '

        if self['constrained_orders']:
            mystr = mystr + " ".join('%s%s%d' % (key, operator, value) for 
                                     (key,(value, operator)) 
                                   in self['constrained_orders'].items()) + ' '

        # Add perturbation_couplings
        if self['perturbation_couplings']:
            mystr = mystr + '[ '
            if self['NLO_mode']!='tree':
                if self['NLO_mode']=='virt' and not self['has_born']:
                    mystr = mystr + 'sqrvirt = '
                else:
                    mystr = mystr + self['NLO_mode'] + ' = '
            for order in self['perturbation_couplings']:
                mystr = mystr + order + ' '
            mystr = mystr + '] '

        if self['squared_orders']:
            mystr = mystr + " ".join([key + '^2%s%d'%\
                (self.get_squared_order_type(key),self['squared_orders'][key]) \
              for key in self['squared_orders'].keys() \
                                    if print_weighted or key!='WEIGHTED']) + ' '

        # Remove last space
        mystr = mystr[:-1]

        if self.get('id') or self.get('overall_orders'):
            mystr += " @%d" % self.get('id')
            if self.get('overall_orders'):
                mystr += " " + " ".join([key + '=' + repr(self['orders'][key]) \
                       for key in sorted(self['orders'])]) + ' '
        
        if not self.get('decay_chains'):
            return mystr

        for decay in self['decay_chains']:
            mystr = mystr + '\n' + \
                    decay.nice_string(indent + 2).replace('Process', 'Decay')

        return mystr

    def get_process_with_legs(self, LegList):
        """ Return a Process object which has the same properties of this 
            ProcessDefinition but with the specified LegList as legs attribute. 
            """
            
        return Process({\
            'legs': LegList,
            'model':self.get('model'),
            'id': self.get('id'),
            'orders': self.get('orders'),
            'sqorders_types': self.get('sqorders_types'),
            'squared_orders': self.get('squared_orders'),
            'constrained_orders': self.get('constrained_orders'),
            'has_born': self.get('has_born'),
            'required_s_channels': self.get('required_s_channels'),
            'forbidden_onsh_s_channels': self.get('forbidden_onsh_s_channels'),            
            'forbidden_s_channels': self.get('forbidden_s_channels'),
            'forbidden_particles': self.get('forbidden_particles'),
            'perturbation_couplings': self.get('perturbation_couplings'),
            'is_decay_chain': self.get('is_decay_chain'),
            'overall_orders': self.get('overall_orders'),
            'split_orders': self.get('split_orders'),
            'born_orders': self.get('born_orders'),
            'NLO_mode': self.get('NLO_mode')
            })
            
    def get_process(self, initial_state_ids, final_state_ids):
        """ Return a Process object which has the same properties of this 
            ProcessDefinition but with the specified given leg ids. """
        
        # First make sure that the desired particle ids belong to those defined
        # in this process definition.
        my_isids = [leg.get('ids') for leg in self.get('legs') \
              if not leg.get('state')]
        my_fsids = [leg.get('ids') for leg in self.get('legs') \
             if leg.get('state')]            
        for i, is_id in enumerate(initial_state_ids):
            assert is_id in my_isids[i]
        for i, fs_id in enumerate(final_state_ids):
            assert fs_id in my_fsids[i]
        
        return self.get_process_with_legs(LegList(\
               [Leg({'id': id, 'state':False}) for id in initial_state_ids] + \
               [Leg({'id': id, 'state':True}) for id in final_state_ids]))

    def __eq__(self, other):
        """Overloading the equality operator, so that only comparison
        of process id and legs is being done, using compare_for_sort."""

        return super(Process, self).__eq__(other)

#===============================================================================
# ProcessDefinitionList
#===============================================================================
class ProcessDefinitionList(PhysicsObjectList):
    """List of ProcessDefinition objects
    """

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

        return isinstance(obj, ProcessDefinition)

#===============================================================================
# Global helper functions
#===============================================================================

def make_unique(doubletlist):
    """Make sure there are no doublets in the list doubletlist.
    Note that this is a slow implementation, so don't use if speed 
    is needed"""

    assert isinstance(doubletlist, list), \
           "Argument to make_unique must be list"
    

    uniquelist = []
    for elem in doubletlist:
        if elem not in uniquelist:
            uniquelist.append(elem)

    doubletlist[:] = uniquelist[:]