No OneTemporary
Actions

Size

60 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,430 +1,432 @@
	#! /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);
	"""

	# sort out the couplings
	couplingDefns = { "QED" : 99, "QCD" : 99 }
	try :
	for coupling in FR.all_orders:
	name = coupling.name.upper()
	couplingDefns[name]= coupling.expansion_order
	except:
	- pass
	-
	+ for coupling in FR.all_couplings:
	+ for name,value in coupling.order.iteritems():
	+ if(name not in couplingDefns) :
	+ couplingDefns[name]=99

	# sort out the particles

	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,couplingDefns)
	vertexConverter.readArgs(args)
	# convert the vertices
	-couplingDefns = vertexConverter.convert()
	+vertexConverter.convert()
	cdefs=""
	couplingOrders=""
	ncount=2
	for name,value in couplingDefns.iteritems() :
	if(name=="QED") :
	couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" %(name,1,value)
	elif (name=="QCD") :
	couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" %(name,2,value)
	else :
	ncount+=1
	cdefs +=" const T %s = %s;\n" % (name,ncount)
	couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" % (name,ncount,value)
	# coupling definitions
	couplingTemplate= """\
	namespace ThePEG {{
	namespace Helicity {{
	namespace CouplingType {{
	typedef unsigned T;
	/**
	* Enums for couplings
	*/
	{coup}
	}}
	}}
	}}
	"""
	if(cdefs!="") :
	cdefs = couplingTemplate.format(coup=cdefs)
	parmtextsubs['couplings'] = cdefs
	parmtextsubs['couplingOrders'] = couplingOrders

	# particles

	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/helpers.py b/Models/Feynrules/python/ufo2peg/helpers.py
	--- a/Models/Feynrules/python/ufo2peg/helpers.py
	+++ b/Models/Feynrules/python/ufo2peg/helpers.py
	@@ -1,274 +1,266 @@
	from string import Template
	from os import path
	import sys,cmath
	import re

	"""
	Helper functions for the Herwig Feynrules converter
	"""

	class CheckUnique:
	"""Uniqueness checker.

	An object of this class remembers the value it was called with first.
	Any subsequent call to it will only succeed if the same value is passed again.
	For example,

	>>> f = CheckUnique()
	>>> f(5)
	>>> f(5)
	>>> f(4)
	Traceback (most recent call last):
	...
	AssertionError
	"""
	def __init__(self):
	self.val = None

	def __call__(self,val):
	"""Store value on first call, then compare."""
	if self.val is None:
	self.val = val
	else:
	assert( val == self.val )



	def is_number(s):
	"""Check if a value is a number."""
	try:
	float(s)
	except ValueError:
	return False
	return True



	def getTemplate(name):
	"""Create a python string template from a file."""
	templatename = '{name}.template'.format(name=name)
	# assumes the template files sit next to this script
	moduledir = path.dirname(path.abspath(__file__))
	templatepath = path.join(moduledir,templatename)
	with open(templatepath, 'r') as f:
	templateText = f.read()
	return Template( templateText )

	def writeFile(filename, text):
	"""Write text to a filename."""
	with open(filename,'w') as f:
	f.write(text)

	-def coupling_orders(vertex, coupling,model,defns):
	+def coupling_orders(vertex, coupling, defns):
	# if more than one take QCD over QED and then lowest order in QED
	if type(coupling) is list:
	print 'not sure this happens'
	quit()
	qed = 0
	qcd = 0
	for coup in coupling :
	qed1 = coup.order.get('QED',0)
	qcd1 = coup.order.get('QCD',0)
	if qcd1 != 0:
	if qcd == 0 or (qcd1 != 0 and qcd1 < qcd):
	qcd=qcd1
	qed=qed1
	else:
	if qed == 0 or (qed1 != 0 and qed1 < qed):
	qed=qed1
	else:
	output={}
	- try :
	- for ctype in model.all_orders :
	- output[ctype.name]=coupling.order.get(ctype.name,0)
	- total+=output[ctype.name]
	- except :
	- for ctype,value in coupling.order.iteritems() :
	- output[ctype]=value
	- if(ctype not in defns) : defns[ctype]=99
	- if("QED" not in output) : output["QED"]=0
	- if("QCD" not in output) : output["QCD"]=0
	- return output,defns
	+ for ctype in defns :
	+ output[ctype]=coupling.order.get(ctype,0)
	+ return output

	def def_from_model(FR,s):
	"""Return a C++ line that defines parameter s as coming from the model file."""
	if("hw_kine" in s) :return ""
	stype = typemap(getattr(FR.parameters,s).type)
	return '{t} {s} = model_->{s}();'.format(t=stype,s=s)

	_typemap = {'complex':'Complex',
	'real':'double'}

	def typemap(s):
	return _typemap[s]

	def add_brackets(expr, syms):
	result = expr
	for s in syms:
	pattern = r'({symb})(\W\|$)'.format(symb=s)
	result = re.sub(pattern, r'\1()\2', result)
	return result




	def banner():
	return """\
	===============================================================================================================
	______ ______ _ __ _ _ _
	\| ___\| \| ___ \ \| \| / /\| \| \| \| (_) _ _
	\| \|_ ___ _ _ _ __ \| \|_/ /_ _ \| \| ___ ___ / / \| \|_\| \| ___ _ __ __ __ _ __ _ _\| \|_ _\| \|_
	\| _\|/ _ \\| \| \| \|\| \_ \ \| /\| \| \| \|\| \| / _ \/ __\| / / \| _ \| / _ \\| \__\|\ \ /\ / /\| \| / _` \|\|_ _\|\|_ _\|
	\| \| \| __/\| \|_\| \|\| \| \| \|\| \|\ \\| \|_\| \|\| \|\| __/\__ \ / / \| \| \| \|\| __/\| \| \ V V / \| \|\| (_\| \| \|_\| \|_\|
	\_\| \___\| \__, \|\|_\| \|_\|\_\| \_\|\__,_\|\|_\| \___\|\|___//_/ \_\| \|_/ \___\|\|_\| \_/\_/ \|_\| \__, \|
	__/ \| __/ \|
	\|___/ \|___/
	===============================================================================================================
	generating model/vertex/.model/.in files
	please be patient!
	===============================================================================================================
	"""




	#################### ??? #######################


	# function that replaces alphaS (aS)-dependent variables
	# with their explicit form which also contains strongCoupling
	def aStoStrongCoup(stringin, paramstoreplace, paramstoreplace_expressions):
	#print stringin
	for xx in range(0,len(paramstoreplace)):
	#print paramstoreplace[xx], paramstoreplace_expressions[xx]
	stringout = stringin.replace(paramstoreplace[xx], '(' + PyMathToThePEGMath(paramstoreplace_expressions[xx],allparams) + ')')
	stringout = stringout.replace('aS', '(sqr(strongCoupling(q2))/(4.0*Constants::pi))')
	#print 'resulting string', stringout
	return stringout


	# function that replaces alphaEW (aEW)-dependent variables
	# with their explicit form which also contains weakCoupling
	def aEWtoWeakCoup(stringin, paramstoreplace, paramstoreplace_expressions):
	#print stringin
	for xx in range(0,len(paramstoreplace)):
	#print paramstoreplace[xx], paramstoreplace_expressions[xx]
	stringout = stringin.replace(paramstoreplace[xx], '(' + PyMathToThePEGMath(paramstoreplace_expressions[xx],allparams) + ')')
	stringout = stringout.replace('aEWM1', '(1/(sqr(electroMagneticCoupling(q2))/(4.0*Constants::pi)))')
	#print 'resulting string', stringout
	return stringout




	if __name__ == "__main__":
	import doctest
	doctest.testmod()


	if False:

	# Check if the Vertex is self-conjugate or not
	pdgcode = [0,0,0,0]
	notsmvertex = False
	vhasw = 0
	vhasz = 0
	vhasf = 0
	vhasg = 0
	vhash = 0
	vhasp = 0
	# print 'printing particles in vertex'
	for i in range(len(v.particles)):
	# print v.particles[i].pdg_code
	pdgcode[i] = v.particles[i].pdg_code
	if pdgcode[i] == 23:
	vhasz += 1
	if pdgcode[i] == 22:
	vhasp += 1
	if pdgcode[i] == 25:
	vhash += 1
	if pdgcode[i] == 21:
	vhasg += 1
	if pdgcode[i] == 24:
	vhasw += 1
	if abs(pdgcode[i]) < 7 or (abs(pdgcode[i]) > 10 and abs(pdgcode[i]) < 17):
	vhasf += 1
	if pdgcode[i] not in SMPARTICLES:
	notsmvertex = True


	# treat replacement of SM vertices with BSM vertices?
	if notsmvertex == False:
	if( (vhasf == 2 and vhasz == 1) or (vhasf == 2 and vhasw == 1) or (vhasf == 2 and vhash == 1) or (vhasf == 2 and vhasg == 1) or (vhasf == 2 and vhasp == 0) or (vhasg == 3) or (vhasg == 4) or (vhasw == 2 and vhash == 1) or (vhasw == 3) or (vhasw == 4) or (vhash == 1 and vhasg == 2) or (vhash == 1 and vhasp == 2)):
	#print 'VERTEX INCLUDED IN STANDARD MODEL!'
	v.include = 0
	notincluded += 1
	#continue


	selfconjugate = 0
	for j in range(len(pdgcode)):
	for k in range(len(pdgcode)):
	if j != k and j != 0 and abs(pdgcode[j]) == abs(pdgcode[k]):
	selfconjugate = 1
	#print 'self-conjugate vertex'
	# print pdgcode[j]

	# if the Vertex is not self-conjugate, then add the conjugate vertex
	# automatically
	scfac = [1,1,1,1]
	if selfconjugate == 0:
	#first find the self-conjugate particles
	for u in range(len(v.particles)):
	if v.particles[u].selfconjugate == 0:
	scfac[u] = -1
	# print 'particle ', v.particles[u].pdg_code, ' found not to be self-conjugate'

	if selfconjugate == 0:
	plistarray[1] += str(scfac[1] * v.particles[1].pdg_code) + ',' + str(scfac[0] * v.particles[0].pdg_code) + ',' + str(scfac[2] * v.particles[2].pdg_code)
	if len(v.particles) is 4:
	plistarray[1] += ',' + str(scfac[3] * v.particles[3].pdg_code)
	#print 'Conjugate vertex:', plistarray[1]


	class SkipThisVertex(Exception):
	pass

	def extractAntiSymmetricIndices(instring,funct) :
	terms = instring.strip(funct).strip(")").split(",")
	sign=1.
	for iy in range(0,len(terms)) :
	for ix in range(-1,-len(terms)+iy,-1) :
	swap = False
	if(len(terms[ix])==1 and len(terms[ix-1])==1) :
	swap = int(terms[ix])<int(terms[ix-1])
	elif(len(terms[ix])==2 and len(terms[ix-1])==2) :
	swap = int(terms[ix][1])<int(terms[ix-1][1])
	elif(len(terms[ix])==1 and len(terms[ix-1])==2) :
	swap = True
	if(swap) :
	sign *=-1.
	terms[ix],terms[ix-1] = terms[ix-1],terms[ix]
	return (terms,sign)

	def isGoldstone(p) :
	"""check if particle is a Goldstone"""
	def gstest(name):
	try:
	return getattr(p,name)
	except AttributeError:
	return False
	# names of goldstone bosons
	gsnames = ['goldstone','goldstoneboson','GoldstoneBoson']
	if any(map(gstest, gsnames)):
	return True
	return False

	def isGhost(p) :
	"""Check if particle is a ghost"""
	try:
	getattr(p,'GhostNumber')
	except AttributeError:
	return False
	return p.GhostNumber != 0

	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,771 +1,768 @@
	import sys,pprint

	from .helpers import CheckUnique,getTemplate,writeFile,coupling_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,convertLorentz,generateEvaluateFunction,multipleEvaluate
	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','FFSS','FFVS','VVVV']

	skipped5Point=False

	# 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::FRModel{classname} /Herwig/{modelname}/{classname}
	insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/{classname}
	"""

	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 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

	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.',)

	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.',)

	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.',)

	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.',)

	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.',)

	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.',)

	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.',)

	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!

	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!

	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,)

	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)))
	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.')
	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]
	}
	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.)',)
	label = ('f(-1,{g1},{g2})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,)

	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(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,defns) :
	'Initialize the parameters'
	self.verbose=False
	self.vertex_skipped=False
	self.ignore_skipped=False
	self.model=model
	self.all_vertices= []
	self.vertex_names = {}
	self.modelname=""
	self.no_generic_loop_vertices = False
	self.parmsubs = parmsubs
	self.couplingDefns = defns

	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)
	# 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
	- return self.couplingDefns

	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) :
	global skipped5Point
	# 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)
	# check the order of the vertex and skip 5 points
	if(len(lorentztag)>=5) :
	vertex.herwig_skip_vertex = True
	if(not skipped5Point) :
	skipped5Point = True
	print "Skipping 5 point vertices which aren\'t used in Herwig7"

	if(vertex.herwig_skip_vertex) :
	return (True,"","")
	# get the factor for the vertex
	generic = False
	try:
	lf = lfactors[lorentztag]
	except KeyError:
	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
	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 :
	try :
	return self.extractGeneral(vertex,order,lorentztag,classname,plistarray,pos,cf)
	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,"","")
	else :
	try :
	return self.extractGeneral(vertex,order,lorentztag,classname,plistarray,pos,cf)
	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,"","")


	def extractGeneric(self,vertex,order,lorentztag,classname,plistarray,pos,lf,cf) :
	classes=""
	headers=""
	# identify the maximum colour flow and orders of the couplings
	maxColour=0
	couplingOrders=[]
	self.vertex_names[vertex.name] = [classname]
	for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems():
	maxColour=max(maxColour,color_idx)
	- orders, self.couplingDefns = coupling_orders(vertex, coupling, self.model, self.couplingDefns)
	+ orders = coupling_orders(vertex, coupling, self.couplingDefns)
	if(orders not in couplingOrders) : couplingOrders.append(orders)
	# loop the order of the couplings
	iorder = 0
	for corder in couplingOrders :
	iorder +=1
	cname=classname
	if(iorder!=1) :
	cname= "%s_%s" % (classname,iorder)
	self.vertex_names[vertex.name].append(cname)
	header = ""
	prepend=""
	append=""
	all_couplings=[]
	for ix in range(0,maxColour+1) :
	all_couplings.append([])
	# loop over the colour structures
	for colour in range(0,maxColour+1) :
	for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems() :
	# check colour structure and coupling order
	if(color_idx!=colour) : continue
	- orders, self.couplingDefns = coupling_orders(vertex, coupling, self.model, self.couplingDefns)
	- if(orders!=corder) : continue
	+ if(coupling_orders(vertex, coupling, self.couplingDefns)!=corder) : continue
	# get the prefactor for the lorentz structure
	L = vertex.lorentz[lorentz_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,corder["QCD"],order)
	else:
	raise SkipThisVertex()
	# set the coupling ptrs in the setCoupling call
	couplingptrs = self.setCouplingPtrs(lorentztag,corder["QCD"],append != '',prepend != '')
	# final processing of the couplings
	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 ]
	couplingOrder=""
	for coupName,coupVal in corder.iteritems() :
	couplingOrder+=" orderInCoupling(CouplingType::%s,%s);\n" %(coupName,coupVal)
	### assemble dictionary and fill template
	subs = { 'lorentztag' : lorentztag, # ok
	'classname' : cname, # 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
	'addToPlist' : '\n'.join([ 'addToList(%s);'%s for s in plistarray]),
	'parameters' : '',
	'couplingOrders' : couplingOrder,
	'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 ))
	headers+=VERTEXHEADER.format(**subs)
	classes+=VERTEXCLASS.substitute(subs)
	return (False,classes,headers)

	def extractGeneral(self,vertex,order,lorentztag,classname,plistarray,pos,cf) :
	eps=False
	classes=""
	headers=""
	# check the colour flows, two cases supported either 1 flow or 3 in gggg
	cidx=-1
	gluon4point = (len(pos[8])==4 and vertex.lorentz[0].spins.count(3)==4)
	couplingOrders=[]
	for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems() :
	- orders, self.couplingDefns = coupling_orders(vertex, coupling, self.model, self.couplingDefns)
	+ orders = coupling_orders(vertex, coupling, self.couplingDefns)
	if(orders not in couplingOrders) : couplingOrders.append(orders)
	if(gluon4point) :
	color = vertex.color[color_idx]
	f = color.split("*")
	(o1,s1) = extractAntiSymmetricIndices(f[0],"f(")
	(o2,s2) = extractAntiSymmetricIndices(f[1],"f(")
	if(o2[0]<o1[0]) : o1,o2=o2,o1
	color = "f(%s)*f(%s)" % (",".join(o1),",".join(o2))
	label = 'f(1,3,-1)*f(2,4,-1)'
	if(label==color) :
	cidx=color_idx
	else :
	cidx=color_idx
	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,"","")
	if(cidx<0) :
	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)
	# loop over the different orders in the couplings
	iorder=0
	self.vertex_names[vertex.name]=[classname]
	for corder in couplingOrders :
	iorder +=1
	cname=classname
	if(iorder!=1) :
	cname= "%s_%s" % (classname,iorder)
	self.vertex_names[vertex.name].append(cname)
	defns=[]
	vertexEval=[]
	values=[]
	for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems() :
	if(color_idx != cidx) : continue
	- orders, self.couplingDefns = coupling_orders(vertex, coupling, self.model, self.couplingDefns)
	- if(orders!=corder) : continue
	+ if(coupling_orders(vertex, coupling, self.couplingDefns)!=corder) : continue
	# 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([])
	eps \|= convertLorentz(vertex.lorentz[lorentz_idx],lorentztag,order,i,defns[i],vertexEval[i])
	# we can now generate the evaluate member functions
	header=""
	impls=""
	imax = len(vertex.particles)+1
	if lorentztag in genericVertices :
	imax=1
	spins=vertex.lorentz[0].spins
	mult={}
	for i in range(1,6) :
	if( (spins.count(i)>1 and i!=2) or
	(spins.count(i)>2 and i==2) ) : mult[i] = []
	for i in range(0,imax) :
	(evalHeader,evalCC) = generateEvaluateFunction(self.model,vertex,i,values,defns[i],vertexEval[i],cf)
	if(i!=0 and spins[i-1] in mult) :
	if(len(mult[spins[i-1]])==0) : mult[spins[i-1]].append(evalHeader)
	evalHeader=evalHeader.replace("evaluate(","evaluate%s(" % i)
	evalCC =evalCC .replace("evaluate(","evaluate%s(" % i)
	mult[spins[i-1]].append(evalHeader)
	header+=" "+evalHeader+";\n"
	impls+=evalCC
	# combine the multiple defn if needed
	for (key,val) in mult.iteritems() :
	(evalHeader,evalCC) = multipleEvaluate(vertex,key,val)
	if(evalHeader!="") : header += " "+evalHeader+";\n"
	if(evalCC!="") : impls += evalCC
	impls=impls.replace("evaluate", "FRModel%s::evaluate" % cname)
	couplingOrder=""
	for coupName,coupVal in corder.iteritems() :
	couplingOrder+=" orderInCoupling(CouplingType::%s,%s);\n" %(coupName,coupVal)
	### assemble dictionary and fill template
	subs = { 'lorentztag' : lorentztag,
	'classname' : cname,
	'addToPlist' : '\n'.join([ 'addToList(%s);'%s for s in plistarray]),
	'ModelName' : self.modelname,
	'couplingOrders' : couplingOrder,
	'evaldefs' : header,
	'evalimpls' : impls}
	newHeader = GENERALVERTEXHEADER.format(**subs)
	if(eps) : newHeader +="#include \"ThePEG/Helicity/epsilon.h\"\n"
	headers+=newHeader
	classes+=GENERALVERTEXCLASS.substitute(subs)
	return (False,classes,headers)

	def get_vertices(self,libname):
	vlist = ['library %s\n' % libname]
	for v in self.all_vertices:
	if v.herwig_skip_vertex: continue
	for name in self.vertex_names[v.name] :
	vlist.append( vertexline.format(modelname=self.modelname, classname=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)

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:51 PM (1 d, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805243
Default Alt Text: (60 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions