No OneTemporary
Actions

Size

114 KB

Subscribers

None

View Options

	diff --git a/Models/Feynrules/python/ufo2herwig b/Models/Feynrules/python/ufo2herwig
	--- a/Models/Feynrules/python/ufo2herwig
	+++ b/Models/Feynrules/python/ufo2herwig
	@@ -1,387 +1,392 @@
	#! /usr/bin/env python
	from __future__ import division
	import os, sys, pprint, argparse, re,copy
	from string import strip, Template

	# add path to the ufo conversion modules
	modulepath = os.path.join("@PKGLIBDIR@",'python')
	sys.path.append(modulepath)

	from ufo2peg import *

	# set up the option parser for command line input
	parser = argparse.ArgumentParser(
	description='Create Herwig model files from Feynrules UFO input.'
	)
	parser.add_argument(
	'ufodir',
	metavar='UFO_directory',
	help='the UFO model directory'
	)
	parser.add_argument(
	'-v', '--verbose',
	action="store_true",
	help="print verbose output"
	)
	parser.add_argument(
	'-n','--name',
	default="FRModel",
	help="set custom nametag for the model"
	)
	parser.add_argument(
	'--ignore-skipped',
	action="store_true",
	help="ignore skipped vertices and produce output anyway"
	)
	parser.add_argument(
	'--split-model',
	action="store_true",
	help="Split the model file into pieces to improve compilation for models with many parameters"
	)
	parser.add_argument(
	'--no-generic-loop-vertices',
	action="store_true",
	help="Don't include the automatically generated generic loop vertices for h->gg and h->gamma gamma"
	)
	parser.add_argument(
	+ '--include-generic',
	+ action="store_true",
	+ help="Include support for generic spin structures (still experimental)"
	+)
	+parser.add_argument(
	'--forbidden-particle-name',
	action="append",
	default=["eta","phi"],
	help="Add particle names not allowed as names in UFO models to avoid conflicts with"+\
	"Herwig internal particles, names will have _UFO appended"
	)
	# get the arguments
	args = parser.parse_args()
	# import the model
	import imp
	path,mod = os.path.split(os.path.abspath(args.ufodir))
	FR = imp.load_module(mod,*imp.find_module(mod,[path]))

	##################################################
	##################################################

	# get the Model name from the arguments
	modelname = args.name
	libname = modelname + '.so'



	# define arrays and variables
	#allplist = ""
	parmdecls = []
	parmgetters = []
	parmconstr = []
	doinit = []

	paramstoreplace_ = []
	paramstoreplace_expressions_ = []

	# get external parameters for printing
	parmsubs = dict( [ (p.name, p.value)
	for p in FR.all_parameters
	if p.nature == 'external' ] )

	# evaluate python cmath
	def evaluate(x):
	import cmath
	return eval(x,
	{'cmath':cmath,
	'complexconjugate':FR.function_library.complexconjugate},
	parmsubs)


	## get internal params into arrays
	internal = ( p for p in FR.all_parameters
	if p.nature == 'internal' )

	#paramstoreplaceEW_ = []
	#paramstoreplaceEW_expressions_ = []

	# calculate internal parameters
	for p in internal:
	parmsubs.update( { p.name : evaluate(p.value) } )
	# if 'aS' in p.value and p.name != 'aS':
	# paramstoreplace_.append(p.name)
	# paramstoreplace_expressions_.append(p.value)
	# if 'aEWM1' in p.value and p.name != 'aEWM1':
	# paramstoreplaceEW_.append(p.name)
	# paramstoreplaceEW_expressions_.append(p.value)
	parmvalues=copy.copy(parmsubs)


	# more arrays used for substitution in templates
	paramsforstream = []
	parmmodelconstr = []

	# loop over parameters and fill in template stuff according to internal/external and complex/real
	# WARNING: Complex external parameter input not tested!
	if args.verbose:
	print 'verbose mode on: printing all parameters'
	print '-'*60
	paramsstuff = ('name', 'expression', 'default value', 'nature')
	pprint.pprint(paramsstuff)



	interfacedecl_T = """\
	static Parameter<{modelname}, {type}> interface{pname}
	("{pname}",
	"The interface for parameter {pname}",
	&{modelname}::{pname}_, {value}, 0, 0,
	false, false, Interface::nolimits);
	"""

	massnames = {}
	widthnames = {}

	for particle in FR.all_particles:
	# skip ghosts and goldstones
	if(isGhost(particle) or isGoldstone(particle)) :
	continue
	if particle.mass != 'ZERO':
	massnames[particle.mass] = abs(particle.pdg_code)

	if particle.width != 'ZERO':
	widthnames[particle.width] = abs(particle.pdg_code)

	interfaceDecls = []
	modelparameters = {}

	for p in FR.all_parameters:
	value = parmsubs[p.name]
	if p.type == 'real':
	assert( value.imag < 1.0e-16 )
	value = value.real
	if p.nature == 'external':
	if p not in massnames and p not in widthnames:
	interfaceDecls.append(
	interfacedecl_T.format(modelname=modelname,
	pname=p.name,
	value=value,
	type=typemap(p.type))
	)
	else:
	interfaceDecls.append('\n// no interface for %s. Use particle definition instead.\n' % p.name)
	if hasattr(p,'lhablock'):
	lhalabel = '{lhablock}_{lhacode}'.format( lhablock=p.lhablock.upper(), lhacode='_'.join(map(str,p.lhacode)) )
	if p not in massnames and p not in widthnames:
	parmmodelconstr.append('set %s:%s ${%s}' % (modelname, p.name, lhalabel))
	else:
	parmmodelconstr.append('# %s is taken from the particle setup' % p.name)
	modelparameters[lhalabel] = value
	parmsubs[p.name] = lhalabel
	else:
	if p not in massnames and p not in widthnames:
	parmmodelconstr.append('set %s:%s %s' % (modelname, p.name, value))
	else:
	parmmodelconstr.append('# %s is taken from the particle setup' % p.name)
	parmsubs[p.name] = value
	if p not in massnames and p not in widthnames:
	parmconstr.append('%s_(%s)' % (p.name, value))
	else:
	parmconstr.append('%s_()' % p.name)
	else :
	parmconstr.append('%s_()' % p.name)
	parmsubs[p.name] = value
	elif p.type == 'complex':
	value = complex(value)
	if p.nature == 'external':
	#
	# TODO: WE DO NOT HAVE COMPLEX INTERFACES IN THEPEG (yet?)
	#
	# interfaceDecls.append(
	# interfacedecl_T.format(modelname=modelname,
	# pname=p.name,
	# value='Complex(%s,%s)'%(value.real,value.imag),
	# type=typemap(p.type))
	# )
	#
	# parmmodelconstr.append('set %s:%s (%s,%s)' % (modelname, p.name, value.real, value.imag))
	parmconstr.append('%s_(%s,%s)' % (p.name, value.real, value.imag))
	else :
	parmconstr.append('%s_(%s,%s)' % (p.name, 0.,0.))
	parmsubs[p.name] = value
	else:
	raise Exception('Unknown data type "%s".' % p.type)

	parmdecls.append(' %s %s_;' % (typemap(p.type), p.name))
	parmgetters.append(' %s %s() const { return %s_; }' % (typemap(p.type),p.name, p.name))
	paramsforstream.append('%s_' % p.name)
	expression, symbols = 'return %s_' % p.name, None
	if p.nature != 'external':
	expression, symbols = py2cpp(p.value)
	text = add_brackets(expression, symbols)
	text=text.replace('()()','()')
	doinit.append(' %s_ = %s;' % (p.name, text) )
	if p in massnames:
	doinit.append(' resetMass(%s,%s_ * GeV);' % (massnames[p], p.name) )
	elif p in widthnames:
	doinit.append(' getParticleData(%s)->width(%s_ * GeV);' % (widthnames[p], p.name) )
	doinit.append(' getParticleData(%s)->cTau (%s_ == 0.0 ? Length() : hbarc/(%s_*GeV));' % (widthnames[p], p.name, p.name) )
	doinit.append(' getParticleData(%s)->widthCut(10. * %s_ * GeV);' % (widthnames[p], p.name) )


	elif p.nature == 'external':
	if p in massnames:
	doinit.append(' %s_ = getParticleData(%s)->mass() / GeV;' % (p.name, massnames[p]) )
	elif p in widthnames:
	doinit.append(' %s_ = getParticleData(%s)->width() / GeV;' % (p.name, widthnames[p]) )

	if args.verbose:
	pprint.pprint((p.name,p.value, value, p.nature))

	pcwriter = ParamCardWriter(FR.all_parameters)
	paramcard_output = '\n'.join(pcwriter.output)

	### special treatment
	# if p.name == 'aS':
	# expression = '0.25 * sqr(strongCoupling(q2)) / Constants::pi'
	# elif p.name == 'aEWM1':
	# expression = '4.0 * Constants::pi / sqr(electroMagneticCoupling(q2))'
	# elif p.name == 'Gf':
	# expression = 'generator()->standardModel()->fermiConstant() * GeV2'

	paramconstructor=': '

	for ncount in range(0,len(parmconstr)) :
	paramconstructor += parmconstr[ncount]
	if(ncount != len(parmconstr) -1) :
	paramconstructor += ','
	if(ncount%5 == 0 ) :
	paramconstructor += "\n"

	paramout=""
	paramin =""
	for ncount in range(0,len(paramsforstream)) :
	if(ncount !=0 ) :
	paramout += "<< " + paramsforstream[ncount]
	paramin += ">> " + paramsforstream[ncount]
	else :
	paramout += paramsforstream[ncount]
	paramin += paramsforstream[ncount]
	if(ncount%5 == 0 ) :
	paramout += "\n"
	paramin += "\n"

	parmtextsubs = { 'parmgetters' : '\n'.join(parmgetters),
	'parmdecls' : '\n'.join(parmdecls),
	'parmconstr' : paramconstructor,
	'getters' : '',
	'decls' : '',
	'addVertex' : '',
	'doinit' : '\n'.join(doinit),
	'ostream' : paramout,
	'istream' : paramin ,
	'refs' : '',
	'parmextinter': ''.join(interfaceDecls),
	'ModelName': modelname,
	'calcfunctions': '',
	'param_card_data': paramcard_output
	}

	##################################################
	##################################################
	##################################################




	# set up the conversion of the vertices
	vertexConverter = VertexConverter(FR,parmvalues)
	vertexConverter.readArgs(args)
	# convert the vertices
	vertexConverter.convert()

	##################################################
	##################################################
	##################################################



	plist, names = thepeg_particles(FR,parmsubs,modelname,modelparameters,args.forbidden_particle_name)

	particlelist = [
	"# insert HPConstructor:Outgoing 0 /Herwig/{n}/Particles/{p}".format(n=modelname,p=p)
	for p in names
	]
	# make the first one active to have something runnable in the example .in file
	particlelist[0] = particlelist[0][2:]
	particlelist = '\n'.join(particlelist)

	modelfilesubs = { 'plist' : plist,
	'vlist' : vertexConverter.get_vertices(libname),
	'setcouplings': '\n'.join(parmmodelconstr),
	'ModelName': modelname
	}

	# write the files from templates according to the above subs
	if vertexConverter.should_print():
	MODEL_HWIN = getTemplate('LHC-FR.in')
	if(not args.split_model) :
	MODEL_CC = [getTemplate('Model.cc')]
	else :
	MODEL_EXTRA_CC=getTemplate('Model6.cc')
	extra_names=[]
	extra_calcs=[]
	parmtextsubs['doinit']=""
	for i in range(0,len(doinit)) :
	if( i%20 == 0 ) :
	function_name = "initCalc" +str(int(i/20))
	parmtextsubs['doinit'] += function_name +"();\n"
	parmtextsubs['calcfunctions'] += "void " + function_name + "();\n"
	extra_names.append(function_name)
	extra_calcs.append("")
	extra_calcs[-1] += doinit[i] + "\n"
	for i in range(0,len(extra_names)) :
	ccname = '%s_extra_%s.cc' % (modelname,i)
	writeFile( ccname, MODEL_EXTRA_CC.substitute({'ModelName' : modelname,
	'functionName' : extra_names[i],
	'functionCalcs' : extra_calcs[i] }) )




	MODEL_CC = [getTemplate('Model1.cc'),getTemplate('Model2.cc'),getTemplate('Model3.cc'),
	getTemplate('Model4.cc'),getTemplate('Model5.cc')]
	MODEL_H = getTemplate('Model.h')
	print 'LENGTH',len(MODEL_CC)
	MODELINFILE = getTemplate('FR.model')
	writeFile( 'LHC-%s.in' % modelname,
	MODEL_HWIN.substitute({ 'ModelName' : modelname,
	'Particles' : particlelist })
	)

	modeltemplate = Template( MODELINFILE.substitute(modelfilesubs) )

	writeFile( '%s.h' % modelname, MODEL_H.substitute(parmtextsubs) )
	for i in range(0,len(MODEL_CC)) :
	if(len(MODEL_CC)==1) :
	ccname = '%s.cc' % modelname
	else :
	ccname = '%s.cc' % (modelname + str(i))
	writeFile( ccname, MODEL_CC[i].substitute(parmtextsubs) )
	writeFile( modelname +'.template', modeltemplate.template )
	writeFile( modelname +'.model', modeltemplate.substitute( modelparameters ) )
	# copy the Makefile-FR to current directory,
	# replace with the modelname for compilation
	with open(os.path.join(modulepath,'Makefile-FR'),'r') as orig:
	with open('Makefile','w') as dest:
	dest.write(orig.read().replace("FeynrulesModel.so", libname))


	print 'finished generating model:\t', modelname
	print 'model directory:\t\t', args.ufodir
	print 'generated:\t\t\t', len(FR.all_vertices), 'vertices'
	print '='*60
	print 'library:\t\t\t', libname
	print 'input file:\t\t\t', 'LHC-' + modelname +'.in'
	print 'model file:\t\t\t', modelname +'.model'
	print '='*60
	print """\
	To complete the installation, compile by typing "make".
	An example input file is provided as LHC-FRModel.in,
	you'll need to change the required particles in there.
	"""
	print 'DONE!'
	print '='*60

	diff --git a/Models/Feynrules/python/ufo2peg/Vertex_class.template b/Models/Feynrules/python/ufo2peg/GeneralVertex_class.template
	copy from Models/Feynrules/python/ufo2peg/Vertex_class.template
	copy to Models/Feynrules/python/ufo2peg/GeneralVertex_class.template
	--- a/Models/Feynrules/python/ufo2peg/Vertex_class.template
	+++ b/Models/Feynrules/python/ufo2peg/GeneralVertex_class.template
	@@ -1,43 +1,39 @@
	-class ${ModelName}${classname}: public ${lorentztag}Vertex {
	+class ${ModelName}${classname}: public Abstract${lorentztag}Vertex {
	public:
	${ModelName}${classname}() {
	- ${header}
	${addToPlist}
	}
	- void setCoupling(Energy2 ${couplingptrs}) {
	- ${symbolrefs}
	- ${prepend}
	- // getParams(q2);
	- ${norm}
	- ${left}
	- ${right}
	- ${append}
	- }
	+
	+ ${evaldefs}
	+
	void persistentOutput(PersistentOStream & os) const { os << model_; }
	void persistentInput(PersistentIStream & is, int) { is >> model_; }
	- // static void Init();
	+
	+ virtual void setCoupling(Energy2, tcPDPtr,
	+ tcPDPtr, tcPDPtr) {assert(false);}
	+
	+ virtual void setCoupling(Energy2,tcPDPtr,tcPDPtr,tcPDPtr,
	+ tcPDPtr) {assert(false);}
	+
	protected:
	+
	IBPtr clone() const { return new_ptr(*this); }
	IBPtr fullclone() const { return new_ptr(*this); }
	void doinit() {
	model_ = dynamic_ptr_cast<tcHw${ModelName}Ptr>
	(generator()->standardModel());
	assert(model_);
	- // getParams(q2);
	- ${parameters}
	- ${setCouplings}
	- orderInGem(${qedorder});
	- orderInGs(${qcdorder});
	- ${lorentztag}Vertex::doinit();
	+ Abstract${lorentztag}Vertex::doinit();
	}
	- // void getParams(Energy2);
	+
	private:
	+
	${ModelName}${classname} & operator=(const ${ModelName}${classname} &);
	- // Complex leftval, rightval, normval;
	+
	tcHw${ModelName}Ptr model_;
	};
	-DescribeClass<${ModelName}${classname},Helicity::${lorentztag}Vertex>
	+DescribeClass<${ModelName}${classname},Helicity::Abstract${lorentztag}Vertex>
	describeHerwig${ModelName}${classname}("Herwig::${ModelName}${classname}",
	"${ModelName}.so");
	-// void ${ModelName}${classname}::getParams(Energy2 ) {
	-// }
	+
	+${evalimpls}
	diff --git a/Models/Feynrules/python/ufo2peg/check_lorentz.py b/Models/Feynrules/python/ufo2peg/check_lorentz.py
	--- a/Models/Feynrules/python/ufo2peg/check_lorentz.py
	+++ b/Models/Feynrules/python/ufo2peg/check_lorentz.py
	@@ -1,864 +1,1435 @@
	import itertools,cmath,re,sys
	-from .helpers import SkipThisVertex,extractAntiSymmetricIndices
	+from .helpers import SkipThisVertex,extractAntiSymmetricIndices,def_from_model
	from .converter import py2cpp
	from .lorentzparser import parse_lorentz
	+import sympy,string
	+from string import Template
	+from sympy import Matrix,Symbol

	def compare(a,b) :
	num=abs(a-b)
	den=abs(a+b)
	if(den == 0. and 1e-10) :
	return True
	return num/den<1e-10

	def evaluate(x,model,parmsubs):
	import cmath
	return eval(x,
	{'cmath':cmath,
	'complexconjugate':model.function_library.complexconjugate},
	parmsubs)

	# ordering for EW VVV vertices
	def VVVordering(vertex) :
	pattern = "if((p1->id()==%s&&p2->id()==%s&&p3->id()==%s)"+\
	"\|\|(p1->id()==%s&&p2->id()==%s&&p3->id()==%s)\|\|"+\
	"(p1->id()==%s&&p2->id()==%s&&p3->id()==%s)) {norm(-norm());}"
	ordering = pattern%(vertex.particles[1].pdg_code,
	vertex.particles[0].pdg_code,
	vertex.particles[2].pdg_code,
	vertex.particles[0].pdg_code,
	vertex.particles[2].pdg_code,
	vertex.particles[1].pdg_code,
	vertex.particles[2].pdg_code,
	vertex.particles[1].pdg_code,
	vertex.particles[0].pdg_code)
	return ordering

	def tensorCouplings(vertex,value,prefactors,L,lorentztag,pos,all_couplings,order) :
	# split the structure into its different terms for analysis
	ordering=""
	structures = extractStructures(L)
	if(lorentztag == 'SST') :
	terms=[['P(1003,2)','P(2003,1)'],
	['P(1003,1)','P(2003,2)'],
	['P(-1,1)','P(-1,2)','Metric(1003,2003)'],
	['Metric(1003,2003)']]
	signs=[1.,1.,-1.,-1.]
	new_couplings=[False]*len(terms)
	elif(lorentztag == 'FFT' ) :
	terms=[['P(2003,1)','Gamma(1003,2,1)'],
	['P(2003,2)','Gamma(1003,2,1)'],
	['P(1003,1)','Gamma(2003,2,1)'],
	['P(1003,2)','Gamma(2003,2,1)'],
	['P(-1,1)','Gamma(-1,2,1)','Metric(1003,2003)'],
	['P(-1,2)','Gamma(-1,2,1)','Metric(1003,2003)'],
	['Metric(1003,2003)']]
	signs=[1.,-1.,1.,-1.,-0.5,0.5,1.]
	new_couplings=[False]3len(terms)
	elif(lorentztag == 'VVT' ) :
	terms=[['P(-1,1)','P(-1,2)','Metric(1,2003)','Metric(2,1003)'], # from C term
	['P(-1,1)','P(-1,2)','Metric(1,1003)','Metric(2,2003)'], # from C term
	['P(-1,1)','P(-1,2)','Metric(1,2)','Metric(1003,2003)'], # from C term
	['P(1,2)','P(2,1)','Metric(1003,2003)'], # from D term (sym)
	['P(1,2)','P(2003,1)','Metric(2,1003)'], # 1st term
	['P(1,2)','P(1003,1)','Metric(2,2003)'], # 1st swap
	['P(2,1)','P(2003,2)','Metric(1,1003)'], # 2nd term
	['P(2,1)','P(1003,2)','Metric(1,2003)'], # 2nd swap
	['P(1003,2)','P(2003,1)','Metric(1,2)'], # 3rd term
	['P(1003,1)','P(2003,2)','Metric(1,2)'], # 3rd swap
	['Metric(1,2003)','Metric(2,1003)'], # from mass term
	['Metric(1,1003)','Metric(2,2003)'], # from mass term
	['Metric(1,2)','Metric(1003,2003)'], # from mass term
	['P(1,1)','P(2,1)','Metric(1003,2003)'], # gauge terms
	['P(1,2)','P(2,2)','Metric(1003,2003)'], # gauge terms
	['P(1,1)','P(2,2)','Metric(1003,2003)'], # gauge terms
	['P(1003,1)','P(1,1)','Metric(2,2003)'], # gauge terms
	['P(1003,2)','P(2,2)','Metric(1,2003)'], # gauge terms
	['P(2003,1)','P(1,1)','Metric(2,1003)'], # gauge terms
	['P(2003,2)','P(2,2)','Metric(1,1003)'], # gauge terms
	]
	signs=[1.,1.,-1.,1.,-1.,-1.,-1.,-1.,1.,1.,1.,1.,-1.,1.,1.,0.25,-1.,-1.,-1.,-1.]
	new_couplings=[False]*len(terms)
	elif(lorentztag == 'FFVT' ) :
	terms = [['Gamma(2004,2,1)','Metric(3,1004)'],
	['Gamma(1004,2,1)','Metric(3,2004)'],
	['Gamma(3,2,1)','Metric(1004,2004)'],
	['Gamma(2004,2,-1)','Metric(3,1004)'],
	['Gamma(1004,2,-1)','Metric(3,2004)'],
	['Gamma(3,2,-1)','Metric(1004,2004)']]
	signs=[1.,1.,-0.5,1.,1.,-0.5]
	new_couplings=[False]3len(terms)
	elif(lorentztag == 'VVVT' ) :
	# the F(mu nu,rho sigma lambda) terms first
	terms = [['P(2004,2)','Metric(1,1004)','Metric(2,3)'],['P(2004,3)','Metric(1,1004)','Metric(2,3)'],
	['P(1004,2)','Metric(1,2004)','Metric(2,3)'],['P(1004,3)','Metric(1,2004)','Metric(2,3)'],
	['P(2004,3)','Metric(1,3)','Metric(2,1004)'],['P(2004,1)','Metric(1,3)','Metric(2,1004)'],
	['P(1004,3)','Metric(1,3)','Metric(2,2004)'],['P(1004,1)','Metric(1,3)','Metric(2,2004)'],
	['P(2004,1)','Metric(1,2)','Metric(3,1004)'],['P(2004,2)','Metric(1,2)','Metric(3,1004)'],
	['P(1004,1)','Metric(1,2)','Metric(3,2004)'],['P(1004,2)','Metric(1,2)','Metric(3,2004)'],
	['P(3,1)','Metric(1,2004)','Metric(2,1004)'],['P(3,2)','Metric(1,2004)','Metric(2,1004)'],
	['P(3,1)','Metric(1,1004)','Metric(2,2004)'],['P(3,2)','Metric(1,1004)','Metric(2,2004)'],
	['P(3,1)','Metric(1,2)','Metric(1004,2004)'],['P(3,2)','Metric(1,2)','Metric(1004,2004)'],
	['P(2,3)','Metric(1,2004)','Metric(3,1004)'],['P(2,1)','Metric(1,2004)','Metric(3,1004)'],
	['P(2,3)','Metric(1,1004)','Metric(3,2004)'],['P(2,1)','Metric(1,1004)','Metric(3,2004)'],
	['P(2,3)','Metric(1,3)','Metric(1004,2004)'],['P(2,1)','Metric(1,3)','Metric(1004,2004)'],
	['P(1,2)','Metric(2,2004)','Metric(3,1004)'],['P(1,3)','Metric(2,2004)','Metric(3,1004)'],
	['P(1,2)','Metric(2,1004)','Metric(3,2004)'],['P(1,3)','Metric(2,1004)','Metric(3,2004)'],
	['P(1,2)','Metric(2,3)','Metric(1004,2004)'],['P(1,3)','Metric(2,3)','Metric(1004,2004)']]
	signs = [1.,-1.,1.,-1.,1.,-1.,1.,-1.,1.,-1.,1.,-1.,
	1.,-1.,1.,-1.,-1.,1.,1.,-1.,1.,-1.,-1.,1.,1.,-1.,1.,-1.,-1.,1.]
	new_couplings=[False]*len(terms)
	l = lambda c: len(pos[c])
	if l(8)!=3 :
	ordering = VVVordering(vertex)
	# unknown
	else :
	raise Exception('Unknown data type "%s".' % lorentztag)
	iterm=0
	try :
	for term in terms:
	for perm in itertools.permutations(term):
	label = '*'.join(perm)
	for istruct in range(0,len(structures)) :
	if label in structures[istruct] :
	reminder = structures[istruct].replace(label,'1.',1)
	loc=iterm
	if(reminder.find("ProjM")>=0) :
	reminder=re.sub("\ProjM$.,.$","",reminder)
	loc+=len(terms)
	elif(reminder.find("ProjP")>=0) :
	reminder=re.sub("\ProjP$.,.$","",reminder)
	loc+=2*len(terms)
	structures[istruct] = "Done"
	val = eval(reminder, {'cmath':cmath} )*signs[iterm]
	if(new_couplings[loc]) :
	new_couplings[loc] += val
	else :
	new_couplings[loc] = val
	iterm+=1
	except :
	SkipThisVertex()
	# check we've handled all the terms
	for val in structures:
	if(val!="Done") :
	raise SkipThisVertex()
	# special for FFVT
	if(lorentztag=="FFVT") :
	t_couplings=new_couplings
	new_couplings=[False]*9
	for i in range(0,9) :
	j = i+3*(i/3)
	k = i+3+3*(i/3)
	if( not t_couplings[j]) :
	new_couplings[i] = t_couplings[k]
	else :
	new_couplings[i] = t_couplings[j]
	# set the couplings
	for icoup in range(0,len(new_couplings)) :
	if(new_couplings[icoup]) :
	new_couplings[icoup] = '(%s) * (%s) *(%s)' % (new_couplings[icoup],prefactors,value)
	if(len(all_couplings)==0) :
	all_couplings=new_couplings
	else :
	for icoup in range(0,len(new_couplings)) :
	if(new_couplings[icoup] and all_couplings[icoup]) :
	all_couplings[icoup] = '(%s) + (%s) ' % (new_couplings[icoup],all_couplings[icoup])
	elif(new_couplings[icoup]) :
	all_couplings[icoup] = new_couplings[icoup]
	# return the results
	return (ordering,all_couplings)

	def processTensorCouplings(lorentztag,vertex,model,parmsubs,all_couplings) :
	# check for fermion vertices (i.e. has L/R couplings)
	fermions = "FF" in lorentztag
	# test and process the values of the couplings
	tval = [False]*3
	value = [False]*3
	# loop over the colours
	for icolor in range(0,len(all_couplings)) :
	lmax = len(all_couplings[icolor])
	if(fermions) : lmax /=3
	# loop over the different terms
	for ix in range(0,lmax) :
	test = [False]*3
	# normal case
	if( not fermions ) :
	test[0] = all_couplings[icolor][ix]
	else :
	# first case vector but no L/R couplings
	if( not all_couplings[icolor][lmax+ix] and
	not all_couplings[icolor][2*lmax+ix] ) :
	test[0] = all_couplings[icolor][ix]
	# special for mass terms and massless particles
	if(not all_couplings[icolor][ix]) :
	code = abs(vertex.particles[0].pdg_code)
	if(ix==6 and code ==12 or code ==14 or code==16) :
	continue
	else :
	raise SkipThisVertex()
	# second case L/R couplings
	elif( not all_couplings[icolor][ix] ) :
	# L=R, replace with vector
	if(all_couplings[icolor][lmax+ix] ==
	all_couplings[icolor][2*lmax+ix]) :
	test[0] = all_couplings[icolor][lmax+ix]
	else :
	test[1] = all_couplings[icolor][lmax+ix]
	test[2] = all_couplings[icolor][2*lmax+ix]
	else :
	raise SkipThisVertex()
	# special handling of mass terms
	# scalar divide by m**2
	if((ix==3 and lorentztag=="SST") or
	( ix>=10 and ix<=12 and lorentztag=="VVT" )) :
	for i in range(0,len(test)) :
	if(test[i]) :
	test[i] = '(%s)/%s**2' % (test[i],vertex.particles[0].mass.value)
	# fermions divide by 4*m
	elif(ix==6 and lorentztag=="FFT" and
	float(vertex.particles[0].mass.value) != 0. ) :
	for i in range(0,len(test)) :
	if(test[i]) :
	test[i] = '-(%s)/%s/4' % (test[i],vertex.particles[0].mass.value)
	# set values on first pass
	if(not tval[0] and not tval[1] and not tval[2]) :
	value = test
	for i in range(0,len(test)) :
	if(test[i]) : tval[i] = evaluate(test[i],model,parmsubs)
	else :
	for i in range(0,len(test)) :
	if(not test[i] and not tval[i]) :
	continue
	if(not test[i] or not tval[i]) :
	# special for mass terms and vectors
	if(lorentztag=="VVT" and ix >=10 and ix <=12 and
	float(vertex.particles[0].mass.value) == 0. ) :
	continue
	# special for vector gauge terms
	if(lorentztag=="VVT" and ix>=13) :
	continue
	raise SkipThisVertex()
	tval2 = evaluate(test[i],model,parmsubs)
	if(abs(tval[i]-tval2)>1e-6) :
	# special for fermion mass term if fermions massless
	if(lorentztag=="FFT" and ix ==6 and tval2 == 0. and
	float(vertex.particles[0].mass.value) == 0. ) :
	continue
	raise SkipThisVertex()
	# simple clean up
	for i in range(0,len(value)):
	if(value[i]) :
	value[i] = value[i].replace("(1.0) * ","").replace(" * (1)","")
	# put everything together
	coup_left = 0.
	coup_right = 0.
	coup_norm = 0.
	if(lorentztag == "SST" or lorentztag == "VVT" or
	lorentztag == "VVVT" or lorentztag == "FFT" ) :
	coup_norm = value[0]
	if(value[1] or value[2]) :
	raise SkipThisVertex()
	elif(lorentztag=="FFVT") :
	if(not value[1] and not value[2]) :
	coup_norm = value[0]
	coup_left = "1."
	coup_right = "1."
	elif(not value[0]) :
	coup_norm = "1."
	if(value[1] and value[2]) :
	coup_left = value[1]
	coup_right = value[2]
	elif(value[1]) :
	coup_left = value[1]
	coup_right = "0."
	elif(value[2]) :
	coup_left = "0."
	coup_right = value[2]
	else :
	raise SkipThisVertex()
	else :
	raise SkipThisVertex()
	else :
	raise SkipThisVertex()
	# return the answer
	return (coup_left,coup_right,coup_norm)

	def extractStructures(L) :
	structure1 = L.structure.split()
	structures =[]
	sign=''
	for struct in structure1 :
	if(struct=='+') :
	continue
	elif(struct=='-') :
	sign='-'
	else :
	structures.append(sign+struct.strip())
	sign=''
	return structures



	def changeSign(sign1,sign2) :
	if((sign1=="+" and sign2=="+") or\
	(sign1=="-" and sign2=="-")) :
	return "+"
	else :
	return "-"

	def epsilonOrder(eps) :
	terms,sign = extractAntiSymmetricIndices(eps,"Epsilon(")
	return (sign,"Epsilon(%s,%s,%s,%s)" % (terms[0],terms[1],terms[2],terms[3]))


	def VVSEpsilon(couplings,struct) :
	if(struct.find("Epsilon")<0) :
	return
	fact=""
	sign="+"
	if(struct[-1]==")") :
	fact=struct.split("(")[0]
	if(fact.find("Epsilon")>=0) :
	fact=""
	else :
	struct=struct.split("(",1)[1][:-1]
	if(fact[0]=="-") :
	sign="-"
	fact=fact[1:]
	split = struct.split("*")
	# find the epsilon
	eps=""
	for piece in split :
	if(piece.find("Epsilon")>=0) :
	eps=piece
	split.remove(piece)
	break
	# and any prefactors
	for piece in split :
	if(piece.find("P(")<0) :
	split.remove(piece)
	if(piece[0]=="+" or piece[0]=="-") :
	sign=changeSign(sign,piece[0])
	piece=piece[1:]
	if(fact=="") :
	fact=piece
	else :
	fact = "( %s ) * ( %s )" % (fact , piece)
	# now sort out the momenta
	for piece in split :
	terms=piece.split(",")
	terms[0]=terms[0].strip("P(")
	terms[1]=terms[1].strip(")")
	eps=eps.replace(terms[0],"P%s"%terms[1])
	(nsign,eps)=epsilonOrder(eps)
	if(nsign>0) : sign=changeSign(sign,"-")
	if(fact=="") : fact="1."
	if(eps!="Epsilon(1,2,P1,P2)") :
	return
	if(couplings[6]==0.) :
	couplings[6] = "( %s%s )" % (sign,fact)
	else :
	couplings[6] = "( %s ) + ( %s%s )" % (couplings[6],sign,fact)


	def scalarVectorCouplings(value,prefactors,L,lorentztag,all_couplings,order) :
	# set up the types of term we are looking for
	if(lorentztag=="VVS") :
	couplings=[0.,0.,0.,0.,0.,0.,0.]
	terms=[['P(-1,%s)' % order[0],
	'P(-1,%s)' % order[1],
	'Metric(%s,%s)' %(order[0],order[1])],
	['P(1,%s)' % order[0],
	'P(2,%s)' % order[0]],
	['P(1,%s)' % order[0],
	'P(2,%s)' % order[1]],
	['P(1,%s)' % order[1],
	'P(2,%s)' % order[0]],
	['P(1,%s)' % order[1],
	'P(2,%s)' % order[1]],
	['Metric(%s,%s)'%(order[0],order[1])]]
	elif(lorentztag=="VVSS") :
	couplings=[0.]
	terms=[['Metric(%s,%s)' % (order[0],order[1])]]
	elif(lorentztag=="VSS"):
	couplings=[0.,0.]
	terms=[['P(%s,%s)' % (order[0],order[2])],
	['P(%s,%s)' % (order[0],order[1])]]
	# extract the lorentz structures
	structures = extractStructures(L)
	# handle the scalar couplings
	itype=-1
	try :
	for term in terms:
	itype+=1
	for perm in itertools.permutations(term):
	label = '*'.join(perm)
	for istruct in range(0,len(structures)) :
	if label in structures[istruct] :
	reminder = structures[istruct].replace(label,'1.',1)
	couplings[itype]+=eval(reminder, {'cmath':cmath} )
	structures[istruct]='Done'
	except :
	raise SkipThisVertex()
	# special for VVS and epsilon
	# handle the pseudoscalar couplings
	for struct in structures :
	if(struct != "Done" ) :
	if(lorentztag=="VVS") :
	VVSEpsilon(couplings,struct)
	else :
	raise SkipThisVertex()
	# put it all together
	if(len(all_couplings)==0) :
	for ic in range(0,len(couplings)) :
	if(couplings[ic]!=0.) :
	all_couplings.append('(%s) * (%s) * (%s)' % (prefactors,value,couplings[ic]))
	else :
	all_couplings.append(False)
	else :
	for ic in range(0,len(couplings)) :
	if(couplings[ic]!=0. and all_couplings[ic]) :
	all_couplings[ic] = '(%s) * (%s) * (%s) + (%s) ' % (prefactors,value,
	couplings[ic],all_couplings[ic])
	elif(couplings[ic]!=0) :
	all_couplings[ic] = '(%s) * (%s) * (%s) ' % (prefactors,value,couplings[ic])
	return all_couplings

	def processScalarVectorCouplings(lorentztag,vertex,model,parmsubs,all_couplings,header,order) :
	# check the values
	tval = [False]*len(all_couplings[0])
	value =[False]*len(all_couplings[0])
	for icolor in range(0,len(all_couplings)) :
	for ix in range(0,len(all_couplings[icolor])) :
	if(not value[ix]) :
	value[ix] = all_couplings[icolor][ix]
	if(value[ix] and not tval[ix]) :
	tval[ix] = evaluate(value[ix],model,parmsubs)
	elif(value[ix]) :
	tval2 = evaluate(all_couplings[icolor][0],model,parmsubs)
	if(abs(tval[ix]-tval2)>1e-6) :
	raise SkipThisVertex()

	append = ""
	symbols = set()
	coup_norm=0.
	if(lorentztag=="VVS") :
	if(not value[0] and not value[1] and not value[2] and \
	not value[3] and not value[4] and not value[6] and value[5]) :
	coup_norm=value[5]
	else :
	for ix in range(0,len(value)) :
	if(value[ix]) :
	value[ix], sym = py2cpp(value[ix])
	symbols \|= sym
	else :
	value[ix]=0.
	lorentztag = 'GeneralVVS'
	header="kinematics(true);"
	# g_mu,nv piece of coupling
	if(value[5]!=0.) :
	append +='a00( %s + Complex(( %s )* GeV2/invariant(1,2)));\n' % ( value[0],value[5])
	else :
	append +='a00( %s );\n' % value[0]
	# other couplings
	append += 'a11( %s );\n a12( %s );\n a21( %s );\n a22( %s );\n aEp( %s );\n' % \
	( value[1],value[2],value[3],value[4],value[6] )
	coup_norm="1."
	elif(lorentztag=="VVSS") :
	coup_norm = value[0]
	elif(lorentztag=="VSS") :
	if(abs(tval[0]+tval[1])>1e-6) :
	for ix in range(0,len(value)) :
	if(value[ix]) :
	value[ix], sym = py2cpp(value[ix])
	symbols \|= sym
	else :
	value[ix]=0.
	coup_norm = "1."
	append = 'if(p2->id()==%s) { a( %s ); b( %s);}\n else { a( %s ); b( %s);}' \
	% (vertex.particles[order[1]-1].pdg_code,
	value[0],value[1],value[1],value[0])
	else :
	coup_norm = value[1]
	append = 'if(p2->id()!=%s){norm(-norm());}' \
	% vertex.particles[order[1]-1].pdg_code
	# return the answer
	return (coup_norm,append,lorentztag,header,symbols)

	def getIndices(term) :
	if(term[0:2]=="P(") :
	indices = term.strip(")").strip("P(").split(",")
	mom = int(indices[1])
	index = int(indices[0])
	return (True,mom,index)
	else :
	return (False,0,0)


	def lorentzScalar(vertex,L) :
	dotProduct = """(invariant( i[{i1}], i[{i2}] )/GeV2)"""
	structures=L.structure.split()
	output="("
	for struct in structures:
	if(struct=="+" or struct=="-") :
	output+=struct
	continue
	structure = struct.split("*")
	worked = False
	mom=-1
	newTerm=""
	while True :
	term = structure[-1]
	structure.pop()
	(momentum,mom,index) = getIndices(term)
	if( not momentum) : break
	# look for the matching momenta
	for term in structure :
	(momentum,mom2,index2) = getIndices(term)
	if(index2==index) :
	structure.remove(term)
	dot = dotProduct.format(i1=mom-1,i2=mom2-1)
	if(newTerm=="") :
	newTerm = dot
	else :
	newTerm = " ( %s) * ( %s ) " % (newTerm,dot)
	if(len(structure)==0) :
	worked = True
	break
	if(not worked) :
	return False
	else :
	output+=newTerm
	output+=")"
	return output

	kinematicsline = """\
	long id [3]={{{id1},{id2},{id3}}};
	long id2[3]={{p1->id(),p2->id(),p3->id()}};
	unsigned int i[3];
	for(unsigned int ix=0;ix<3;++ix) {{
	for(unsigned int iy=0;iy<3;++iy) {{
	if(id[ix]==id2[iy]) {{
	i[ix] = iy;
	id2[iy]=0;
	break;
	}}
	}}
	}}
	double hw_kine1 = {kine};
	"""

	kinematicsline2 = """\
	long id [4]={{{id1},{id2},{id3},{id4}}};
	long id2[4]={{p1->id(),p2->id(),p3->id(),p4->id()}};
	unsigned int i[4];
	for(unsigned int ix=0;ix<4;++ix) {{
	for(unsigned int iy=0;iy<4;++iy) {{
	if(id[ix]==id2[iy]) {{
	i[ix] = iy;
	id2[iy]=0;
	break;
	}}
	}}
	}}
	double hw_kine1 = {kine};
	"""

	kinematicsline3 ="""\
	double hw_kine{i} = {kine};
	"""

	def scalarCouplings(vertex,value,prefactors,L,lorentztag,
	all_couplings,prepend,header) :
	try :
	val = int(L.structure)
	except :
	output = lorentzScalar(vertex,L)
	if( not output ) :
	raise SkipThisVertex()
	else :
	if(prepend=="") :
	if(lorentztag=="SSS") :
	prepend = kinematicsline.format(id1=vertex.particles[0].pdg_code,
	id2=vertex.particles[1].pdg_code,
	id3=vertex.particles[2].pdg_code,
	kine=output)
	else :
	prepend = kinematicsline2.format(id1=vertex.particles[0].pdg_code,
	id2=vertex.particles[1].pdg_code,
	id3=vertex.particles[2].pdg_code,
	id4=vertex.particles[2].pdg_code,
	kine=output)
	value = "(%s) *(hw_kine1)" % value
	else :
	osplit=prepend.split("\n")
	i=-1
	while osplit[i]=="":
	i=i-1
	ikin=int(osplit[i].split("=")[0].replace("double hw_kine",""))+1
	prepend +=kinematicsline3.format(kine=output,i=ikin)
	value = "(%s) *(hw_kine%s)" % (value,ikin)
	header="kinematics(true);"
	if(len(all_couplings)==0) :
	all_couplings.append('(%s) * (%s)' % (prefactors,value))
	else :
	all_couplings[0] = '(%s) * (%s) + (%s)' % (prefactors,value,all_couplings[0])
	return (prepend, header,all_couplings)

	def processScalarCouplings(model,parmsubs,all_couplings) :
	tval = False
	value = False
	for icolor in range(0,len(all_couplings)) :
	if(len(all_couplings[icolor])!=1) :
	raise SkipThisVertex()
	if(not value) :
	value = all_couplings[icolor][0]
	m = re.findall('hw_kine[0-9]*', all_couplings[icolor][0])
	if m:
	for kine in m:
	# bizarre number for checks, must be a better option
	parmsubs[kine] = 987654321.
	if(not tval) :
	tval = evaluate(value,model,parmsubs)
	else :
	tval2 = evaluate(all_couplings[icolor][0],model,parmsubs)
	if(abs(tval[i]-tval2)>1e-6) :
	raise SkipThisVertex()
	# cleanup and return the answer
	return value.replace("(1.0) * ","").replace(" * (1)","")

	def vectorCouplings(vertex,value,prefactors,L,lorentztag,pos,
	all_couplings,append,qcd,order) :
	structures=extractStructures(L)
	terms=[]
	signs=[]
	if(lorentztag=="VVV") :
	terms=[['P(%s,%s)' % (order[2],order[0]),'Metric(%s,%s)' % (order[0],order[1])],
	['P(%s,%s)' % (order[2],order[1]),'Metric(%s,%s)' % (order[0],order[1])],
	['P(%s,%s)' % (order[1],order[0]),'Metric(%s,%s)' % (order[0],order[2])],
	['P(%s,%s)' % (order[1],order[2]),'Metric(%s,%s)' % (order[0],order[2])],
	['P(%s,%s)' % (order[0],order[1]),'Metric(%s,%s)' % (order[1],order[2])],
	['P(%s,%s)' % (order[0],order[2]),'Metric(%s,%s)' % (order[1],order[2])]]
	signs=[1.,-1.,-1.,1.,1.,-1.]
	elif(lorentztag=="VVVV") :
	terms=[['Metric(%s,%s)' % (order[0],order[3]),'Metric(%s,%s)' % (order[1],order[2])],
	['Metric(%s,%s)' % (order[0],order[2]),'Metric(%s,%s)' % (order[1],order[3])],
	['Metric(%s,%s)' % (order[0],order[1]),'Metric(%s,%s)' % (order[2],order[3])]]
	signs=[1.,1.,1.]
	elif(lorentztag=="VVVS") :
	terms=[['P(%s,%s)' % (order[2],order[0]),'Metric(%s,%s)' % (order[0],order[1])],
	['P(%s,%s)' % (order[2],order[1]),'Metric(%s,%s)' % (order[0],order[1])],
	['P(%s,%s)' % (order[1],order[0]),'Metric(%s,%s)' % (order[0],order[2])],
	['P(%s,%s)' % (order[1],order[2]),'Metric(%s,%s)' % (order[0],order[2])],
	['P(%s,%s)' % (order[0],order[1]),'Metric(%s,%s)' % (order[1],order[2])],
	['P(%s,%s)' % (order[0],order[2]),'Metric(%s,%s)' % (order[1],order[2])],
	['Epsilon(1,2,3,-1)','P(-1,1)'],['Epsilon(1,2,3,-1)','P(-1,2)'],
	['Epsilon(1,2,3,-1)','P(-1,3)']]
	signs=[1.,-1.,-1.,1.,1.,-1.,1.,1.,1.]

	# extract the couplings
	new_couplings = [False]*len(terms)
	iterm=0
	try :
	for term in terms:
	for perm in itertools.permutations(term):
	label = '*'.join(perm)
	for istruct in range(0,len(structures)) :
	if label in structures[istruct] :
	reminder = structures[istruct].replace(label,'1.',1)
	structures[istruct] = "Done"
	val = eval(reminder, {'cmath':cmath} )*signs[iterm]
	if(new_couplings[iterm]) :
	new_couplings[iterm] += val
	else :
	new_couplings[iterm] = val
	iterm += 1
	except :
	raise SkipThisVertex()
	# check we've handled all the terms
	for val in structures:
	if(val!="Done") :
	raise SkipThisVertex()
	# set the couplings
	for icoup in range(0,len(new_couplings)) :
	if(new_couplings[icoup]) :
	new_couplings[icoup] = '(%s) * (%s) *(%s)' % (new_couplings[icoup],prefactors,value)
	if(len(all_couplings)==0) :
	all_couplings=new_couplings
	else :
	for icoup in range(0,len(new_couplings)) :
	if(new_couplings[icoup] and all_couplings[icoup]) :
	all_couplings[icoup] = '(%s) * (%s) *(%s) + (%s) ' % (new_couplings[icoup],prefactors,value,all_couplings[icoup])
	elif(new_couplings[icoup]) :
	all_couplings[icoup] = new_couplings[icoup]
	# ordering for VVV type vertices
	if(len(pos[8]) != 3 and (lorentztag=="VVV" or lorentztag=="VVVS")) :
	append = VVVordering(vertex)
	return all_couplings,append

	def processVectorCouplings(lorentztag,vertex,model,parmsubs,all_couplings,append,header) :
	value = False
	tval = False
	if(lorentztag=="VVV") :
	for icolor in range(0,len(all_couplings)) :
	# loop over the different terms
	for ix in range(0,len(all_couplings[icolor])) :
	if(not value) :
	value = all_couplings[icolor][ix]
	tval = evaluate(value,model,parmsubs)
	else :
	tval2 = evaluate(all_couplings[icolor][ix],model,parmsubs)
	if(abs(tval-tval2)>1e-6) :
	raise SkipThisVertex()
	elif(lorentztag=="VVVV") :
	order=[]
	colours = vertex.color
	if(len(colours)==1) :
	tval=[]
	for i in range(0,3) :
	tval.append(evaluate(all_couplings[0][i],model,parmsubs))
	if(compare(tval[2],-2.*tval[1]) and
	compare(tval[2],-2.*tval[0]) ) :
	order=[0,1,2,3]
	value = "0.5*(%s)" % all_couplings[0][2]
	elif(compare(tval[1],-2.*tval[2]) and
	compare(tval[1],-2.*tval[0]) ) :
	order=[0,2,1,3]
	value = "0.5*(%s)" % all_couplings[0][1]
	elif(compare(tval[0],-2.*tval[2]) and
	compare(tval[0],-2.*tval[1]) ) :
	order=[0,3,1,2]
	value = "0.5*(%s)" % all_couplings[0][0]
	else:
	sys.stderr.write(
	'Warning: unsupported {tag} ( {ps} ) Lorentz structure in {name}:\n'
	.format(tag="VVVV", name=vertex.name, ps=' '.join(map(str,vertex.particles)))
	)
	raise SkipThisVertex()
	pattern = \
	"bool done[4]={false,false,false,false};\n" + \
	" tcPDPtr part[4]={p1,p2,p3,p4};\n" + \
	" unsigned int iorder[4]={0,0,0,0};\n" + \
	" for(unsigned int ix=0;ix<4;++ix) {\n" + \
	" if(!done[0] && part[ix]->id()==%s) {done[0]=true; iorder[%s] = ix; continue;}\n" + \
	" if(!done[1] && part[ix]->id()==%s) {done[1]=true; iorder[%s] = ix; continue;}\n" + \
	" if(!done[2] && part[ix]->id()==%s) {done[2]=true; iorder[%s] = ix; continue;}\n" + \
	" if(!done[3] && part[ix]->id()==%s) {done[3]=true; iorder[%s] = ix; continue;}\n" + \
	" }\n" + \
	" setType(2);\n" + \
	" setOrder(iorder[0],iorder[1],iorder[2],iorder[3]);"
	append = pattern % ( vertex.particles[0].pdg_code,order[0],
	vertex.particles[1].pdg_code,order[1],
	vertex.particles[2].pdg_code,order[2],
	vertex.particles[3].pdg_code,order[3] )
	else :
	for icolor in range(0,len(all_couplings)) :
	col=colours[icolor].split("*")
	if(len(col)==2 and "f(" in col[0] and "f(" in col[1]) :
	sign = 1
	for i in range(0,2) :
	col[i],stemp = extractAntiSymmetricIndices(col[i],"f(")
	for ix in range(0,len(col[i])): col[i][ix]=int(col[i][ix])
	sign *=stemp
	if(col[0][0]>col[1][0]) : col[0],col[1] = col[1],col[0]
	# first flow
	if(col[0][0]==1 and col[0][1]==2 and col[1][0] ==3 and col[1][1] == 4) :
	if(all_couplings[icolor][2] or not all_couplings[icolor][0] or
	not all_couplings[icolor][1]) :
	raise SkipThisVertex()
	if(not value) :
	value = all_couplings[icolor][0]
	tval = evaluate(value,model,parmsubs)
	tval2 = evaluate(all_couplings[icolor][0],model,parmsubs)
	tval3 = -evaluate(all_couplings[icolor][1],model,parmsubs)
	elif(col[0][0]==1 and col[0][1]==3 and col[1][0] ==2 and col[1][1] == 4) :
	if(all_couplings[icolor][1] or not all_couplings[icolor][0] or
	not all_couplings[icolor][2]) :
	raise SkipThisVertex()
	if(not value) :
	value = all_couplings[icolor][0]
	tval = evaluate(value,model,parmsubs)
	tval2 = evaluate(all_couplings[icolor][0],model,parmsubs)
	tval3 = -evaluate(all_couplings[icolor][2],model,parmsubs)
	elif(col[0][0]==1 and col[0][1]==4 and col[1][0] ==2 and col[1][1] == 3) :
	if(all_couplings[icolor][0] or not all_couplings[icolor][1] or
	not all_couplings[icolor][2]) :
	raise SkipThisVertex()
	if(not value) :
	value = all_couplings[icolor][1]
	tval = evaluate(value,model,parmsubs)
	tval2 = evaluate(all_couplings[icolor][1],model,parmsubs)
	tval3 = -evaluate(all_couplings[icolor][2],model,parmsubs)
	else :
	raise SkipThisVertex()
	if(abs(tval-tval2)>1e-6 or abs(tval-tval3)>1e-6 ) :
	raise SkipThisVertex()
	append = 'setType(1);\nsetOrder(0,1,2,3);'
	else :
	print 'unknown colour structure for VVVV vertex'
	raise SkipThisVertex()
	elif(lorentztag=="VVVS") :
	try :
	# two distinct cases 0-5 = , 6-8=
	if(all_couplings[0][0]) :
	imin=0
	imax=6
	header="scalar(true);"
	else :
	imin=6
	imax=9
	header="scalar(false);"
	for icolor in range(0,len(all_couplings)) :
	# loop over the different terms
	for ix in range(imin,imax) :
	if(not value) :
	value = all_couplings[icolor][ix]
	tval = evaluate(value,model,parmsubs)
	else :
	tval2 = evaluate(value,model,parmsubs)
	if(abs(tval-tval2)>1e-6) :
	raise SkipThisVertex()
	except :
	SkipThisVertex()
	# cleanup and return the answer
	value = value.replace("(1.0) * ","").replace(" * (1)","")
	return (value,append,header)

	def fermionCouplings(value,prefactors,L,all_couplings,order) :
	new_couplings=[False,False]
	try :
	new_couplings[0],new_couplings[1] = parse_lorentz(L.structure)
	except :
	raise SkipThisVertex()
	for i in range(0,2) :
	if new_couplings[i]:
	new_couplings[i] = '(%s) * (%s) * (%s)' % (prefactors,new_couplings[i],value)
	if(len(all_couplings)==0) :
	all_couplings=new_couplings
	else :
	for i in range(0,len(new_couplings)) :
	if(new_couplings[i] and all_couplings[i]) :
	all_couplings[i] = '(%s) + (%s) ' % (new_couplings[i],all_couplings[i])
	elif(new_couplings[i]) :
	all_couplings[i] = new_couplings[i]
	return all_couplings

	def processFermionCouplings(lorentztag,vertex,model,parmsubs,all_couplings) :
	leftcontent = all_couplings[0][0] if all_couplings[0][0] else "0."
	rightcontent = all_couplings[0][1] if all_couplings[0][1] else "0."
	tval=[evaluate( leftcontent,model,parmsubs),
	evaluate(rightcontent,model,parmsubs)]
	for icolor in range(0,len(all_couplings)) :
	# loop over the different terms
	for ix in range(0,len(all_couplings[icolor])) :
	tval2 = evaluate(all_couplings[icolor][ix],model,parmsubs) if all_couplings[icolor][ix] else 0.
	if(abs(tval[ix]-tval2)>1e-6) :
	raise SkipThisVertex()
	normcontent = "1."
	append=""
	if lorentztag == 'FFV':
	append = ('if(p1->id()!=%s) {Complex ltemp=left(), rtemp=right(); left(-rtemp); right(-ltemp);}'
	% vertex.particles[0].pdg_code)
	return normcontent,leftcontent,rightcontent,append

	def RSCouplings(value,prefactors,L,all_couplings,order) :
	raise SkipThisVertex()
	+
	+class LorentzStructure:
	+ """A simple example class to store a Lorentz structures"""
	+ name=""
	+ value=0.
	+ lorentz=[]
	+ spin=[]
	+ def __repr__(self):
	+ output = self.name
	+ if(self.name=="int" or self.name=="sign") :
	+ output += "=%s" % self.value
	+ elif(len(self.spin)==0) :
	+ output += "("
	+ for val in self.lorentz :
	+ output += "%s," % val
	+ output=output.rstrip(",")
	+ output+=")"
	+ elif(len(self.lorentz)==0) :
	+ output += "("
	+ for val in self.spin :
	+ output += "%s," % val
	+ output=output.rstrip(",")
	+ output+=")"
	+ else :
	+ output += "("
	+ for val in self.lorentz :
	+ output += "%s," % val
	+ for val in self.spin :
	+ output += "%s," % val
	+ output=output.rstrip(",")
	+ output+=")"
	+ return output
	+
	+def LorentzCompare(a,b) :
	+ if(a.name=="int" and b.name=="int") :
	+ return b.value-a.value
	+ elif(a.name=="int") :
	+ return -1
	+ elif(b.name=="int") :
	+ return 1
	+ elif(len(a.spin)==0) :
	+ if(len(b.spin)==0) :
	+ return len(b.lorentz)-len(a.lorentz)
	+ else :
	+ return -1
	+ elif(len(b.spin)==0) :
	+ return 1
	+ else :
	+ if(len(a.spin)==0 or len(b.spin)==0) :
	+ print 'index problem',a.name,b.name
	+ print a.spin,b.spin
	+ quit()
	+ if(a.spin[0]>0 or b.spin[1]>0 ) : return -1
	+ if(a.spin[1]>0 or b.spin[0]>0 ) : return 1
	+ if(a.spin[1]==b.spin[0]) : return -1
	+ if(b.spin[1]==a.spin[0]) : return 1
	+ return 0
	+
	+def extractIndices(struct) :
	+ if(struct.find("(")<0) : return []
	+ temp=struct.split("(")[1].split(")")[0].split(",")
	+ output=[]
	+ for val in temp :
	+ output.append(int(val))
	+ return output
	+
	+def parse_structure(structure) :
	+ output=[]
	+ # signs between terms
	+ if(structure=="+" or structure=="-") :
	+ output.append(LorentzStructure())
	+ output[0].name="sign"
	+ output[0].value=structure[0]+"1."
	+ output[0].value=float(output[0].value)
	+ return output
	+ # simple numeric pre/post factors
	+ elif((structure[0]=="-" or structure[0]=="+") and
	+ structure[-1]==")" and structure[1]=="(") :
	+ output.append(LorentzStructure())
	+ output[0].name="int"
	+ output[0].value=structure[0]+"1."
	+ output[0].value=float(output[0].value)
	+ structure=structure[2:-1]
	+ elif(structure[0]=="(") :
	+ temp=structure.rsplit(")",1)
	+ structure=temp[0][1:]
	+ output.append(LorentzStructure())
	+ output[0].name="int"
	+ output[0].value="1."+temp[1]
	+ output[0].value=float(eval(output[0].value))
	+ # split up the structure
	+ for struct in structure.split("*"):
	+ ind=extractIndices(struct)
	+ # different types of object
	+ # object with only spin indices
	+ if(struct.find("Identity")==0 or
	+ struct.find("Proj")==0 or
	+ struct.find("Gamma5")==0) :
	+ output.append(LorentzStructure())
	+ output[-1].spin=ind
	+ output[-1].name=struct.split("(")[0]
	+ output[-1].value=1.
	+ if(len(struct.replace("%s(%s,%s)" % (output[-1].name,ind[0],ind[1]),""))!=0) :
	+ print "problem A"
	+ quit()
	+ # objects with 2 lorentz indices
	+ elif(struct.find("Metric")==0 or struct.find("P(")==0) :
	+ output.append(LorentzStructure())
	+ output[-1].lorentz=ind
	+ output[-1].name=struct.split("(")[0]
	+ output[-1].value=1.
	+ if(len(struct.replace("%s(%s,%s)" % (output[-1].name,ind[0],ind[1]),""))!=0) :
	+ print "problem B"
	+ quit()
	+ # 1 lorentz and 1 spin index
	+ elif(struct.find("Gamma")==0) :
	+ output.append(LorentzStructure())
	+ output[-1].lorentz=[ind[0]]
	+ output[-1].spin=[ind[1],ind[2]]
	+ output[-1].name=struct.split("(")[0]
	+ output[-1].value=1.
	+ if(len(struct.replace("%s(%s,%s,%s)" % (output[-1].name,ind[0],ind[1],ind[2]),""))!=0) :
	+ print "problem C",struct
	+ quit()
	+ # objects with 4 lorentz indices
	+ elif(struct.find("Epsilon")==0) :
	+ output.append(LorentzStructure())
	+ output[-1].lorentz=ind
	+ output[-1].name=struct.split("(")[0]
	+ output[-1].value=1.
	+ if(len(struct.replace("%s(%s,%s,%s,%s)" % (output[-1].name,ind[0],ind[1],ind[2],ind[3]),""))!=0) :
	+ print "problem D"
	+ quit()
	+ # scalars
	+ else :
	+ try :
	+ output.append(LorentzStructure())
	+ output[-1].value=float(struct)
	+ output[0].name="int"
	+ except :
	+ print struct
	+ quit()
	+ # now do the sorting
	+ if(len(output)==1) : return output
	+ return sorted(output,cmp=LorentzCompare)
	+
	+def contractLorentz(L,parsed,lorentztag,order) :
	+ for l in range(0,len(parsed)) :
	+ for j in range(0,len(parsed[l])) :
	+ # replace indices with polarization vectors
	+ ll = len(parsed[l][j].lorentz)
	+ if(parsed[l][j].name=="P") :
	+ ll=1
	+ found=False
	+ for k in range(0,len(parsed[l])) :
	+ if(j==k) : continue
	+ for i in range(0,len(parsed[l][k].lorentz)) :
	+ if(parsed[l][k].lorentz[i]==parsed[l][j].lorentz[0]) :
	+ parsed[l][k].lorentz[i] = "P%s" % parsed[l][j].lorentz[1]
	+ if(parsed[l][k].name=="P") :
	+ parsed[l][k].lorentz[1] = "P%s" % parsed[l][k].lorentz[1]
	+ parsed[l][k].name="Metric"
	+ found=True
	+ break
	+ if(found) :
	+ parsed[l][j]=''
	+ ll=0
	+ for i in range(0,ll) :
	+ if(parsed[l][j]!="" and parsed[l][j].lorentz[i]>=0
	+ and isinstance(parsed[l][j].lorentz[i],(int,long)) and
	+ L.spins[parsed[l][j].lorentz[i]-1]==3 ) :
	+ parsed[l][j].lorentz[i]= "E%s" % parsed[l][j].lorentz[i]
	+ if(parsed[l][j].name=="P") :
	+ parsed[l][j].lorentz[1] = "P%s" % parsed[l][j].lorentz[1]
	+ parsed[l][j].name="Metric"
	+ parsed[l] = [x for x in parsed[l] if x != ""]
	+
	+def computeUnit(dimension) :
	+ dtemp=dimension[1]+dimension[2]
	+ if(dtemp==0) :
	+ unit="double"
	+ elif(dtemp==1) :
	+ unit="Energy"
	+ elif(dtemp==-1) :
	+ unit="InvEnergy"
	+ elif(dtemp>0) :
	+ unit="Energy%s" % (dtemp)
	+ elif(dtemp<0) :
	+ unit="InvEnergy%s" % (dtemp)
	+ return unit
	+
	+def computeUnit2(dimension,vDim) :
	+ # first correct for any coupling power in vertex
	+ totalDim = int(dimension[0])+dimension[2]+vDim-4
	+ output=""
	+ if(totalDim!=0) :
	+ if(totalDim>0) :
	+ if(totalDim==1) :
	+ output = "1./GeV"
	+ else :
	+ output = "1./GeV%s" % totalDim
	+ else :
	+ if(totalDim==-1) :
	+ output = "GeV"
	+ else :
	+ output = "GeV%s" % (-totalDim)
	+ expr=""
	+ # now remove the remaining dimensionality
	+ removal=dimension[1]-int(dimension[0])-vDim+4
	+ if(removal!=0) :
	+ if(removal>0) :
	+ if(removal==1) :
	+ expr = "UnitRemoval::InvE"
	+ else :
	+ expr = "UnitRemoval::InvE%s" % removal
	+ else :
	+ if(removal==-1) :
	+ expr = "UnitRemoval::E"
	+ else :
	+ expr = "UnitRemoval::E%s" % (-dimension)
	+ if(output=="") : return expr
	+ elif(expr=="") : return output
	+ else : return "%s*%s" %(output,expr)
	+
	+def generateVertex(iloc,L,parsed,lorentztag,order,defns) :
	+ # various strings for matrixes
	+ I4 = "Matrix([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])"
	+ G5 = "Matrix([[-1,0,0,0],[0,-1,0,0],[0,0,1,0],[0,0,0,1]])"
	+ PM = "Matrix([[1,0,0,0],[0,1,0,0],[0,0,0,0],[0,0,0,0]])"
	+ PP = "Matrix([[0,0,0,0],[0,0,0,0],[0,0,1,0],[0,0,0,1]])"
	+ vslash = Template("Matrix([[0,0,${v}tmz,-${v}xmy],[0,0,-${v}xpy,${v}tpz],[${v}tpz,${v}xmy,0,0],[${v}xpy,${v}tmz,0,0]])")
	+ vslashS = Template("${v}tpz=Symbol(\"${v}tpz\")\n${v}tmz=Symbol(\"${v}tmz\")\n${v}xpy=Symbol(\"${v}xpy\")\n${v}xmy=Symbol(\"${v}xmy\")\n")
	+ vslashD = Template("complex<${var}> ${v}tpz = ${v}.t()+${v}.z();\n complex<${var}> ${v}tmz = ${v}.t()-${v}.z();\n complex<${var}> ${v}xpy = ${v}.x()+Complex(0.,1.)${v}.y();\n complex<${var}> ${v}xmy = ${v}.x()-Complex(0.,1.)${v}.y();")
	+ vslashM = Template("Matrix([[$m,0,${v}tmz,-${v}xmy],[0,$m,-${v}xpy,${v}tpz],[${v}tpz,${v}xmy,$m,0],[${v}xpy,${v}tmz,0,$m]])")
	+ vslashM2 = Template("Matrix([[$m,0,-${v}tmz,${v}xmy],[0,$m,${v}xpy,-${v}tpz],[-${v}tpz,-${v}xmy,$m,0],[-${v}xpy,-${v}tmz,0,$m]])")
	+ vslashMS = Template("${v}tpz=Symbol(\"${v}tpz\")\n${v}tmz=Symbol(\"${v}tmz\")\n${v}xpy=Symbol(\"${v}xpy\")\n${v}xmy=Symbol(\"${v}xmy\")\n${m}=Symbol(\"${m}\")\n")
	+ dirac=["Matrix([[0,0,1,0],[0,0,0,1],[1,0,0,0],[0,1,0,0]])","Matrix([[0,0,0,1],[0,0,1,0],[0,-1,0,0],[-1,0,0,0]])",
	+ "Matrix([[0,0,0,complex(0, -1)],[0,0,complex(0, 1),0],[0,complex(0, 1),0,0],[complex(0, -1),0,0,0]])",
	+ "Matrix([[0,0,1,0],[0,0,0,-1],[-1,0,0,0],[0,1,0,0]])"]
	+ # order the indices of a dot product
	+ def indSort(a,b) :
	+ if(a[0]==b[0]) :
	+ i1=int(a[1])
	+ i2=int(b[1])
	+ if(i1>i2) :
	+ return 1
	+ elif(i1<i2) :
	+ return -1
	+ else :
	+ return 0
	+ else :
	+ if(a[0]=="E") :
	+ return 1
	+ else :
	+ return -1
	+ # parse the lorentz structures
	+ output = [1.]*len(parsed)
	+ dimension=[]
	+ for i in range(0,len(parsed)) :
	+ dimension.append([0,0,0])
	+ for i in range (0,len(parsed)) :
	+ # replace signs
	+ if(len(parsed[i])==1 and parsed[i][0].name=="sign") :
	+ if(parsed[i][0].value>0) :
	+ output[i]="+"
	+ else :
	+ output[i]="-"
	+ parsed[i][0]=''
	+ continue
	+ # replace integers
	+ for j in range(0,len(parsed[i])) :
	+ if(parsed[i][j]!="" and parsed[i][j].name=="int") :
	+ output[i] *= parsed[i][j].value
	+ parsed[i][j]=""
	+ continue
	+ output[i] = "(%s)" % output[i]
	+ for j in range(0,len(parsed[i])) :
	+ if(parsed[i][j]!="" and parsed[i][j].name=="Metric") :
	+ (ind1,ind2) = sorted((parsed[i][j].lorentz[0],parsed[i][j].lorentz[1]),cmp=indSort)
	+ # this product already dealt with ?
	+ if((ind1,ind2) in defns) :
	+ output[i] += "*(%s)" % defns[(ind1,ind2)][0]
	+ parsed[i][j]=""
	+ if(ind1[0]=="P") : dimension[i][2] +=1
	+ if(ind1[1]=="P") : dimension[i][2] +=1
	+ continue
	+ # handle the product
	+ name = "dot%s" % (len(defns)+1)
	+ parsed[i][j]=""
	+ if(ind1[0]=="P") :
	+ # dot product of two momenta
	+ if(ind2[0]=="P") :
	+ dimension[i][2] +=2
	+ defns[(ind1,ind2)] = [name,"Energy2 %s = %s*%s;" % (name,ind1,ind2)]
	+ output[i] += "*(%s)" % name
	+ elif(ind2[0]=="E") :
	+ dimension[i][2] +=1
	+ if(int(ind2[1])==iloc) :
	+ output[i] += "*(%s)" % ind1
	+ else :
	+ defns[(ind1,ind2)] = [name,"complex<Energy> %s = %s*%s;" % (name,ind1,ind2)]
	+ output[i] += "*(%s)" % name
	+ elif(ind1[0]=="E") :
	+ if(ind2[0]!="E") :
	+ print "problem"
	+ quit()
	+ if(int(ind1[1])==iloc) :
	+ output[i] += "*(%s)" % ind2
	+ elif(int(ind2[1])==iloc) :
	+ output[i] += "*(%s)" % ind1
	+ else :
	+ defns[(ind1,ind2)] = [name,"complex<double> %s = %s*%s;" % (name,ind1,ind2)]
	+ output[i] += "*(%s)" % name
	+ # remove any (now) empty elements
	+ for i in range (0,len(parsed)) :
	+ parsed[i] = [x for x in parsed[i] if x != ""]
	+ # now for gamma matrix strings
	+ if(lorentztag[0]=="F") :
	+ result=""
	+ for i in range(0,len(parsed)):
	+ if(len(parsed[i])==0) :
	+ continue
	+ if(len(parsed[i])==1 and parsed[i][0]=="") :
	+ result+=parsed[i][0]
	+ continue
	+ # first and last lorentz indices
	+ sind=parsed[i][0].spin[0]
	+ lind=0
	+ # first piece of the expression we need to evaluate
	+ dtemp=[0,0,0]
	+ if(sind==iloc) :
	+ expr = vslashM2.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind} )
	+ Symbols = vslashMS.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind} )
	+ defns["vvP%s" % sind ] = ["vvP%s" % sind ,
	+ vslashD.substitute({ "var" : "Energy",
	+ "v" : "P%s" % sind })]
	+ dtemp[1]+=1
	+ dtemp[0]+=0.5
	+ else :
	+ Symbols=Template("sbar${s}s1=Symbol(\"sbar${s}s1\")\nsbar${s}s2=Symbol(\"sbar${s}s2\")\nsbar${s}s3=Symbol(\"sbar${s}s3\")\nsbar${s}s4=Symbol(\"sbar${s}s4\")\n").substitute({'s' :parsed[i][0].spin[0]})
	+ expr=Template("Matrix([[sbar${s}s1,sbar${s}s2,sbar${s}s3,sbar${s}s4]])").substitute({'s' : sind })
	+ dtemp[0]+=0.5
	+ # parse the remaining structures
	+ for j in range(0,len(parsed[i])) :
	+ if(parsed[i][j].name=="Identity") :
	+ expr += "*%s" % I4
	+ lind = parsed[i][j].spin[1]
	+ parsed[i][j]=""
	+ elif(parsed[i][j].name=="Gamma5") :
	+ expr += "*%s" % G5
	+ lind = parsed[i][j].spin[1]
	+ parsed[i][j]=""
	+ elif(parsed[i][j].name=="ProjM") :
	+ expr += "*%s" % PM
	+ lind = parsed[i][j].spin[1]
	+ parsed[i][j]=""
	+ elif(parsed[i][j].name=="ProjP") :
	+ expr += "*%s" % PP
	+ lind = parsed[i][j].spin[1]
	+ parsed[i][j]=""
	+ elif(parsed[i][j].name=="Gamma") :
	+ lind = parsed[i][j].spin[1]
	+ # lorentz matrix contracted with the propagator
	+ if(parsed[i][j].lorentz[0] == ("E%s" % iloc ) ) :
	+ expr += "*DUMMY"
	+ else :
	+ expr += "*"+vslash.substitute({ "v" : parsed[i][j].lorentz[0]})
	+ Symbols += vslashS.substitute({ "v" : parsed[i][j].lorentz[0]})
	+ if(parsed[i][j].lorentz[0][0]=="P") :
	+ dtemp[2] += 1
	+ variable="Energy"
	+ else :
	+ variable="double"
	+ defns["vv%s" % parsed[i][j].lorentz[0] ] = \
	+ ["vv%s" % parsed[i][j].lorentz[0],
	+ vslashD.substitute({ "var" : variable,
	+ "v" : parsed[i][j].lorentz[0]})]
	+ parsed[i][j]=""
	+ else :
	+ print 'FFFFFF'
	+ print parsed[i][j]
	+ quit()
	+ # last piece of spin chain
	+ if(lind==iloc) :
	+ expr += "*"+vslashM.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind} )
	+ Symbols += vslashMS.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind} )
	+ defns["vvP%s" % lind ] = ["vvP%s" % lind ,
	+ vslashD.substitute({ "var" : "Energy",
	+ "v" : "P%s" % lind })]
	+ dtemp[1] += 1
	+ dtemp[0] += 0.5
	+ else :
	+ expr+=Template("*Matrix([[s${s}s1],[s${s}s2],[s${s}s3],[s${s}s4]])").substitute({'s' : lind})
	+ Symbols+=Template("s${s}s1=Symbol(\"s${s}s1\")\ns${s}s2=Symbol(\"s${s}s2\")\ns${s}s3=Symbol(\"s${s}s3\")\ns${s}s4=Symbol(\"s${s}s4\")\n").substitute({'s' : lind})
	+ dtemp[0] += 0.5
	+ parsed[i] = [x for x in parsed[i] if x != ""]
	+ if(len(parsed[i])!=0) :
	+ print "ERROR"
	+ quit()
	+ # off-shell vector
	+ if(expr.find("DUMMY")>=0) :
	+ vtemp=[]
	+ defns["I"] = ["I","static Complex I(0.,1.);"]
	+ for matrix in dirac :
	+ temp={}
	+ exec("import sympy\nfrom sympy import Symbol,Matrix\n"+Symbols+"result="+expr.replace("DUMMY",matrix)) in temp
	+ vtemp.append(temp["result"][0,0])
	+ unit = computeUnit(dtemp)
	+ output[i] += "*LorentzVector<complex<%s> >(%s,%s,%s,%s)" % (unit,vtemp[1],vtemp[2],vtemp[3],vtemp[0])
	+ else :
	+ temp={}
	+ if(iloc==0 or (iloc!=sind and iloc!=lind)) :
	+ exec("import sympy\nfrom sympy import Symbol,Matrix\n"+Symbols+"result="+expr) in temp
	+ output[i] += "*(%s)" % temp["result"][0,0]
	+ else :
	+ exec("import sympy\nfrom sympy import Symbol,Matrix\n"+Symbols+"result="+expr) in temp
	+ unit = computeUnit(dtemp)
	+ if(iloc==sind) :
	+ output[i] += "*LorentzSpinor<%s>(%s,%s,%s,%s)" % (unit,temp["result"][0],temp["result"][1],
	+ temp["result"][2],temp["result"][3])
	+ else :
	+ output[i] += "*LorentzSpinorBar<%s >(%s,%s,%s,%s)" % (unit,temp["result"][0],temp["result"][1],
	+ temp["result"][2],temp["result"][3])
	+ dimension[i] = list(map(lambda x, y: x + y, dtemp, dimension[i]))
	+ return (output,dimension)
	+
	+def convertLorentz(Lstruct,lorentztag,order,iloc,defns,evalVertex) :
	+ # split the structure into individual terms
	+ structures=Lstruct.structure.split()
	+ parsed=[]
	+ for struct in structures :
	+ parsed.append(parse_structure(struct))
	+ # convert lorentz contractions to dot products
	+ contractLorentz(Lstruct,parsed,lorentztag,order)
	+ # now in a position to generate the code
	+ evalVertex.append(generateVertex(iloc,Lstruct,parsed,lorentztag,order,defns))
	+
	+evaluateTemplate = """\
	+{decl} {{
	+ {momenta}
	+ {waves}
	+ {defns}
	+ {symbols}
	+ {couplings}
	+ {result}
	+}}
	+"""
	+
	+def generateEvaluateFunction(model,vertex,iloc,values,defns,vertexEval,cf) :
	+ # first construct the signature of the function
	+ iferm=0
	+ decls=[]
	+ offType="Complex"
	+ momenta=[]
	+ waves=[]
	+ poff=""
	+ fermionReplace=[]
	+ for i in range(0,len(vertex.lorentz[0].spins)) :
	+ spin = vertex.lorentz[0].spins[i]
	+ if(i+1==iloc) :
	+ if(spin==1) :
	+ offType="ScalarWaveFunction"
	+ elif(spin==2) :
	+ if(iferm==0) :
	+ offType="SpinorBarWaveFunction"
	+ else :
	+ offType="SpinorWaveFunction"
	+ iferm+=1
	+ elif(spin==3) :
	+ offType="VectorWaveFunction"
	+ elif(spin==4) :
	+ offType="TensorWaveFunction"
	+ else :
	+ print 'unknown spin',spin
	+ quit()
	+ poff = ("Lorentz5Momentum P%s = " % (i+1) ) + poff
	+ else :
	+ if(spin==1) :
	+ decls.append("ScalarWaveFunction & scaW%s" % (i+1))
	+ momenta.append("Lorentz5Momentum P%s = scaW%s.momentum();" % (i+1,i+1))
	+ waves.append("Complex sca%s = scaW%s.wave();" % (i+1,i+1))
	+ elif(spin==2) :
	+ if(iferm==0) :
	+ decls.append("SpinorWaveFunction & sW%s" % (i+1))
	+ momenta.append("Lorentz5Momentum P%s = sW%s.momentum();" % (i+1,i+1))
	+ waves.append("LorentzSpinor<double> s%s = sW%s.wave();" % (i+1,i+1))
	+ fermionReplace.append("s%s"%(i+1))
	+ else :
	+ decls.append("SpinorBarWaveFunction & sbarW%s" % (i+1))
	+ momenta.append("Lorentz5Momentum P%s = sbarW%s.momentum();" % (i+1,i+1))
	+ waves.append("LorentzSpinorBar<double> sbar%s = sbarW%s.wave();" % (i+1,i+1))
	+ fermionReplace.append("sbar%s"%(i+1))
	+ iferm +=1
	+ elif(spin==3) :
	+ decls.append("VectorWaveFunction & vW%s" % (i+1))
	+ momenta.append("Lorentz5Momentum P%s = vW%s.momentum();" % (i+1,i+1))
	+ waves.append("LorentzPolarizationVector E%s = vW%s.wave();" % (i+1,i+1))
	+ elif(spin==4) :
	+ decls.append("TensorWaveFunction & tW%s" % (i+1))
	+ momenta.append("Lorentz5Momentum P%s = tW%s.momentum();" % (i+1,i+1))
	+ waves.append("LorentzTensor t%s = tW%s.wave();" % (i+1,i+1))
	+ else :
	+ print 'unknown spin',spin
	+ quit()
	+ poff += "-P%s" % (i+1)
	+ sig=""
	+ if(iloc==0) :
	+ sig="%s evaluate(Energy2, const %s)" % (offType,", const ".join(decls))
	+ else :
	+ sig="%s evaluate(Energy2, int iopt, tcPDPtr out, const %s, complex<Energy> mass=-GeV, complex<Energy> width=-GeV)" % (offType,", const ".join(decls))
	+ momenta.append(poff+";")
	+ # cat the definitions
	+ defString=""
	+ for (key,value) in defns.iteritems() :
	+ defString+=" %s\n" %value[1]
	+
	+ oval=""
	+ symbols=set()
	+ localCouplings=[]
	+ result=""
	+ for j in range(0,len(vertexEval)) :
	+ (vals,dim) = vertexEval[j]
	+ expr=""
	+ dimCheck=dim[0]
	+ for i in range(0,len(vals)) :
	+ if(vals[i]=="+" or vals[i]=="-") :
	+ expr +=vals[i]
	+ else :
	+ if(dimCheck[0]!=dim[i][0] or dimCheck[1]!=dim[i][1] or
	+ dimCheck[2]!=dim[i][2]) :
	+ print "DIMENSION PROBLEM",dimCheck,dim[i],vertex
	+ quit()
	+ expr += "(%s)" % vals[i]
	+ for rep in fermionReplace :
	+ for i in range(1,5) :
	+ oldVal = "%ss%s" % (rep,i)
	+ newVal = "%s.s%s()" % (rep,i)
	+ expr=expr.replace(oldVal,newVal)
	+ vDim = len(vertex.lorentz[0].spins)
	+
	+ unit = computeUnit2(dimCheck,vDim)
	+ if(unit!="") :
	+ expr = "(%s)*(%s)" % (expr,unit)
	+ val, sym = py2cpp(values[j])
	+ localCouplings.append("Complex local_C%s = %s;\n" % (j,val))
	+ symbols \|=sym
	+ if(result!="") :
	+ result += " + (local_C%s)*(%s)" % (j,expr)
	+ else :
	+ result += " (local_C%s)*(%s) " % (j,expr)
	+ if(iloc==0) :
	+ result = "return (%s)*(%s);\n" % (result,py2cpp(cf[0])[0])
	+ else :
	+ if(vertex.lorentz[0].spins[iloc-1] == 3 ) :
	+ result = "LorentzPolarizationVector vtemp = %s;\n" % result
	+ result +=" Energy2 p2 = P%s.m2();\nComplex fact = -Complex(0.,1.)(%s)propagator(iopt,p2,out,mass,width);\n if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass();\n complex<Energy2> mass2 = sqr(mass);\n if(mass.real()==ZERO) {\n vtemp =factvtemp;\n}\n else {\ncomplex<Energy> dot = P%svtemp;\n vtemp = fact(vtemp-dot/mass2P%s);\n}\nreturn VectorWaveFunction(-P%s,out,vtemp.x(),vtemp.y(),vtemp.z(),vtemp.t());\n" % (iloc,py2cpp(cf[0])[0],iloc,iloc,iloc)
	+ elif(vertex.lorentz[0].spins[iloc-1] == 2 ) :
	+ result = "if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass();\n Energy2 p2 = P%s.m2();\n Complex fact = Complex(0.,1.)(%s)propagator(iopt,p2,out,mass,width);\n Lorentz%s<double> newSpin = fact*(%s);\n return %s(-P%s,out,newSpin.s1(),newSpin.s2(),newSpin.s3(),newSpin.s4());" % \
	+ (iloc,py2cpp(cf[0])[0],offType.replace("WaveFunction",""),result.replace( "M%s" % iloc, "mass" ),offType,iloc)
	+
	+ header="virtual %s" % sig
	+ sig=sig.replace("=-GeV","")
	+ symboldefs = [ def_from_model(model,s) for s in symbols ]
	+ function = evaluateTemplate.format(decl=sig,momenta="\n ".join(momenta),defns=defString,
	+ waves="\n ".join(waves),symbols='\n '.join(symboldefs),
	+ couplings="\n ".join(localCouplings),
	+ result=result)
	+ return (header,function)
	diff --git a/Models/Feynrules/python/ufo2peg/vertices.py b/Models/Feynrules/python/ufo2peg/vertices.py
	--- a/Models/Feynrules/python/ufo2peg/vertices.py
	+++ b/Models/Feynrules/python/ufo2peg/vertices.py
	@@ -1,663 +1,683 @@
	import sys,pprint

	from .helpers import CheckUnique,getTemplate,writeFile,qcd_qed_orders,def_from_model
	from .converter import py2cpp
	from .collapse_vertices import collapse_vertices
	from .check_lorentz import tensorCouplings,VVVordering,lorentzScalar,\
	processTensorCouplings,scalarCouplings,processScalarCouplings,scalarVectorCouplings,\
	processScalarVectorCouplings,vectorCouplings,processVectorCouplings,fermionCouplings,processFermionCouplings,\
	- RSCouplings
	+ RSCouplings,convertLorentz,generateEvaluateFunction
	from .helpers import SkipThisVertex,extractAntiSymmetricIndices,isGoldstone

	# prefactors for vertices
	lfactors = {
	'FFV' : '-complex(0,1)', # ok
	'VVV' : 'complex(0,1)', # changed to fix ttbar
	'VVVS' : 'complex(0,1)', # should be as VVV
	'VVVV' : 'complex(0,1)',
	'VVS' : '-complex(0,1)',
	'VSS' : '-complex(0,1)', # changed to minus to fix dL ->x1 W- d
	'SSS' : '-complex(0,1)', # ok
	'VVSS' : '-complex(0,1)', # ok
	'VVT' : 'complex(0,2)',
	'VVVT' : '-complex(0,2)',
	'SSSS' : '-complex(0,1)', # ok
	'FFS' : '-complex(0,1)', # ok
	'SST' : 'complex(0,2)',
	'FFT' : '-complex(0,8)',
	'FFVT' : '-complex(0,4)',
	'RFS' : 'complex(0,1)',
	'RFV' : 'complex(0,1)',
	}

	+genericVertices=['FFVV']
	+
	# template for the header for a vertex
	VERTEXHEADER = """\
	#include "ThePEG/Helicity/Vertex/{spindirectory}/{lorentztag}Vertex.h"
	"""
	+GENERALVERTEXHEADER = """\
	+#include "ThePEG/Helicity/Vertex/Abstract{lorentztag}Vertex.h"
	+"""

	# template for the implmentation for a vertex
	VERTEXCLASS = getTemplate('Vertex_class')
	+GENERALVERTEXCLASS = getTemplate('GeneralVertex_class')

	# template for the .cc file for vertices
	VERTEX = getTemplate('Vertex.cc')

	vertexline = """\
	create Herwig::{modelname}V_{vname} /Herwig/{modelname}/V_{vname}
	insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_{vname}
	"""


	def get_lorentztag(spin):
	"""Produce a ThePEG spin tag for the given numeric FR spins."""
	spins = { 1 : 'S', 2 : 'F', 3 : 'V', 4 : 'R', 5 : 'T', -1 : 'U' }
	result=[]
	for i in range(0,len(spin)) :
	result.append((spins[spin[i]],i+1))
	def spinsort(a,b):
	"""Helper function for ThePEG's FVST spin tag ordering."""
	(a1,a2) = a
	(b1,b2) = b
	if a1 == b1: return 0
	for letter in 'URFVST':
	if a1 == letter: return -1
	if b1 == letter: return 1

	result = sorted(result, cmp=spinsort)
	order=[]
	output=""
	for i in range(0,len(result)) :
	(a,b) = result[i]
	order.append(b)
	output+=a
	return (output,order)

	def unique_lorentztag(vertex):
	"""Check and return the Lorentz tag of the vertex."""
	unique = CheckUnique()
	for l in vertex.lorentz:
	(lorentztag,order) = get_lorentztag(l.spins)
	unique( lorentztag )
	lname = l.name[:len(lorentztag)]
	if sorted(lorentztag) != sorted(lname):
	raise Exception("Lorentztags: %s is not %s in %s"
	% (lorentztag,lname,vertex))
	return (lorentztag,order)

	def colors(vertex) :
	try:
	unique = CheckUnique()
	for pl in vertex.particles_list:
	struct = [ p.color for p in pl ]
	unique(struct)
	except:
	struct = [ p.color for p in vertex.particles ]
	pos = colorpositions(struct)
	L = len(struct)
	return (L,pos)

	def colorfactor(vertex,L,pos,lorentztag):
	def match(patterns,color=vertex.color):
	result = [ p == t
	for p,t in zip(patterns,color) ]
	return all(result)

	label = None
	l = lambda c: len(pos[c])
	if l(1) == L:
	label = ('1',)
	- if match(label): return ('1',)
	+ if match(label): return ('1.',)

	elif l(3) == l(-3) == 1 and l(1) == L-2:
	nums = [pos[3][0], pos[-3][0]]
	label = ('Identity({0},{1})'.format(*sorted(nums)),)
	- if match(label): return ('1',)
	+ if match(label): return ('1.',)

	elif l(6) == l(-6) == 1 and l(1) == L-2:
	nums = [pos[6][0], pos[-6][0]]
	label = ('Identity({0},{1})'.format(*sorted(nums)),)
	- if match(label): return ('1',)
	+ if match(label): return ('1.',)

	elif l(6) == l(-6) == 2 and L==4:
	sys.stderr.write(
	'Warning: Unknown colour structure 6 6 6~ 6~ ( {ps} ) in {name}.\n'
	.format(name=vertex.name, ps=' '.join(map(str,vertex.particles)))
	)
	raise SkipThisVertex()

	elif l(8) == l(3) == l(-3) == 1 and l(1) == L-3:
	label = ('T({g},{q},{qb})'.format(g=pos[8][0],q=pos[3][0],qb=pos[-3][0]),)
	- if match(label): return ('1',)
	+ if match(label): return ('1.',)

	elif l(8) == l(6) == l(-6) == 1 and l(1) == L-3:
	label = ('T6({g},{s},{sb})'.format(g=pos[8][0],s=pos[6][0],sb=pos[-6][0]),)
	- if match(label): return ('1',)
	+ if match(label): return ('1.',)

	elif l(6) == 1 and l(-3) == 2 and L==3:
	label = ('K6({s},{qb1},{qb2})'.format(s=pos[6][0],qb1=pos[-3][0],qb2=pos[-3][1]),)
	- if match(label): return ('1',)
	+ if match(label): return ('1.',)

	elif l(-6) == 1 and l(3) == 2 and L==3:
	label = ('K6Bar({sb},{q1},{q2})'.format(sb=pos[-6][0],q1=pos[3][0],q2=pos[3][1]),)
	- if match(label): return ('1',)
	+ if match(label): return ('1.',)

	elif l(3) == L == 3:
	colors=[]
	for color in vertex.color :
	order,sign = extractAntiSymmetricIndices(color,"Epsilon(")
	colors.append("Epsilon(%s,%s,%s)" % (order[0],order[1],order[2]))
	label = ('Epsilon(1,2,3)',)
	- if match(label,colors): return ('1',) # TODO check factor!
	+ if match(label,colors): return ('1.',) # TODO check factor!

	elif l(-3) == L == 3:
	colors=[]
	for color in vertex.color :
	order,sign = extractAntiSymmetricIndices(color,"EpsilonBar(")
	colors.append("Epsilon(%s,%s,%s)" % (order[0],order[1],order[2]))
	label = ('EpsilonBar(1,2,3)',)
	- if match(label): return ('1',) # TODO check factor!
	+ if match(label): return ('1.',) # TODO check factor!

	elif l(8) == L == 3:
	colors=[]
	for color in vertex.color :
	order,sign = extractAntiSymmetricIndices(color,"f(")
	colors.append("f(%s,%s,%s)" % (order[0],order[1],order[2]))
	# if lorentz is FFV get extra minus sign
	if lorentztag in ['FFV'] : sign *=-1
	label = ('f(1,2,3)',)
	- if match(label,colors): return ('-complex(0,1)*(%s)'%sign,)
	+ if match(label,colors): return ('-complex(0,1.)*(%s)'%sign,)

	elif l(8) == L == 4:
	colors=[]
	for color in vertex.color :
	f = color.split("*")
	(o1,s1) = extractAntiSymmetricIndices(f[0],"f(")
	(o2,s2) = extractAntiSymmetricIndices(f[1],"f(")
	if(o2[0]<o1[0]) : o1,o2=o2,o1
	colors.append("f(%s)*f(%s)" % (",".join(o1),",".join(o2)))
	def coloursort(a,b) :
	if a == b: return 0
	i1=int(a[4])
	i2=int(b[4])
	if(i1==i2) : return 0
	elif(i1<i2) : return -1
	else : return 1
	colors=sorted(colors,cmp=coloursort)
	label = ('f(1,2,-1)*f(3,4,-1)',
	'f(1,3,-1)*f(2,4,-1)',
	'f(1,4,-1)*f(2,3,-1)')
	nmatch=0
	for c1 in label:
	for c2 in colors :
	if(c1==c2) : nmatch+=1
	if(nmatch==2 and lorentztag=="VVSS") :
	- return ('1','1')
	+ return ('1.','1.')
	elif(nmatch==3 and lorentztag=="VVVV") :
	return ('-1.','-1.','-1.')

	elif l(8) == 2 and l(3) == l(-3) == 1 and L==4:
	subs = {
	'g1' : pos[8][0],
	'g2' : pos[8][1],
	'qq' : pos[3][0],
	'qb' : pos[-3][0]
	}
	- label = ('T({g1},-1,{qb})T({g2},{qq},-1)'.format(*subs),
	- 'T({g1},{qq},-1)T({g2},-1,{qb})'.format(*subs))
	- if match(label): return ('1.','1.')
	+ if(vertex.lorentz[0].spins.count(1)==2) :
	+ label = ('T({g1},-1,{qb})T({g2},{qq},-1)'.format(*subs),
	+ 'T({g1},{qq},-1)T({g2},-1,{qb})'.format(*subs))
	+ if match(label): return ('1.','1.')
	+ elif(vertex.lorentz[0].spins.count(2)==2) :
	+ label = ('f({g1},{g2},-1)T(-1,{qq},{qb})'.format(*subs),)
	+ if match(label): return ('complex(0.,1.)',)

	elif l(8) == 2 and l(6) == l(-6) == 1 and L==4:
	subs = {
	'g1' : pos[8][0],
	'g2' : pos[8][1],
	'qq' : pos[6][0],
	'qb' : pos[-6][0]
	}
	label = ('T6({g1},-1,{qb})T6({g2},{qq},-1)'.format(*subs),
	'T6({g1},{qq},-1)T6({g2},-1,{qb})'.format(*subs))
	if match(label): return ('1.','1.')

	elif l(8) == 2 and l(8)+l(1)==L :
	subs = { 'g1' : pos[8][0], 'g2' : pos[8][1] }
	label = ('Identity({g1},{g2})'.format(**subs),)
	if match(label) : return ('1.',)

	elif l(8) == 3 and l(1)==1 and L==4 :
	colors=[]
	for color in vertex.color :
	order,sign = extractAntiSymmetricIndices(color,"f(")
	colors.append("f(%s,%s,%s)" % (order[0],order[1],order[2]))
	label = ('f(1,2,3)',)
	- if match(label,colors): return ('-complex(0,1)*(%s)'%sign,)
	+ if match(label,colors): return ('-complex(0.,1.)*(%s)'%sign,)

	sys.stderr.write(
	"Warning: Unknown colour structure {color} ( {ps} ) in {name}.\n"
	.format(color = ' '.join(vertex.color), name = vertex.name,
	ps = ' '.join(map(str,vertex.particles)))
	)
	raise SkipThisVertex()

	def colorpositions(struct):
	positions = {
	1 : [],
	3 : [],
	-3 : [],
	6 : [],
	-6 : [],
	8 : [],
	}
	for i,s in enumerate(struct,1):
	positions[s].append(i)
	return positions

	def spindirectory(lt):
	"""Return the spin directory name for a given Lorentz tag."""
	if 'T' in lt:
	spin_directory = 'Tensor'
	elif 'S' in lt:
	spin_directory = 'Scalar'
	elif 'V' in lt:
	spin_directory = 'Vector'
	else:
	raise Exception("Unknown Lorentz tag {lt}.".format(lt=lt))
	return spin_directory

	def write_vertex_file(subs):
	'Write the .cc file for some vertices'
	newname = '%s_Vertices_%03d.cc' % (subs['ModelName'],subs['vertexnumber'])
	subs['newname'] = newname
	writeFile( newname, VERTEX.substitute(subs) )

	def checkGhostGoldstoneVertex(lorentztag,vertex) :
	'check if vertex has ghosts or goldstones'
	# remove vertices involving ghost fields
	if 'U' in lorentztag:
	return True
	# remove vertices involving goldstones
	for p in vertex.particles:
	if(isGoldstone(p)):
	return True
	return False

	-def calculatePrefactor(globalsign,lorentztag,lf,cf) :
	- if(globalsign!=1.) :
	- prefactors = '(%s) * (%s) * (%s)' \
	- % (globalsign**(len(lorentztag)-2),lf,cf)
	- else :
	- prefactors = '(%s) * (%s)' \
	- % (lf,cf)
	- return prefactors
	- # fact=[]
	- # if(globalsign!=1.) :
	- # fact.append(globalsign**(len(lorentztag)-2))
	- # if(lf!="1") :
	- # fact.append(lf)
	- # if(cf!="1") :
	- # fact.append(cf)
	- # if(len(fact)==0) : return "1"
	- # prefactor = '(%s)' % fact[0]
	- # for ix in range(1,len(fact)) :
	- # prefactor = '%s * (%s)' % (prefactor,fact[ix])
	- # return prefactor
	+def calculatePrefactor(lf,cf) :
	+ prefactor = '(%s) * (%s)' % (lf,cf)
	+ return prefactor

	def couplingValue(coupling) :
	if type(coupling) is not list:
	value = coupling.value
	else:
	value = "("
	for coup in coupling :
	value += '+(%s)' % coup.value
	value +=")"
	return value

	def epsilonSign(vertex,couplingptrs,append) :
	EPSSIGN = """\
	double sign = {epssign};
	if((p1->id()=={id1} && p2->id()=={id3} && p3->id()=={id2}) \|\|
	(p1->id()=={id2} && p2->id()=={id1} && p3->id()=={id3}) \|\|
	(p1->id()=={id3} && p2->id()=={id2} && p3->id()=={id1})) {{
	sign *= -1.;
	}}
	norm(norm()*sign);
	"""
	if(not "p1" in couplingptrs[0]) :
	couplingptrs[0] += ' p1'
	if(not "p2" in couplingptrs[1]) :
	couplingptrs[1] += ' p2'
	if(not "p3" in couplingptrs[2]) :
	couplingptrs[2] += ' p3'
	if("Bar" not in vertex.color[0]) :
	order,sign = extractAntiSymmetricIndices(vertex.color[0],"Epsilon(")
	else :
	order,sign = extractAntiSymmetricIndices(vertex.color[0],"EpsilonBar(")
	subs = {"id1" : vertex.particles[int(order[0])-1].pdg_code,
	"id2" : vertex.particles[int(order[1])-1].pdg_code,
	"id3" : vertex.particles[int(order[2])-1].pdg_code,
	"epssign" : sign }
	append+=EPSSIGN.format(**subs)
	return couplingptrs,append

	class VertexConverter:
	'Convert the vertices in a FR model to extract the information ThePEG needs.'
	def __init__(self,model,parmsubs) :
	'Initialize the parameters'
	- self.ONE_EACH=True
	self.verbose=False
	self.vertex_skipped=False
	self.ignore_skipped=False
	self.model=model
	self.all_vertices= []
	self.modelname=""
	- self.globalsign=self.global_sign()
	self.no_generic_loop_vertices = False
	self.parmsubs = parmsubs
	-
	- def global_sign(self):
	- 'Initial pass to find global sign at the moment does nothing'
	- return 1.0
	- # for v in self.model.all_vertices:
	- # pids = sorted([ p.pdg_code for p in v.particles ])
	- # if pids != [-11,11,22]: continue
	- # coup = v.couplings
	- # assert( len(coup) == 1 )
	- # val = coup.values()[0].value
	- # val = evaluate(val)
	- # assert( val.real == 0 )
	- # return 1 if val.imag > 0 else -1

	def readArgs(self,args) :
	'Extract the relevant command line arguments'
	self.ignore_skipped = args.ignore_skipped
	self.verbose = args.verbose
	self.modelname = args.name
	self.no_generic_loop_vertices = args.no_generic_loop_vertices
	-
	+ self.include_generic = args.include_generic
	def should_print(self) :
	'Check if we should output the results'
	return not self.vertex_skipped or self.ignore_skipped

	def convert(self) :
	'Convert the vertices'
	if(self.verbose) :
	print 'verbose mode on: printing all vertices'
	print '-'*60
	labels = ('vertex', 'particles', 'Lorentz', 'C_L', 'C_R', 'norm')
	pprint.pprint(labels)
	- # check if we should merge vertices
	- if(self.ONE_EACH) :
	- self.all_vertices = self.model.all_vertices
	- else:
	- self.all_vertices = collapse_vertices(self.model.all_vertices)
	+ # extract the vertices
	+ self.all_vertices = self.model.all_vertices
	# convert the vertices
	vertexclasses, vertexheaders = [], set()
	for vertexnumber,vertex in enumerate(self.all_vertices,1) :
	# process the vertex
	(skip,vertexClass,vertexHeader) = \
	self.processVertex(vertexnumber,vertex)
	# check it can be handled
	if(skip) : continue
	# add to the list
	vertexclasses.append(vertexClass)
	vertexheaders.add(vertexHeader)
	WRAP = 25
	if vertexnumber % WRAP == 0:
	write_vertex_file({'vertexnumber' : vertexnumber//WRAP,
	'vertexclasses' : '\n'.join(vertexclasses),
	'vertexheaders' : ''.join(vertexheaders),
	'ModelName' : self.modelname})
	vertexclasses = []
	vertexheaders = set()
	# exit if there's vertices we can't handle
	if not self.should_print():
	sys.stderr.write(
	"""
	Error: The conversion was unsuccessful, some vertices could not be
	generated. If you think the missing vertices are not important
	and want to go ahead anyway, use --ignore-skipped.
	Herwig may not give correct results, though.
	"""
	)
	sys.exit(1)
	# if still stuff to output it do it
	if vertexclasses:
	write_vertex_file({'vertexnumber' : vertexnumber//WRAP + 1,
	'vertexclasses' : '\n'.join(vertexclasses),
	'vertexheaders' : ''.join(vertexheaders),
	'ModelName' : self.modelname})

	print '='*60

	def setCouplingPtrs(self,lorentztag,qcd,append,prepend) :
	couplingptrs = [',tcPDPtr']*len(lorentztag)
	if lorentztag == 'VSS':
	couplingptrs[1] += ' p2'
	elif lorentztag == 'FFV':
	couplingptrs[0] += ' p1'
	elif (lorentztag == 'VVV' or lorentztag == 'VVVS' or
	lorentztag == "SSS" or lorentztag == "VVVT" ) \
	and (append or prepend ) :
	couplingptrs[0] += ' p1'
	couplingptrs[1] += ' p2'
	couplingptrs[2] += ' p3'
	elif (lorentztag == 'VVVV' and qcd != 2) or\
	(lorentztag == "SSSS" and prepend ):
	couplingptrs[0] += ' p1'
	couplingptrs[1] += ' p2'
	couplingptrs[2] += ' p3'
	couplingptrs[3] += ' p4'
	return couplingptrs

	def processVertex(self,vertexnumber,vertex) :
	# get the Lorentz tag for the vertex
	lorentztag,order = unique_lorentztag(vertex)
	# check if we should skip the vertex
	vertex.herwig_skip_vertex = checkGhostGoldstoneVertex(lorentztag,vertex)
	if(vertex.herwig_skip_vertex) :
	return (True,"","")
	# get the factor for the vertex
	+ generic = False
	try:
	lf = lfactors[lorentztag]
	except KeyError:
	- msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \
	- 'is not supported.\n'.format(tag=lorentztag, name=vertex.name,
	- ps=' '.join(map(str,vertex.particles)))
	- sys.stderr.write(msg)
	- vertex.herwig_skip_vertex = True
	- self.vertex_skipped=True
	- return (True,"","")
	+ if(not self.include_generic) :
	+ msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \
	+ 'is not supported.\n'.format(tag=lorentztag, name=vertex.name,
	+ ps=' '.join(map(str,vertex.particles)))
	+ sys.stderr.write(msg)
	+ vertex.herwig_skip_vertex = True
	+ self.vertex_skipped=True
	+ return (True,"","")
	+ else :
	+ lf=1.
	+ generic=True
	# get the ids of the particles at the vertex
	- if self.ONE_EACH:
	- plistarray = [ ','.join([ str(vertex.particles[o-1].pdg_code) for o in order ]) ]
	- else:
	- plistarray = [ ','.join([ str(p.pdg_code) for p in pl ])
	- for pl in vertex.particles_list ]
	+ plistarray = [ ','.join([ str(vertex.particles[o-1].pdg_code) for o in order ]) ]
	# parse the colour structure for the vertex
	try:
	L,pos = colors(vertex)
	cf = colorfactor(vertex,L,pos,lorentztag)
	except SkipThisVertex:
	msg = 'Warning: Color structure for vertex ( {ps} ) in {name} ' \
	'is not supported.\n'.format(tag=lorentztag, name=vertex.name,
	ps=' '.join(map(str,vertex.particles)))
	sys.stderr.write(msg)
	vertex.herwig_skip_vertex = True
	self.vertex_skipped=True
	return (True,"","")

	### classname
	classname = 'V_%s' % vertex.name
	+ if(not generic) :
	+ try :
	+ return self.extractGeneric(vertex,order,lorentztag,classname,plistarray,pos,lf,cf)
	+ except SkipThisVertex:
	+ if(not self.include_generic) :
	+ msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \
	+ 'is not supported, may have a non-perturbative form.\n'.format(tag=lorentztag, name=vertex.name,
	+ ps=' '.join(map(str,vertex.particles)))
	+
	+ sys.stderr.write(msg)
	+ vertex.herwig_skip_vertex = True
	+ self.vertex_skipped=True
	+ return (True,"","")
	+ else :
	+ return self.extractGeneral(vertex,order,lorentztag,classname,plistarray,pos,cf)
	+ else :
	+ return self.extractGeneral(vertex,order,lorentztag,classname,plistarray,pos,cf)
	+
	+
	+ def extractGeneric(self,vertex,order,lorentztag,classname,plistarray,pos,lf,cf) :
	# try to extract the couplings
	- try:
	- (all_couplings,header,qcd,qed,prepend,append) = \
	+ (all_couplings,header,qcd,qed,prepend,append) = \
	self.extractCouplings(lorentztag,pos,lf,cf,vertex,order)
	- except SkipThisVertex:
	- msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \
	- 'is not supported, may have a non-perturbative form.\n'.format(tag=lorentztag, name=vertex.name,
	- ps=' '.join(map(str,vertex.particles)))
	- sys.stderr.write(msg)
	- vertex.herwig_skip_vertex = True
	- self.vertex_skipped=True
	- return (True,"","")

	# set the coupling ptrs in the setCoupling call
	couplingptrs = self.setCouplingPtrs(lorentztag,qcd,append != '',prepend != '')

	# final processing of the couplings
	- try :
	- symbols = set()
	- if(lorentztag in ['FFS','FFV']) :
	- (normcontent,leftcontent,rightcontent,append) = processFermionCouplings(lorentztag,vertex,
	- self.model,self.parmsubs,
	- all_couplings)
	- elif('T' in lorentztag) :
	- (leftcontent,rightcontent,normcontent) = processTensorCouplings(lorentztag,vertex,self.model,
	- self.parmsubs,all_couplings)
	- elif(lorentztag=="SSS" or lorentztag=="SSSS") :
	- normcontent = processScalarCouplings(self.model,self.parmsubs,all_couplings)
	- elif(lorentztag=="VVS" or lorentztag =="VVSS" or lorentztag=="VSS") :
	- normcontent,append,lorentztag,header,sym = processScalarVectorCouplings(lorentztag,vertex,
	- self.model,self.parmsubs,
	- all_couplings,header,order)
	- symbols \|=sym
	- elif("VVV" in lorentztag) :
	- normcontent,append,header =\
	- processVectorCouplings(lorentztag,vertex,self.model,self.parmsubs,all_couplings,append,header)
	- else :
	- SkipThisVertex()
	- except SkipThisVertex:
	- msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \
	- 'is not supported, may have a non-perturbative form.\n'.format(tag=lorentztag, name=vertex.name,
	- ps=' '.join(map(str,vertex.particles)))
	- sys.stderr.write(msg)
	- vertex.herwig_skip_vertex = True
	- self.vertex_skipped=True
	- return (True,"","")
	-
	+ symbols = set()
	+ if(lorentztag in ['FFS','FFV']) :
	+ (normcontent,leftcontent,rightcontent,append) = processFermionCouplings(lorentztag,vertex,
	+ self.model,self.parmsubs,
	+ all_couplings)
	+ elif('T' in lorentztag) :
	+ (leftcontent,rightcontent,normcontent) = processTensorCouplings(lorentztag,vertex,self.model,
	+ self.parmsubs,all_couplings)
	+ elif(lorentztag=="SSS" or lorentztag=="SSSS") :
	+ normcontent = processScalarCouplings(self.model,self.parmsubs,all_couplings)
	+ elif(lorentztag=="VVS" or lorentztag =="VVSS" or lorentztag=="VSS") :
	+ normcontent,append,lorentztag,header,sym = processScalarVectorCouplings(lorentztag,vertex,
	+ self.model,self.parmsubs,
	+ all_couplings,header,order)
	+ symbols \|=sym
	+ elif("VVV" in lorentztag) :
	+ normcontent,append,header =\
	+ processVectorCouplings(lorentztag,vertex,self.model,self.parmsubs,all_couplings,append,header)
	+ else :
	+ SkipThisVertex()

	### do we need left/right?
	if 'FF' in lorentztag and lorentztag != "FFT":
	#leftcalc = aStoStrongCoup(py2cpp(leftcontent)[0], paramstoreplace_, paramstoreplace_expressions_)
	#rightcalc = aStoStrongCoup(py2cpp(rightcontent)[0], paramstoreplace_, paramstoreplace_expressions_)
	leftcalc, sym = py2cpp(leftcontent)
	symbols \|= sym
	rightcalc, sym = py2cpp(rightcontent)
	symbols \|= sym
	left = 'left(' + leftcalc + ');'
	right = 'right(' + rightcalc + ');'
	else:
	left = ''
	right = ''
	leftcalc = ''
	rightcalc = ''
	#normcalc = aStoStrongCoup(py2cpp(normcontent)[0], paramstoreplace_, paramstoreplace_expressions_)
	normcalc, sym = py2cpp(normcontent)
	symbols \|= sym
	# UFO is GeV by default (?)
	if lorentztag in ['VVS','SSS']:
	normcalc = 'Complex((%s) * GeV / UnitRemoval::E)' % normcalc
	elif lorentztag in ['GeneralVVS']:
	normcalc = 'Complex(-(%s) * UnitRemoval::E / GeV )' % normcalc
	elif lorentztag in ['FFT','VVT', 'SST', 'FFVT', 'VVVT' , 'VVVS' ]:
	normcalc = 'Complex((%s) / GeV * UnitRemoval::E)' % normcalc
	norm = 'norm(' + normcalc + ');'
	# finally special handling for eps tensors
	if(len(vertex.color)==1 and vertex.color[0].find("Epsilon")>=0) :
	couplingptrs, append = epsilonSign(vertex,couplingptrs,append)
	# define unkown symbols from the model
	symboldefs = [ def_from_model(self.model,s) for s in symbols ]
	### assemble dictionary and fill template
	subs = { 'lorentztag' : lorentztag, # ok
	'classname' : classname, # ok
	'symbolrefs' : '\n '.join(symboldefs),
	'left' : left, # doesn't always exist in base
	'right' : right, # doesn't always exist in base
	'norm' : norm, # needs norm, too

	#################### need survey which different setter methods exist in base classes

	'addToPlist' : '\n'.join([ 'addToList(%s);'%s for s in plistarray]),
	'parameters' : '',
	'setCouplings' : '',
	'qedorder' : qed,
	'qcdorder' : qcd,
	'couplingptrs' : ''.join(couplingptrs),
	'spindirectory' : spindirectory(lorentztag),
	'ModelName' : self.modelname,
	'prepend' : prepend,
	'append' : append,
	'header' : header
	} # ok
	-
	# print info if required
	if self.verbose:
	print '-'*60
	pprint.pprint(( classname, plistarray, leftcalc, rightcalc, normcalc ))

	return (False,VERTEXCLASS.substitute(subs),VERTEXHEADER.format(**subs))

	+ def extractGeneral(self,vertex,order,lorentztag,classname,plistarray,pos,cf) :
	+ defns=[]
	+ vertexEval=[]
	+ values=[]
	+ for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems() :
	+ if(color_idx!=0) :
	+ vertex.herwig_skip_vertex = True
	+ self.vertex_skipped=True
	+ msg = 'Warning: General spin structure code currently only '\
	+ 'supports 1 colour structure for {tag} ( {ps} ) in {name}\n'.format(tag=lorentztag, name=vertex.name,
	+ ps=' '.join(map(str,vertex.particles)))
	+ sys.stderr.write(msg)
	+ return (True,"","")
	+ # calculate the value of the coupling
	+ values.append(couplingValue(coupling))
	+ # now to convert the spin structures
	+ for i in range(0,len(vertex.particles)+1) :
	+ if(len(defns)<i+1) :
	+ defns.append({})
	+ vertexEval.append([])
	+ convertLorentz(vertex.lorentz[lorentz_idx],lorentztag,order,i,defns[i],vertexEval[i])
	+
	+ # we can now generate the evaluate member functions
	+ headers=""
	+ impls=""
	+ imax = len(vertex.particles)+1
	+ if lorentztag in genericVertices :
	+ imax=1
	+ for i in range(0,imax) :
	+ (evalHeader,evalCC) = generateEvaluateFunction(self.model,vertex,i,values,defns[i],vertexEval[i],cf)
	+ headers+=" "+evalHeader+";\n"
	+ impls+=evalCC
	+ impls=impls.replace("evaluate(", "FRModel%s::evaluate(" % classname)
	+ ### assemble dictionary and fill template
	+ subs = { 'lorentztag' : lorentztag,
	+ 'classname' : classname,
	+ 'addToPlist' : '\n'.join([ 'addToList(%s);'%s for s in plistarray]),
	+ 'ModelName' : self.modelname,
	+ 'evaldefs' : headers,
	+ 'evalimpls' : impls}
	+ return (False,GENERALVERTEXCLASS.substitute(subs),GENERALVERTEXHEADER.format(**subs))
	+
	+
	+
	def get_vertices(self,libname):
	vlist = ['library %s\n' % libname]
	for v in self.all_vertices:
	if v.herwig_skip_vertex: continue
	vlist.append( vertexline.format(modelname=self.modelname, vname=v.name) )
	if( not self.no_generic_loop_vertices) :
	vlist.append('insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_GenericHPP\n'.format(modelname=self.modelname) )
	vlist.append('insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_GenericHGG\n'.format(modelname=self.modelname) )
	return ''.join(vlist)


	def extractCouplings(self,lorentztag,pos,lf,cf,vertex,order) :
	coup_left = []
	coup_right = []
	coup_norm = []
	header = ""
	qcd=0
	qed=0
	prepend=""
	append=""
	unique_qcd = CheckUnique()
	unique_qed = CheckUnique()
	maxColour=0
	for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems():
	maxColour=max(maxColour,color_idx)
	all_couplings=[]
	for ix in range(0,maxColour+1) :
	all_couplings.append([])
	for colour in range(0,maxColour+1) :
	for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems() :
	if(color_idx!=colour) : continue
	qcd, qed = qcd_qed_orders(vertex, coupling)
	try :
	unique_qcd( qcd )
	unique_qed( qed )
	except :
	msg = 'Different powers of QCD and QED couplings for the same vertex'\
	' is not currently supported for {ps} in {name}.\n'.format(tag=lorentztag, name=vertex.name,
	ps=' '.join(map(str,vertex.particles)))
	sys.stderr.write(msg)
	raise SkipThisVertex()
	L = vertex.lorentz[lorentz_idx]
	- prefactors = calculatePrefactor(self.globalsign,lorentztag,lf,cf[color_idx])
	+ prefactors = calculatePrefactor(lf,cf[color_idx])
	# calculate the value of the coupling
	value = couplingValue(coupling)
	# handling of the different types of couplings
	if lorentztag in ['FFS','FFV']:
	all_couplings[color_idx] = fermionCouplings(value,prefactors,L,all_couplings[color_idx],order)
	elif 'T' in lorentztag :
	append, all_couplings[color_idx] = tensorCouplings(vertex,value,prefactors,L,lorentztag,pos,
	all_couplings[color_idx],order)
	elif 'R' in lorentztag :
	all_couplings[color_idx] = RSCouplings(value,prefactors,L,all_couplings[color_idx],order)
	elif lorentztag == 'VVS' or lorentztag == "VVSS" or lorentztag == "VSS" :
	all_couplings[color_idx] = scalarVectorCouplings(value,prefactors,L,lorentztag,
	all_couplings[color_idx],order)
	elif lorentztag == "SSS" or lorentztag == "SSSS" :
	prepend, header, all_couplings[color_idx] = scalarCouplings(vertex,value,prefactors,L,lorentztag,
	all_couplings[color_idx],prepend,header)
	elif "VVV" in lorentztag :
	all_couplings[color_idx],append = vectorCouplings(vertex,value,prefactors,L,lorentztag,pos,
	all_couplings[color_idx],append,qcd,order)
	else:
	raise SkipThisVertex()

	# return the result
	return (all_couplings,header,qcd,qed,prepend,append)

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 5:19 PM (14 h, 45 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023613
Default Alt Text: (114 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions