No OneTemporary
Actions

Size

31 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,493 +1,521 @@
	#! /usr/bin/env python
	import os, sys, pprint, argparse, re
	from string import strip

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

	args = parser.parse_args()

	modeldir = args.ufodir.rstrip('/')
	modelpath, module = os.path.split(modeldir)
	if modelpath:
	sys.path.append(os.path.abspath(modelpath))

	FR = __import__(module)

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

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

	# copy the Makefile-FR to current directory,
	# replace with the modelname for compilation
	with open('../Makefile-FR','r') as orig:
	with open('Makefile','w') as dest:
	dest.write(orig.read().replace("FeynrulesModel.so", libname))


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

	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)



	# 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);
	"""

	interfaceDecls = []

	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':
	interfaceDecls.append(
	interfacedecl_T.format(modelname=modelname,
	pname=p.name,
	value=value,
	type=typemap(p.type))
	)
	parmconstr.append('%s_(%s)' % (p.name, value))
	parmmodelconstr.append('set %s:%s %s' % (modelname, 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:
	raise Exception('Unknown data type "%s".' % p.type)
	parmsubs[p.name] = value


	if p.nature == 'external':
	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
	else:
	expression, symbols = py2cpp(p.value)
	text = add_brackets(expression, symbols)
	text=text.replace('()()','()')
	parmgetters.append(' %s %s() const { return %s; }' % (typemap(p.type),p.name, text))

	### 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'
	# elif p.name == 'MZ':
	# expression = 'getParticleData(ThePEG::ParticleID::Z0)->mass() / GeV'


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

	parmtextsubs = { 'parmgetters' : '\n'.join(parmgetters),
	'parmdecls' : '\n'.join(parmdecls),
	'parmconstr' : ': %s' % ','.join(parmconstr),
	'getters' : '',
	'decls' : '',
	'addVertex' : '',
	'ostream' : '<< '.join(paramsforstream),
	'istream' : '>> '.join(paramsforstream),
	'refs' : '',
	'parmextinter': ''.join(interfaceDecls),
	'ModelName': modelname
	}

	print '-'*60

	# write the files from templates according to the above subs
	MODEL_H = getTemplate('Model.h')
	MODEL_CC = getTemplate('Model.cc')
	MODEL_HWIN = getTemplate('LHC-FR.in')
	writeFile( '%s.h' % modelname, MODEL_H.substitute(parmtextsubs) )
	writeFile( '%s.cc' % modelname, MODEL_CC.substitute(parmtextsubs) )
	writeFile( 'LHC-%s.in' % modelname, MODEL_HWIN.substitute(parmtextsubs) )


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


	# get vertex template
	VERTEX = getTemplate('Vertex.cc')
	VERTEXCLASS = getTemplate('Vertex_class')
	VERTEXHEADER = """\
	#include "ThePEG/Helicity/Vertex/{spindirectory}/{lorentztag}Vertex.h"
	"""

	def produce_vertex_file(subs):
	newname = '%s_Vertices_%03d.cc' % (subs['ModelName'],subs['vertexnumber'])
	writeFile( newname, VERTEX.substitute(subs) )

	if args.verbose:
	print 'verbose mode on: printing all vertices'
	print '-'*60
	labels = ('vertex', 'particles', 'Lorentz', 'C_L', 'C_R', 'norm')
	pprint.pprint(labels)

	### initial pass to find global sign
	def global_sign():
	return 1.0
	for v in FR.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

	vertexclasses, vertexheaders = [], set()

	ONE_EACH = True
	if ONE_EACH:
	all_vertices = FR.all_vertices
	else:
	all_vertices = collapse_vertices(FR.all_vertices)

	globalsign = global_sign()
	for vertexnumber,vertex in enumerate(all_vertices,1):

	lorentztag = unique_lorentztag(vertex)

	# remove vertices involving ghost fields
	if 'U' in lorentztag:
	vertex.herwig_skip_vertex = True
	continue
	else:
	vertex.herwig_skip_vertex = False

	lfactors = {
	'FFV' : '-complex(0,1)', # ok
	'VVV' : 'complex(0,1)', # changed to fix ttbar
	'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,1)',
	'VVVT' : '-complex(0,1)',
	'SSSS' : '-complex(0,1)', # ok
	'FFS' : '-complex(0,1)', # ok
	'SST' : 'complex(0,1)',
	'FFT' : 'complex(0,1)',
	'FFVT' : '-complex(0,1)',
	}

	lf = lfactors[lorentztag]

	### Particle ids
	if ONE_EACH:
	plistarray = [ ','.join([ str(p.pdg_code) for p in vertex.particles ]) ]
	else:
	plistarray = [ ','.join([ str(p.pdg_code) for p in pl ])
	for pl in vertex.particles_list ]

	try:
	- cf = colorfactor(vertex)
	+ L,pos = colors(vertex)
	+ cf = colorfactor(vertex,L,pos)
	except SkipThisVertex:
	vertex.herwig_skip_vertex = True
	continue

	### classname
	classname = 'V_%s' % vertex.name

	### parse couplings
	unique_qcd = CheckUnique()
	unique_qed = CheckUnique()

	coup_left = []
	coup_right = []
	coup_norm = []

	if ONE_EACH:
	items = vertex.couplings.iteritems()
	else:
	items = vertex.couplings
	for (color_idx,lorentz_idx),coupling in items:
	# if more than one take QCD over QED and then lowest order in QED
	if (type(coupling) is list) :
	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:
	qed = coupling.order.get('QED',0)
	qcd = coupling.order.get('QCD',0)
	# WARNING: FIX FOR CASES WHEN BOTH ARE ZERO
	# Is there a better way to treat this?
	if qed == qcd == 0: qed = 1
	unique_qcd( qcd )
	unique_qed( qed )
	L = vertex.lorentz[lorentz_idx]
	prefactors = '(%s) * (%s) * (%s)' % (globalsign**(len(lorentztag)-2),lf,cf[color_idx])

	ordering = ''
	if lorentztag in ['FFS','FFV']:
	left,right = parse_lorentz(L.structure)
	if left:
	coup_left.append('(%s) * (%s) * (%s)' % (prefactors,left,coupling.value))
	if right:
	coup_right.append('(%s) * (%s) * (%s)' % (prefactors,right,coupling.value))
	if lorentztag == 'FFV':
	ordering = ('if(p1->id()!=%s) {Complex ltemp=left(), rtemp=right(); left(-rtemp); right(-ltemp);}'
	% vertex.particles[0].pdg_code)
	else:
	if lorentztag == 'VSS':
	if L.structure == 'P(1,3) - P(1,2)':
	prefactors += ' * -1'
	ordering = 'if(p2->id()!=%s){norm(-norm());}'%vertex.particles[1].pdg_code
	+ elif lorentztag == 'VVVV':
	+ if(qcd==2) :
	+ ordering = 'setType(1);\nsetOrder(0,1,2,3);'
	+ # NEED TO THINK MORE ABOUT THIS
	+ else :
	+ ordering = 'setType(2);'
	+ elif lorentztag == 'VVV' :
	+ l = lambda c: len(pos[c])
	+ if l(8)==3 :
	+ pass
	+ else :
	+ 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)
	if(type(coupling) is not list ):
	value = coupling.value
	else :
	value = "("
	for coup in coupling :
	value += '+(%s)' % coup.value
	value +=")"
	coup_norm.append('(%s) * (%s)' % (prefactors,value))


	#print 'Colour :',vertex.color[color_idx]
	#print 'Lorentz %s:'%L.name, L.spins, L.structure
	#print 'Coupling %s:'%C.name, C.value, '\nQED=%s'%qed, 'QCD=%s'%qcd
	#print '---------------'


	leftcontent = ' + '.join(coup_left) if coup_left else '0'
	rightcontent = ' + '.join(coup_right) if coup_right else '0'
	normcontent = ' + '.join(coup_norm) if coup_norm else '1'

	# #print 'Left:',leftcontent
	# #print 'Right:',rightcontent
	# #print 'Norm:',normcontent
	# #print '---------------'


	# ### do we need left/right?
	symbols = set()
	if 'FF' in lorentztag:
	# 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
	if lorentztag in ['VVS','SSS']:
	# UFO is GeV by default (?)
	normcalc = '(%s) * GeV / UnitRemoval::E' % normcalc
	norm = 'norm(' + normcalc + ');'

	# define unkown symbols from the model
	symboldefs = [ def_from_model(FR,s) for s in symbols ]

	couplingptrs = [',tcPDPtr']*len(lorentztag)
	if lorentztag == 'VSS':
	couplingptrs[1] += ' p2'
	elif lorentztag == 'FFV':
	couplingptrs[0] += ' p1'
	+ elif lorentztag == 'VVV':
	+ couplingptrs[0] += ' p1'
	+ couplingptrs[1] += ' p2'
	+ couplingptrs[2] += ' p3'

	### 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' : modelname,
	'ordering' : ordering,
	} # ok

	if args.verbose:
	print '-'*60
	pprint.pprint(( classname, plistarray, leftcalc, rightcalc, normcalc ))

	vertexclasses.append(VERTEXCLASS.substitute(subs))
	vertexheaders.add(VERTEXHEADER.format(**subs))

	WRAP = 25
	if vertexnumber % WRAP == 0:
	produce_vertex_file({'classname':classname,
	'vertexnumber' : vertexnumber//WRAP,
	'vertexclasses' : '\n'.join(vertexclasses),
	'vertexheaders' : ''.join(vertexheaders),
	'ModelName' : modelname})
	vertexclasses = []
	vertexheaders = set()

	if vertexclasses:
	produce_vertex_file({'classname':classname,
	'vertexnumber' : vertexnumber//WRAP + 1,
	'vertexclasses' : '\n'.join(vertexclasses),
	'vertexheaders' : ''.join(vertexheaders),
	'ModelName' : modelname})

	print '='*60

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

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

	def get_vertices():
	vlist = ['library %s\n' % libname]
	for v in all_vertices:
	if v.herwig_skip_vertex: continue
	vlist.append( vertexline.format(modelname=modelname, vname=v.name) )
	return ''.join(vlist)


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


	MODELINFILE = getTemplate('FR.model')

	writeFile( modelname +'.model', MODELINFILE.substitute(modelfilesubs) )

	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",
	copy the generated .so file, into the Herwig++ lib directory
	and the .model and .in files in the directory that you wish to run in.
	"""
	print 'DONE!'
	print '='*60

	diff --git a/Models/Feynrules/python/ufo2peg/__init__.py b/Models/Feynrules/python/ufo2peg/__init__.py
	--- a/Models/Feynrules/python/ufo2peg/__init__.py
	+++ b/Models/Feynrules/python/ufo2peg/__init__.py
	@@ -1,10 +1,10 @@
	from .particles import thepeg_particles

	from .helpers import CheckUnique, getTemplate, writeFile, def_from_model
	-from .helpers import unique_lorentztag, colorfactor, SkipThisVertex
	+from .helpers import unique_lorentztag, colors, colorfactor, SkipThisVertex
	from .helpers import spindirectory, add_brackets, typemap

	from .converter import py2cpp
	from .lorentzparser import parse_lorentz

	from .collapse_vertices import collapse_vertices
	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,353 +1,355 @@
	from string import Template
	from os import path
	import sys

	"""
	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 template from a file."""
	templatename = '{}.template'.format(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 get_lorentztag(spin):
	"""Produce a ThePEG spin tag for the given numeric FR spins."""
	spins = { 1 : 'S', 2 : 'F', 3 : 'V', -1 : 'U', 5 : 'T' }
	result = [ spins[s] for s in spin ]

	def spinsort(a,b):
	"""Helper function for ThePEG's FVST spin tag ordering."""
	if a == b: return 0
	for letter in 'UFVST':
	if a == letter: return -1
	if b == letter: return 1

	result = sorted(result, cmp=spinsort)
	return ''.join(result)

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


	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 {}.".format(lt))
	return spin_directory


	class SkipThisVertex(Exception):
	pass

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

	-def colorfactor(vertex):
	+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):
	def match(patterns):
	result = [ p == t
	for p,t in zip(patterns,vertex.color) ]
	return all(result)

	label = None
	- L = len(struct)
	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({},{})'.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({},{})'.format(*sorted(nums)),)
	if match(label): return ('1',)

	elif l(6) == l(-6) == 2 and L==4:
	sys.stderr.write('Warning: skipping colour structure 6 6 6~ 6~ in %s.\n'
	% vertex)
	raise SkipThisVertex()

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

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

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

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

	elif l(8) == L == 3:
	# if lorentz is FFV get extra minus sign
	lorentztag = unique_lorentztag(vertex)
	factor = '*(-1)' if lorentztag in ['FFV'] else ''
	label = ('f(1,2,3)',)
	if match(label): return ('-complex(0,1)%s'%factor,)
	label = ('f(3,2,1)',)
	if match(label): return ('complex(0,1)%s'%factor,)
	label = ('f(2,1,3)',)
	if match(label): return ('complex(0,1)%s'%factor,)

	elif l(8) == L == 4:
	label = ('f(-1,1,2)*f(3,4,-1)',
	'f(-1,1,3)*f(2,4,-1)',
	'f(-1,1,4)*f(2,3,-1)',
	)
	- if match(label): return ('1','1','1')
	+ if match(label): return ('-1./3.','-1./3.','-1./3.')

	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 ('0.5','0.5')

	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 ('0.5','0.5')

	print vertex
	raise Exception("Unknown colour tag {}.".format(vertex.color))


	def def_from_model(FR,s):
	"""Return a C++ line that defines parameter s as coming from the model file."""
	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:
	result = result.replace(s,s+'()')
	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]

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 6:14 PM (1 d, 19 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805545
Default Alt Text: (31 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions