Page MenuHomeHEPForge
Files F8308619

No OneTemporary
Actions

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,538 +1,565 @@
#! /usr/bin/env python
from __future__ import division
import os, sys, pprint, argparse, re
from string import strip
# 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"
+)
+
+vertex_skipped = False
args = parser.parse_args()
+def should_print():
+ return not vertex_skipped or args.ignore_skipped
+
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(os.path.join(modulepath,'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
+print
-# write the files from templates according to the above subs
-MODEL_H = getTemplate('Model.h')
-MODEL_CC = getTemplate('Model.cc')
-writeFile( '%s.h' % modelname, MODEL_H.substitute(parmtextsubs) )
-writeFile( '%s.cc' % modelname, MODEL_CC.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):
+def write_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)
vertex.herwig_skip_vertex = False
# remove vertices involving ghost fields
if 'U' in lorentztag:
vertex.herwig_skip_vertex = True
continue
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,2)',
'VVVT' : '-complex(0,2)',
'SSSS' : '-complex(0,1)', # ok
'FFS' : '-complex(0,1)', # ok
'SST' : 'complex(0,2)',
'FFT' : '-complex(0,8)',
'FFVT' : '-complex(0,4)',
}
try:
lf = lfactors[lorentztag]
except KeyError:
- msg = 'Warning: Lorentz structure {} ' \
- 'is not supported in {}.\n'.format(lorentztag, vertex.name)
+ msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \
+ 'is not supported.\n'.format(tag=lorentztag, name=vertex.name,
+ ps=' '.join(map(str,vertex.particles)))
sys.stderr.write(msg)
vertex.herwig_skip_vertex = True
+ vertex_skipped = True
continue
### 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:
L,pos = colors(vertex)
cf = colorfactor(vertex,L,pos)
except SkipThisVertex:
vertex.herwig_skip_vertex = True
+ vertex_skipped = 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
try:
for (color_idx,lorentz_idx),coupling in items:
qcd, qed = qcd_qed_orders(vertex, coupling)
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)
elif 'T' in lorentztag :
all_coup, left_coup, right_coup, ordering = \
tensorCouplings(vertex,coupling,prefactors,L,lorentztag,pos)
coup_norm += all_coup
coup_left += left_coup
coup_right += right_coup
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);'
else:
ordering, factor = EWVVVVCouplings(vertex,L)
prefactors += ' * (%s)' % factor
elif lorentztag == 'VVV':
if len(pos[8]) != 3:
ordering = VVVordering(vertex)
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 '---------------'
except SkipThisVertex:
vertex.herwig_skip_vertex = True
+ vertex_skipped = True
continue
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 and lorentztag != "FFT":
# 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
# UFO is GeV by default (?)
if lorentztag in ['VVS','SSS']:
normcalc = '(%s) * GeV / UnitRemoval::E' % normcalc
elif lorentztag in ['FFT','VVT', 'SST', 'FFVT', 'VVVT' ]:
normcalc = 'Complex((%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' or lorentztag == 'VVVT' :
couplingptrs[0] += ' p1'
couplingptrs[1] += ' p2'
couplingptrs[2] += ' p3'
elif lorentztag == 'VVVV' :
couplingptrs[0] += ' p1'
couplingptrs[1] += ' p2'
couplingptrs[2] += ' p3'
couplingptrs[3] += ' p4'
### 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})
+ write_vertex_file({'classname':classname,
+ 'vertexnumber' : vertexnumber//WRAP,
+ 'vertexclasses' : '\n'.join(vertexclasses),
+ 'vertexheaders' : ''.join(vertexheaders),
+ 'ModelName' : modelname})
vertexclasses = []
vertexheaders = set()
+
+
+if not should_print():
+ sys.stderr.write(
+"""
+Error: The conversion was unsuccessful, some vertices could not be
+ generated. If you think the missing vertices are not important
+ and want to go ahead anyway, use --ignore-skipped.
+ Herwig++ may not give correct results, though.
+
+"""
+ )
+
+ sys.exit(1)
if vertexclasses:
- produce_vertex_file({'classname':classname,
- 'vertexnumber' : vertexnumber//WRAP + 1,
- 'vertexclasses' : '\n'.join(vertexclasses),
- 'vertexheaders' : ''.join(vertexheaders),
- 'ModelName' : modelname})
+ write_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)
plist, names = thepeg_particles(FR,parmsubs)
particlelist = [
"# insert HPConstructor:Outgoing 0 /Herwig/{}/Particles/{}".format(modelname,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)
-MODEL_HWIN = getTemplate('LHC-FR.in')
-writeFile( 'LHC-%s.in' % modelname,
- MODEL_HWIN.substitute({ 'ModelName' : modelname,
- 'Particles' : particlelist,
-}))
-
modelfilesubs = { 'plist' : plist,
'vlist' : get_vertices(),
'setcouplings': '\n'.join(parmmodelconstr),
'ModelName': modelname
}
+# write the files from templates according to the above subs
+if should_print():
+ MODEL_HWIN = getTemplate('LHC-FR.in')
+ MODEL_H = getTemplate('Model.h')
+ MODEL_CC = getTemplate('Model.cc')
+ MODELINFILE = getTemplate('FR.model')
+ writeFile( 'LHC-%s.in' % modelname,
+ MODEL_HWIN.substitute({ 'ModelName' : modelname,
+ 'Particles' : particlelist })
+ )
+ writeFile( '%s.h' % modelname, MODEL_H.substitute(parmtextsubs) )
+ writeFile( '%s.cc' % modelname, MODEL_CC.substitute(parmtextsubs) )
+ writeFile( modelname +'.model', MODELINFILE.substitute(modelfilesubs) )
+ # 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))
-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".
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/helpers.py b/Models/Feynrules/python/ufo2peg/helpers.py
--- a/Models/Feynrules/python/ufo2peg/helpers.py
+++ b/Models/Feynrules/python/ufo2peg/helpers.py
@@ -1,589 +1,600 @@
from string import Template
from os import path
import sys,itertools,cmath
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 python string 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 qcd_qed_orders(vertex, coupling):
# 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 + 2 != len(vertex.particles):
qed = len(vertex.particles) - qcd - 2
return qcd, qed
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 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 = 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)
+ sys.stderr.write(
+ 'Warning: Unknown colour structure 6 6 6~ 6~ ( {ps} ) in {name}.\n'
+ .format(name=vertex.name, ps=' '.join(map(str,vertex.particles)))
+ )
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./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')
elif l(8) == 2 and l(8)+l(1)==L :
subs = { 'g1' : pos[8][0], 'g2' : pos[8][1] }
label = ('Identity({g1},{g2})'.format(**subs),)
if match(label) : return ('1.',)
elif l(8) == 3 and l(1)==1 and L==4 :
label = ('f(1,2,3)',)
if match(label): return ('-complex(0,1)',)
label = ('f(3,2,1)',)
if match(label): return ('complex(0,1)',)
label = ('f(2,1,3)',)
if match(label): return ('complex(0,1)',)
- sys.stderr.write("Warning: Unknown colour structure {v.color} in {v.name}.\n"
- .format(v=vertex))
+ sys.stderr.write(
+ "Warning: Unknown colour structure {color} ( {ps} ) in {name}.\n"
+ .format(color = ' '.join(vertex.color), name = vertex.name,
+ ps = ' '.join(map(str,vertex.particles)))
+ )
raise SkipThisVertex()
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:
pattern = r'({})(\W|$)'.format(s)
result = re.sub(pattern, r'\1()\2', result)
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]
# ordering for EW VVV vertices
def VVVordering(vertex) :
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)
return ordering
def tensorCouplings(vertex,coupling,prefactors,L,lorentztag,pos) :
# split the structure into its different terms for analysis
ordering=""
structure1 = L.structure.split()
structures =[]
sign=''
for struct in structure1 :
if(struct=='+') :
continue
elif(struct=='-') :
sign='-'
else :
structures.append(sign+struct.strip())
sign=''
lterms=[]
rterms=[]
if(lorentztag == 'SST') :
terms=[['P(1003,2)','P(2003,1)'],
['P(1003,1)','P(2003,2)'],
['P(-1,1)','P(-1,2)','Metric(1003,2003)']]
signs=[1.,1.,-1.]
elif(lorentztag == 'FFT' ) :
terms=[['P(2003,1)','Gamma(1003,2,1)'],
['P(2003,2)','Gamma(1003,2,1)'],
['P(1003,1)','Gamma(2003,2,1)'],
['P(1003,2)','Gamma(2003,2,1)'],
['P(-1,1)','Gamma(-1,2,1)','Metric(1003,2003)'],
['P(-1,2)','Gamma(-1,2,1)','Metric(1003,2003)']]
signs=[1.,-1.,1.,-1.,-0.5,0.5]
elif(lorentztag == 'VVT' ) :
terms=[['P(-1,1)','P(-1,2)','Metric(1,2003)','Metric(2,1003)'],
['P(-1,1)','P(-1,2)','Metric(1,1003)','Metric(2,2003)'],
['P(-1,1)','P(-1,2)','Metric(1,2)','Metric(1003,2003)'],
['P(1,2)','P(2,1)','Metric(1003,2003)'],
['P(1,2)','P(2003,1)','Metric(2,1003)'],
['P(1,2)','P(1003,1)','Metric(2,2003)'],
['P(2,1)','P(2003,2)','Metric(1,1003)'],
['P(2,1)','P(1003,2)','Metric(1,2003)'],
['P(1003,2)','P(2003,1)','Metric(1,2)'],
['P(1003,1)','P(2003,2)','Metric(1,2)']]
signs=[1.,1.,-1.,1.,-1.,-1.,-1.,-1.,1.,1.]
elif(lorentztag == 'FFVT' ) :
terms = [['Gamma(2004,2,1)','Metric(3,1004)'],
['Gamma(1004,2,1)','Metric(3,2004)'],
['Gamma(3,2,1)','Metric(1004,2004)']]
lterms=[['Gamma(2004,2,-1)','Metric(3,1004)','ProjM(-1,1)'],
['Gamma(1004,2,-1)','Metric(3,2004)','ProjM(-1,1)'],
['Gamma(3,2,-1)','Metric(1004,2004)','ProjM(-1,1)']]
rterms=[['Gamma(2004,2,-1)','Metric(3,1004)','ProjP(-1,1)'],
['Gamma(1004,2,-1)','Metric(3,2004)','ProjP(-1,1)'],
['Gamma(3,2,-1)','Metric(1004,2004)','ProjP(-1,1)']]
signs=[1.,1.,-0.5]
elif(lorentztag == 'VVVT' ) :
# the F(mu nu,rho sigma lambda) terms first
terms = [['P(2004,2)','Metric(1,1004)','Metric(2,3)'],['P(2004,3)','Metric(1,1004)','Metric(2,3)'],
['P(1004,2)','Metric(1,2004)','Metric(2,3)'],['P(1004,3)','Metric(1,2004)','Metric(2,3)'],
['P(2004,3)','Metric(1,3)','Metric(2,1004)'],['P(2004,1)','Metric(1,3)','Metric(2,1004)'],
['P(1004,3)','Metric(1,3)','Metric(2,2004)'],['P(1004,1)','Metric(1,3)','Metric(2,2004)'],
['P(2004,1)','Metric(1,2)','Metric(3,1004)'],['P(2004,2)','Metric(1,2)','Metric(3,1004)'],
['P(1004,1)','Metric(1,2)','Metric(3,2004)'],['P(1004,2)','Metric(1,2)','Metric(3,2004)'],
['P(3,1)','Metric(1,2004)','Metric(2,1004)'],['P(3,2)','Metric(1,2004)','Metric(2,1004)'],
['P(3,1)','Metric(1,1004)','Metric(2,2004)'],['P(3,2)','Metric(1,1004)','Metric(2,2004)'],
['P(3,1)','Metric(1,2)','Metric(1004,2004)'],['P(3,2)','Metric(1,2)','Metric(1004,2004)'],
['P(2,3)','Metric(1,2004)','Metric(3,1004)'],['P(2,1)','Metric(1,2004)','Metric(3,1004)'],
['P(2,3)','Metric(1,1004)','Metric(3,2004)'],['P(2,1)','Metric(1,1004)','Metric(3,2004)'],
['P(2,3)','Metric(1,3)','Metric(1004,2004)'],['P(2,1)','Metric(1,3)','Metric(1004,2004)'],
['P(1,2)','Metric(2,2004)','Metric(3,1004)'],['P(1,3)','Metric(2,2004)','Metric(3,1004)'],
['P(1,2)','Metric(2,1004)','Metric(3,2004)'],['P(1,3)','Metric(2,1004)','Metric(3,2004)'],
['P(1,2)','Metric(2,3)','Metric(1004,2004)'],['P(1,3)','Metric(2,3)','Metric(1004,2004)']]
signs = [1.,-1.,1.,-1.,1.,-1.,1.,-1.,1.,-1.,1.,-1.,
1.,-1.,1.,-1.,-1.,1.,1.,-1.,1.,-1.,-1.,1.,1.,-1.,1.,-1.,-1.,1.]
signs = [1.,-1.,1.,-1.,1.,-1.,1.,-1.,1.,-1.,1.,-1.,
1.,-1.,1.,-1.,-1.,1.,1.,-1.,1.,-1.,-1.,1.,1.,-1.,1.,-1.,-1.,1.]
l = lambda c: len(pos[c])
if l(8)==3 :
pass
else :
ordering = VVVordering(vertex)
# unknown
else :
raise Exception('Unknown data type "%s".' % lorentztag)
sum = [0.,0.,0.]
itype = 0
for types in (terms,lterms,rterms) :
i=0
for term in types:
for perm in itertools.permutations(term):
for j in range(0,len(perm)) :
if(j==0) :
label=perm[j]
else :
label+='*'+perm[j]
for struct in structures :
if label in struct :
reminder = struct.replace(label,'1.',1)
sum[itype] += eval(reminder, {'cmath':cmath} )*signs[i]
i+=1
itype += 1
all_coup = []
left_coup = []
right_coup = []
if(len(lterms)==0) :
all_coup.append('(%s) *(%s) * (%s)' % (sum[0]/float(len(signs)), prefactors,coupling.value))
else :
sum[1] += sum[0]
sum[2] += sum[0]
left_coup .append('(%s) * (%s) * (%s)' % (prefactors,sum[1]/float(len(signs)),coupling.value))
right_coup.append('(%s) * (%s) * (%s)' % (prefactors,sum[2]/float(len(signs)),coupling.value))
return (all_coup,left_coup,right_coup,ordering)
def EWVVVVCouplings(vertex,L) :
terms=['Metric(1,2)*Metric(3,4)',
'Metric(1,3)*Metric(2,4)',
'Metric(1,4)*Metric(2,3)']
structure1 = L.structure.split()
structures =[]
sign=''
for struct in structure1 :
if(struct=='+') :
continue
elif(struct=='-') :
sign='-'
else :
structures.append(sign+struct.strip())
sign=''
factors=[]
for term in terms:
for struct in structures :
if term in struct :
reminder = struct.replace(term,'1.',1)
factors.append(eval(reminder, {'cmath':cmath} ))
if len(factors) != 3:
- sys.stderr.write('Warning: unsupported {} Lorentz structure in {}:\n{}\n'
- .format(unique_lorentztag(vertex), vertex.name, L.structure))
+ sys.stderr.write(
+ 'Warning: unsupported {tag} ( {ps} ) Lorentz structure in {name}:\n{lstr}\n'
+ .format(tag=unique_lorentztag(vertex), name=vertex.name,
+ lstr=L.structure, ps=' '.join(map(str,vertex.particles)))
+ )
raise SkipThisVertex()
factor=0.
order=[]
if(factors[0]==-2.*factors[1] and factors[0]==-2.*factors[2] ) :
order=[0,1,2,3]
factor = factors[0]/2.
elif(factors[1]==-2.*factors[0] and factors[1]==-2.*factors[2] ) :
order=[0,2,1,3]
factor = factors[1]/2.
elif(factors[2]==-2.*factors[0] and factors[2]==-2.*factors[1] ) :
order=[0,3,1,2]
factor = factors[2]/2.
else:
- sys.stderr.write('Warning: unsupported {} Lorentz structure in {}:\n{}\n'
- .format(unique_lorentztag(vertex), vertex.name, L.structure))
+ sys.stderr.write(
+ 'Warning: unsupported {tag} ( {ps} ) Lorentz structure in {name}:\n{lstr}\n'
+ .format(tag=unique_lorentztag(vertex), name=vertex.name,
+ lstr=L.structure, ps=' '.join(map(str,vertex.particles)))
+ )
raise SkipThisVertex()
pattern = \
"bool done[4]={false,false,false,false};\n" + \
" tcPDPtr part[4]={p1,p2,p3,p4};\n" + \
" unsigned int iorder[4]={0,0,0,0};\n" + \
" for(unsigned int ix=0;ix<4;++ix) {\n" + \
" if(!done[0] && part[ix]->id()==%s) {done[0]=true; iorder[%s] = ix; continue;}\n" + \
" if(!done[1] && part[ix]->id()==%s) {done[1]=true; iorder[%s] = ix; continue;}\n" + \
" if(!done[2] && part[ix]->id()==%s) {done[2]=true; iorder[%s] = ix; continue;}\n" + \
" if(!done[3] && part[ix]->id()==%s) {done[3]=true; iorder[%s] = ix; continue;}\n" + \
" }\n" + \
" setType(2);\n" + \
" setOrder(iorder[0],iorder[1],iorder[2],iorder[3]);"
ordering=pattern % ( vertex.particles[0].pdg_code,order[0],
vertex.particles[1].pdg_code,order[1],
vertex.particles[2].pdg_code,order[2],
vertex.particles[3].pdg_code,order[3] )
return (ordering,factor)

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 12:38 PM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022796
Default Alt Text
(44 KB)

Event Timeline