diff --git a/Models/Feynrules/python/ufo2herwig b/Models/Feynrules/python/ufo2herwig
--- a/Models/Feynrules/python/ufo2herwig
+++ b/Models/Feynrules/python/ufo2herwig
@@ -1,473 +1,474 @@
#! /usr/bin/env python
from __future__ import division
from __future__ import print_function
import os, sys, pprint, argparse, re,copy
from string import 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"
)
parser.add_argument(
'--use-Herwig-Higgs',
action="store_true",
help="Use the internal Herwig SM modeling and vertices for Higgs interactions, there may be sign issues but some UFO models have very minimal Higgs interactions"
)
parser.add_argument(
'--use-generic-for-tensors',
action="store_true",
help="Use the generic machinery for all tensor vertices (debugging only)"
)
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"
)
parser.add_argument(
'--convert',
action="store_true",
default=False,
help="Convert the UFO model for python2 to python3 before loading it."
)
# get the arguments
args = parser.parse_args()
-# import the model
-import importlib,importlib.util
+# convert model to python 3 if needed
if(args.convert) :
convertToPython3(args.ufodir)
-# try :
-# path,mod = os.path.split(os.path.abspath(args.ufodir))
-# spec = importlib.util.spec_from_file_location(mod,os.path.join(os.path.abspath(args.ufodir),"__init__.py"))
-# FR = importlib.util.module_from_spec(spec)
-# spec.loader.exec_module(FR)
-# except :
-# print ("Could not load the UFO python module.\n",
-# "This is usually because you are using python3 and the UFO model is in python2.\n",
-# "The --convert option can be used to try and convert it, but there is no guarantee.\n",
-# "As this modifies the model you should make sure you have a backup copy.\n")
-# quit()
-
-
-path,mod = os.path.split(os.path.abspath(args.ufodir))
-spec = importlib.util.spec_from_file_location(mod,os.path.join(os.path.abspath(args.ufodir),"__init__.py"))
-FR = spec.loader.load_module()
+# import the model
+try :
+ import importlib,importlib.util
+ try :
+ path,mod = os.path.split(os.path.abspath(args.ufodir))
+ spec = importlib.util.spec_from_file_location(mod,os.path.join(os.path.abspath(args.ufodir),"__init__.py"))
+ FR = spec.loader.load_module()
+ except :
+ print ("Could not load the UFO python module.\n",
+ "This is usually because you are using python3 and the UFO model is in python2.\n",
+ "The --convert option can be used to try and convert it, but there is no guarantee.\n",
+ "As this modifies the model you should make sure you have a backup copy.\n")
+ quit()
+except :
+ print ("Newer python import did not work, reverting to python 2 imp module")
+ 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,
'im':FR.function_library.im,
're':FR.function_library.re},
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:
for coupling in FR.all_couplings:
for name,value in coupling.order.items():
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' 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' 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':
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.type == 'complex') :
doinit.append(' if(std::isnan(%s_.real()) || std::isnan(%s_.imag()) || std::isinf(%s_.real()) || std::isinf(%s_.imag())) {throw InitException() << "Calculated parameter %s is nan or inf in Feynrules model. Check your input parameters.";}' % (p.name,p.name,p.name,p.name,p.name) )
else :
doinit.append(' if(std::isnan(%s_) || std::isinf(%s_)) {throw InitException() << "Calculated parameter %s is nan or inf in Feynrules model. Check your input parameters.";}' % (p.name,p.name,p.name) )
if p in massnames:
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:
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,couplingDefns)
vertexConverter.readArgs(args)
# convert the vertices
vertexConverter.convert()
cdefs=""
couplingOrders=""
ncount=2
for name,value in couplingDefns.items() :
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,args.use_Herwig_Higgs)
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)