Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7879468
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
15 KB
Subscribers
None
View Options
diff --git a/Models/Feynrules/python/ufo2herwig b/Models/Feynrules/python/ufo2herwig
--- a/Models/Feynrules/python/ufo2herwig
+++ b/Models/Feynrules/python/ufo2herwig
@@ -1,387 +1,399 @@
#! /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(
'--forbidden-particle-name',
action="append",
default=["eta","phi"],
help="Add particle names not allowed as names in UFO models to avoid conflicts with"+\
"Herwig internal particles, names will have _UFO appended"
)
# get the arguments
args = parser.parse_args()
# import the model
import imp
path,mod = os.path.split(os.path.abspath(args.ufodir))
FR = imp.load_module(mod,*imp.find_module(mod,[path]))
##################################################
##################################################
# get the Model name from the arguments
modelname = args.name
libname = modelname + '.so'
# define arrays and variables
#allplist = ""
parmdecls = []
parmgetters = []
parmconstr = []
doinit = []
paramstoreplace_ = []
paramstoreplace_expressions_ = []
# get external parameters for printing
parmsubs = dict( [ (p.name, p.value)
for p in FR.all_parameters
if p.nature == 'external' ] )
# evaluate python cmath
def evaluate(x):
import cmath
return eval(x,
{'cmath':cmath,
'complexconjugate':FR.function_library.complexconjugate},
parmsubs)
## get internal params into arrays
internal = ( p for p in FR.all_parameters
if p.nature == 'internal' )
#paramstoreplaceEW_ = []
#paramstoreplaceEW_expressions_ = []
# calculate internal parameters
for p in internal:
parmsubs.update( { p.name : evaluate(p.value) } )
# if 'aS' in p.value and p.name != 'aS':
# paramstoreplace_.append(p.name)
# paramstoreplace_expressions_.append(p.value)
# if 'aEWM1' in p.value and p.name != 'aEWM1':
# paramstoreplaceEW_.append(p.name)
# paramstoreplaceEW_expressions_.append(p.value)
parmvalues=copy.copy(parmsubs)
# more arrays used for substitution in templates
paramsforstream = []
parmmodelconstr = []
# loop over parameters and fill in template stuff according to internal/external and complex/real
# WARNING: Complex external parameter input not tested!
if args.verbose:
print 'verbose mode on: printing all parameters'
print '-'*60
paramsstuff = ('name', 'expression', 'default value', 'nature')
pprint.pprint(paramsstuff)
interfacedecl_T = """\
static Parameter<{modelname}, {type}> interface{pname}
("{pname}",
"The interface for parameter {pname}",
&{modelname}::{pname}_, {value}, 0, 0,
false, false, Interface::nolimits);
"""
massnames = {}
widthnames = {}
for particle in FR.all_particles:
# skip ghosts and goldstones
if(isGhost(particle) or isGoldstone(particle)) :
continue
- if particle.mass != 'ZERO':
- massnames[particle.mass] = abs(particle.pdg_code)
+ if particle.mass != 'ZERO' and particle.mass.name != 'ZERO':
+ if(particle.mass in massnames) :
+ if(abs(particle.pdg_code) not in massnames[particle.mass]) :
+ massnames[particle.mass].append(abs(particle.pdg_code))
+ else :
+ massnames[particle.mass] = [abs(particle.pdg_code)]
- if particle.width != 'ZERO':
- widthnames[particle.width] = abs(particle.pdg_code)
+ if particle.width != 'ZERO' and particle.width.name != 'ZERO':
+ if(particle.width in widthnames) :
+ if(abs(particle.pdg_code) not in widthnames[particle.width]) :
+ widthnames[particle.width].append(abs(particle.pdg_code))
+ else :
+ 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) )
+ for idCode in massnames[p] :
+ doinit.append(' resetMass(%s,%s_ * GeV);' % (idCode, p.name) )
+ if p in widthnames:
+ for idCode in widthnames[p] :
+ doinit.append(' getParticleData(%s)->width(%s_ * GeV);' % (idCode, p.name) )
+ doinit.append(' getParticleData(%s)->cTau (%s_ == 0.0 ? Length() : hbarc/(%s_*GeV));' % (idCode, p.name, p.name) )
+ doinit.append(' getParticleData(%s)->widthCut(10. * %s_ * GeV);' % (idCode, 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]) )
+ for idCode in massnames[p] :
+ doinit.append(' %s_ = getParticleData(%s)->mass() / GeV;' % (p.name, idCode) )
+ if p in widthnames:
+ for idCode in widthnames[p] :
+ doinit.append(' %s_ = getParticleData(%s)->width() / GeV;' % (p.name, idCode) )
if args.verbose:
pprint.pprint((p.name,p.value, value, p.nature))
pcwriter = ParamCardWriter(FR.all_parameters)
paramcard_output = '\n'.join(pcwriter.output)
### special treatment
# if p.name == 'aS':
# expression = '0.25 * sqr(strongCoupling(q2)) / Constants::pi'
# elif p.name == 'aEWM1':
# expression = '4.0 * Constants::pi / sqr(electroMagneticCoupling(q2))'
# elif p.name == 'Gf':
# expression = 'generator()->standardModel()->fermiConstant() * GeV2'
paramconstructor=': '
for ncount in range(0,len(parmconstr)) :
paramconstructor += parmconstr[ncount]
if(ncount != len(parmconstr) -1) :
paramconstructor += ','
if(ncount%5 == 0 ) :
paramconstructor += "\n"
paramout=""
paramin =""
for ncount in range(0,len(paramsforstream)) :
if(ncount !=0 ) :
paramout += "<< " + paramsforstream[ncount]
paramin += ">> " + paramsforstream[ncount]
else :
paramout += paramsforstream[ncount]
paramin += paramsforstream[ncount]
if(ncount%5 == 0 ) :
paramout += "\n"
paramin += "\n"
parmtextsubs = { 'parmgetters' : '\n'.join(parmgetters),
'parmdecls' : '\n'.join(parmdecls),
'parmconstr' : paramconstructor,
'getters' : '',
'decls' : '',
'addVertex' : '',
'doinit' : '\n'.join(doinit),
'ostream' : paramout,
'istream' : paramin ,
'refs' : '',
'parmextinter': ''.join(interfaceDecls),
'ModelName': modelname,
'calcfunctions': '',
'param_card_data': paramcard_output
}
##################################################
##################################################
##################################################
# set up the conversion of the vertices
vertexConverter = VertexConverter(FR,parmvalues)
vertexConverter.readArgs(args)
# convert the vertices
vertexConverter.convert()
##################################################
##################################################
##################################################
plist, names = thepeg_particles(FR,parmsubs,modelname,modelparameters,args.forbidden_particle_name)
particlelist = [
"# insert HPConstructor:Outgoing 0 /Herwig/{n}/Particles/{p}".format(n=modelname,p=p)
for p in names
]
# make the first one active to have something runnable in the example .in file
particlelist[0] = particlelist[0][2:]
particlelist = '\n'.join(particlelist)
modelfilesubs = { 'plist' : plist,
'vlist' : vertexConverter.get_vertices(libname),
'setcouplings': '\n'.join(parmmodelconstr),
'ModelName': modelname
}
# write the files from templates according to the above subs
if vertexConverter.should_print():
MODEL_HWIN = getTemplate('LHC-FR.in')
if(not args.split_model) :
MODEL_CC = [getTemplate('Model.cc')]
else :
MODEL_EXTRA_CC=getTemplate('Model6.cc')
extra_names=[]
extra_calcs=[]
parmtextsubs['doinit']=""
for i in range(0,len(doinit)) :
if( i%20 == 0 ) :
function_name = "initCalc" +str(int(i/20))
parmtextsubs['doinit'] += function_name +"();\n"
parmtextsubs['calcfunctions'] += "void " + function_name + "();\n"
extra_names.append(function_name)
extra_calcs.append("")
extra_calcs[-1] += doinit[i] + "\n"
for i in range(0,len(extra_names)) :
ccname = '%s_extra_%s.cc' % (modelname,i)
writeFile( ccname, MODEL_EXTRA_CC.substitute({'ModelName' : modelname,
'functionName' : extra_names[i],
'functionCalcs' : extra_calcs[i] }) )
MODEL_CC = [getTemplate('Model1.cc'),getTemplate('Model2.cc'),getTemplate('Model3.cc'),
getTemplate('Model4.cc'),getTemplate('Model5.cc')]
MODEL_H = getTemplate('Model.h')
print 'LENGTH',len(MODEL_CC)
MODELINFILE = getTemplate('FR.model')
writeFile( 'LHC-%s.in' % modelname,
MODEL_HWIN.substitute({ 'ModelName' : modelname,
'Particles' : particlelist })
)
modeltemplate = Template( MODELINFILE.substitute(modelfilesubs) )
writeFile( '%s.h' % modelname, MODEL_H.substitute(parmtextsubs) )
for i in range(0,len(MODEL_CC)) :
if(len(MODEL_CC)==1) :
ccname = '%s.cc' % modelname
else :
ccname = '%s.cc' % (modelname + str(i))
writeFile( ccname, MODEL_CC[i].substitute(parmtextsubs) )
writeFile( modelname +'.template', modeltemplate.template )
writeFile( modelname +'.model', modeltemplate.substitute( modelparameters ) )
# copy the Makefile-FR to current directory,
# replace with the modelname for compilation
with open(os.path.join(modulepath,'Makefile-FR'),'r') as orig:
with open('Makefile','w') as dest:
dest.write(orig.read().replace("FeynrulesModel.so", libname))
print 'finished generating model:\t', modelname
print 'model directory:\t\t', args.ufodir
print 'generated:\t\t\t', len(FR.all_vertices), 'vertices'
print '='*60
print 'library:\t\t\t', libname
print 'input file:\t\t\t', 'LHC-' + modelname +'.in'
print 'model file:\t\t\t', modelname +'.model'
print '='*60
print """\
To complete the installation, compile by typing "make".
An example input file is provided as LHC-FRModel.in,
you'll need to change the required particles in there.
"""
print 'DONE!'
print '='*60
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:29 PM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806040
Default Alt Text
(15 KB)
Attached To
R563 testingHerwigHG
Event Timeline
Log In to Comment