diff --git a/Models/Feynrules/python/ufo2herwig b/Models/Feynrules/python/ufo2herwig --- a/Models/Feynrules/python/ufo2herwig +++ b/Models/Feynrules/python/ufo2herwig @@ -1,399 +1,400 @@ #! /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}, + 'complexconjugate':FR.function_library.complexconjugate, + 'im':FR.function_library.im}, 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' 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': # # 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: 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) 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 diff --git a/Models/Feynrules/python/ufo2peg/converter.py b/Models/Feynrules/python/ufo2peg/converter.py --- a/Models/Feynrules/python/ufo2peg/converter.py +++ b/Models/Feynrules/python/ufo2peg/converter.py @@ -1,162 +1,165 @@ """ AST visitor class to convert Python expressions into C++ as used by ThePEG """ import ast def py2cpp(expr): """Convert expr to C++ form. Wraps the converter class.""" result = PyToCpp().parse(expr) return result class PyToCppException(Exception): """Base class for all PyToCpp exceptions.""" class PyToCpp(ast.NodeVisitor): """Convert Python math expressions into C++. Returns a tuple (expr,syms): expr -- C++-compatible expression syms -- set of all free variables Usage: >>> expr = '3+2**a*b' >>> PyToCpp().parse(expr) ('(3.0+(pow(2.0,a)*b))', set(['a', 'b'])) Note: The converter is currently not generic, it relies on the conventions of Feynrules' UFO format on the one hand and ThePEG's C++ types on the other. """ def parse(self,expression): """Convert expression to C++ format.""" self.result = [] self.symbols = set() + expression=expression.replace("abs(","cmath.abs(") tree = ast.parse(expression) #print ast.dump(tree) return self.visit(tree) ################################## def visit_Module(self,node): self.generic_visit(node) return ''.join(self.result), self.symbols def generic_visit(self,node): typename = type(node).__name__ harmless = ['Module','Expr'] if typename not in harmless: raise PyToCppException('Missing implementation for %s' % typename) super(PyToCpp,self).generic_visit(node) def visit_UnaryOp(self,node): self.result.append('(') self.visit(node.op) self.visit(node.operand) self.result.append(')') def visit_BinOp(self,node): if type(node.op) == type(ast.Pow()): return self.pow_node(node) self.result.append('(') self.visit(node.left) self.visit(node.op) self.visit(node.right) self.result.append(')') def pow_node(self,node): if is_square(node): self.result.append('sqr(') self.visit(node.left) self.result.append(')') else: self.result.append('pow(') self.visit(node.left) self.result.append(',') self.visit(node.right) self.result.append(')') def visit_Call(self,node): if is_ii(node): self.result.append('ii') else: self.visit(node.func) self.result.append('(') for a in node.args: self.visit(a) self.result.append(',') if self.result[-1] == ',': del self.result[-1] self.result.append(')') def visit_Attribute(self,node): if node.value.id != 'cmath': err = "Don't know how to convert %s module." % node.value.id raise PyToCppException(err) self.result.append(node.attr) def visit_Num(self,node): # some zeros are encoded as 0j if node.n == 0: text = '0.0' else: text = str(float(node.n)) self.result.append(text) def visit_Name(self,node): text = str(node.id) if text == 'complex': text = 'Complex' elif text == 'complexconjugate': text = 'conj' + elif text == 'im': + text = 'imag' elif text not in []: self.symbols.add(text) self.result.append(text) def visit_Mult(self,node): self.result.append('*') def visit_Add(self,node): self.result.append('+') def visit_Sub(self,node): self.result.append('-') def visit_USub(self,node): self.result.append('-') def visit_UAdd(self,node): self.result.append('+') def visit_Div(self,node): self.result.append('/') def visit_Pow(self,node): err = "Shold never get here. BinaryOp catches Pow calls." raise PyToCppException(err) ### Helpers def is_square(node): """Check if a Pow object is just a square.""" try: return node.right.n == 2.0 except: return False def is_ii(node): """Check if a Call object is just the imaginary unit.""" try: return ( node.func.id == 'complex' and node.args[0].n == 0 and node.args[1].n == 1 ) except: return False if __name__ == "__main__": import doctest doctest.testmod()