No OneTemporary
Actions

Size

28 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,389 +1,442 @@
	#! /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, float(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.nature == 'external':
	interfaceDecls.append(
	interfacedecl_T.format(modelname=modelname,
	pname=p.name,
	value=value,
	type=typemap(p.type))
	)

	if p.type == 'real':
	assert( value.imag < 1.0e-16 )
	value = value.real
	if p.nature == 'external':
	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':
	parmconstr.append('%s_(%s,%s)' % (p.name, value.real, value.imag))
	parmmodelconstr.append('set %s:%s (%s,%s)' % (modelname, 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)
	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():
	+ 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()
	-all_vertices = collapse_vertices(FR.all_vertices)
	+
	+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)

	+ lfactors = {
	+ 'FFV' : '-complex(0,1)', # ok
	+ 'VVV' : '-complex(0,1)',
	+ 'VVVV' : 'complex(0,1)',
	+ 'VVS' : '-complex(0,1)',
	+ 'VSS' : 'complex(0,1)',
	+ '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
	- plistarray = [ ','.join([ str(p.pdg_code) for p in pl ])
	- for pl in vertex.particles_list ]
	-
	- colortag = color_ok(vertex)
	+ 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 ]
	+
	+ cf = colorfactor(vertex)
	+# print 'Factor:',cf

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

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

	coup_left = []
	coup_right = []
	coup_norm = []
	-
	- for (color_idx,lorentz_idx),coupling in vertex.couplings:
	+
	+ if ONE_EACH:
	+ items = vertex.couplings.iteritems()
	+ else:
	+ items = vertex.couplings
	+ for (color_idx,lorentz_idx),coupling in items:
	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 == 0 and qcd == 0:
	- qed = 1
	+ if qed == qcd == 0: qed = 1
	unique_qcd( qcd )
	unique_qed( qed )
	L = vertex.lorentz[lorentz_idx]
	- C = vertex.color[color_idx]
	- print '-->',C
	+ 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' % (left,coupling.value))
	+ coup_left.append('%s * %s * %s' % (prefactors,left,coupling.value))
	if right:
	- coup_right.append('%s * %s' % (right,coupling.value))
	+ coup_right.append('%s * %s * %s' % (prefactors,right,coupling.value))
	else:
	- coup_norm.append(coupling.value)
	+ 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
	+ coup_norm.append('%s * %s' % (prefactors,coupling.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
	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'
	### 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' : ',tcPDPtr'*len(lorentztag),
	+ '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 'U' in unique_lorentztag(v): 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/Vertex_class.template b/Models/Feynrules/python/ufo2peg/Vertex_class.template
	--- a/Models/Feynrules/python/ufo2peg/Vertex_class.template
	+++ b/Models/Feynrules/python/ufo2peg/Vertex_class.template
	@@ -1,40 +1,41 @@
	class ${ModelName}${classname}: public ${lorentztag}Vertex {
	public:
	${ModelName}${classname}() {
	${addToPlist}
	}
	void setCoupling(Energy2 ${couplingptrs}) {
	${symbolrefs}
	// getParams(q2);
	${norm}
	${left}
	${right}
	+ ${ordering}
	}
	void persistentOutput(PersistentOStream & os) const { os << model_; }
	void persistentInput(PersistentIStream & is, int) { is >> model_; }
	// static void Init();
	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();
	}
	// void getParams(Energy2);
	private:
	${ModelName}${classname} & operator=(const ${ModelName}${classname} &);
	// Complex leftval, rightval, normval;
	tcHw${ModelName}Ptr model_;
	};
	DescribeClass<${ModelName}${classname},Helicity::${lorentztag}Vertex>
	describeHerwig${ModelName}${classname}("Herwig::${ModelName}${classname}",
	"${ModelName}.so");
	// void ${ModelName}${classname}::getParams(Energy2 ) {
	// }
	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,9 +1,10 @@
	from .particles import thepeg_particles

	from .helpers import CheckUnique, getTemplate, writeFile, def_from_model
	-from .helpers import unique_lorentztag, color_ok, spindirectory, add_brackets, typemap
	+from .helpers import unique_lorentztag, colorfactor
	+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,289 +1,309 @@
	from string import Template
	from os import path
	+#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 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 'FVST':
	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 )
	assert( lorentztag == l.name[:len(lorentztag)] )
	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




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

	-def unique_colortag(vertex):
	- unique = CheckUnique()
	- for pl in vertex.particles_list:
	- struct = [ p.color for p in pl ]
	- unique(struct)
	+def colorfactor(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)

	- pos = colorpositions(struct)
	+ 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(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 == 3:
	- label = ('f(3,2,1)',) ###### check if sign as expected
	+ label = ('f(1,2,3)',)
	+ if match(label): return ('-complex(0,1)',)
	+ label = ('f(3,2,1)',)
	+ if match(label): return ('complex(0,1)',)
	+
	+ 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')
	+
	elif l(8) == 2 and l(3) == l(-3) == 1:
	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))
	- return struct, label
	+ if match(label): return ('0.5','0.5')

	-def color_ok(vertex):
	- colortag = unique_colortag(vertex)
	- print
	- print colortag, vertex.color
	- assert( vertex.color == colortag[1] or colortag[1] is None )
	- return colortag[1]
	+ 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:11 PM (1 d, 19 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805538
Default Alt Text: (28 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions