diff --git a/Models/Feynrules/python/ufo2herwig b/Models/Feynrules/python/ufo2herwig --- a/Models/Feynrules/python/ufo2herwig +++ b/Models/Feynrules/python/ufo2herwig @@ -1,449 +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 strip, Template +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 imp +import importlib,importlib.util +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)) -FR = imp.load_module(mod,*imp.find_module(mod,[path])) +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) + ################################################## ################################################## # 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 + 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.iteritems(): + 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.iteritems() : +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) + 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 """\ +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 +""") +print('DONE!') +print('='*60) 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,9 @@ from .particles import thepeg_particles - +from .helpers import convertToPython3 from .helpers import CheckUnique, getTemplate, writeFile, def_from_model from .helpers import add_brackets, typemap,isGhost,isGoldstone from .converter import py2cpp from .write_paramcard import ParamCardWriter from .vertices import VertexConverter diff --git a/Models/Feynrules/python/ufo2peg/check_lorentz.py b/Models/Feynrules/python/ufo2peg/check_lorentz.py --- a/Models/Feynrules/python/ufo2peg/check_lorentz.py +++ b/Models/Feynrules/python/ufo2peg/check_lorentz.py @@ -1,869 +1,870 @@ +from __future__ import print_function import itertools,cmath,re from .helpers import SkipThisVertex,extractAntiSymmetricIndices from .converter import py2cpp from .lorentzparser import parse_lorentz import string,re def compare(a,b) : num=abs(a-b) den=abs(a+b) if(den == 0. and 1e-10) : return True return num/den<1e-10 def evaluate(x,model,parmsubs): import cmath return eval(x, {'cmath':cmath, 'complexconjugate':model.function_library.complexconjugate}, parmsubs) # ordering for EW VVV vertices (ordering not an issue as all same spin) 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,value,prefactors,L,lorentztag,pos,all_couplings,order) : # split the structure into its different terms for analysis ordering="" structures = extractStructures(L) 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)'], ['Metric(1003,2003)']] signs=[1.,1.,-1.,-1.] new_couplings=[False]*len(terms) 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)'], ['Metric(1003,2003)']] signs=[1.,-1.,1.,-1.,-0.5,0.5,1.] new_couplings=[False]*3*len(terms) elif(lorentztag == 'VVT' ) : terms=[['P(-1,1)','P(-1,2)','Metric(1,2003)','Metric(2,1003)'], # from C term ['P(-1,1)','P(-1,2)','Metric(1,1003)','Metric(2,2003)'], # from C term ['P(-1,1)','P(-1,2)','Metric(1,2)','Metric(1003,2003)'], # from C term ['P(1,2)','P(2,1)','Metric(1003,2003)'], # from D term (sym) ['P(1,2)','P(2003,1)','Metric(2,1003)'], # 1st term ['P(1,2)','P(1003,1)','Metric(2,2003)'], # 1st swap ['P(2,1)','P(2003,2)','Metric(1,1003)'], # 2nd term ['P(2,1)','P(1003,2)','Metric(1,2003)'], # 2nd swap ['P(1003,2)','P(2003,1)','Metric(1,2)'], # 3rd term ['P(1003,1)','P(2003,2)','Metric(1,2)'], # 3rd swap ['Metric(1,2003)','Metric(2,1003)'], # from mass term ['Metric(1,1003)','Metric(2,2003)'], # from mass term ['Metric(1,2)','Metric(1003,2003)'], # from mass term ['P(1,1)','P(2,1)','Metric(1003,2003)'], # gauge terms ['P(1,2)','P(2,2)','Metric(1003,2003)'], # gauge terms ['P(1,1)','P(2,2)','Metric(1003,2003)'], # gauge terms ['P(1003,1)','P(1,1)','Metric(2,2003)'], # gauge terms ['P(1003,2)','P(2,2)','Metric(1,2003)'], # gauge terms ['P(2003,1)','P(1,1)','Metric(2,1003)'], # gauge terms ['P(2003,2)','P(2,2)','Metric(1,1003)'], # gauge terms ] signs=[1.,1.,-1.,1.,-1.,-1.,-1.,-1.,1.,1.,1.,1.,-1.,1.,1.,0.25,-1.,-1.,-1.,-1.] new_couplings=[False]*len(terms) 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)'], ['Gamma(2004,2,-1)','Metric(3,1004)'], ['Gamma(1004,2,-1)','Metric(3,2004)'], ['Gamma(3,2,-1)','Metric(1004,2004)']] signs=[1.,1.,-0.5,1.,1.,-0.5] new_couplings=[False]*3*len(terms) 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.] new_couplings=[False]*len(terms) l = lambda c: len(pos[c]) if l(8)!=3 : ordering = VVVordering(vertex) # unknown else : raise Exception('Unknown data type "%s".' % lorentztag) iterm=0 try : for term in terms: for perm in itertools.permutations(term): label = '*'.join(perm) for istruct in range(0,len(structures)) : if label in structures[istruct] : reminder = structures[istruct].replace(label,'1.',1) loc=iterm if(reminder.find("ProjM")>=0) : reminder=re.sub("\*ProjM\(.*,.\)","",reminder) loc+=len(terms) elif(reminder.find("ProjP")>=0) : reminder=re.sub("\*ProjP\(.*,.\)","",reminder) loc+=2*len(terms) structures[istruct] = "Done" val = eval(reminder, {'cmath':cmath} )*signs[iterm] if(new_couplings[loc]) : new_couplings[loc] += val else : new_couplings[loc] = val iterm+=1 except : SkipThisVertex() # check we've handled all the terms for val in structures: if(val!="Done") : raise SkipThisVertex() # special for FFVT if(lorentztag=="FFVT") : t_couplings=new_couplings new_couplings=[False]*9 for i in range(0,9) : - j = i+3*(i/3) - k = i+3+3*(i/3) + j = i+3*(i//3) + k = i+3+3*(i//3) if( not t_couplings[j]) : new_couplings[i] = t_couplings[k] else : new_couplings[i] = t_couplings[j] # set the couplings for icoup in range(0,len(new_couplings)) : if(new_couplings[icoup]) : new_couplings[icoup] = '(%s) * (%s) *(%s)' % (new_couplings[icoup],prefactors,value) if(len(all_couplings)==0) : all_couplings=new_couplings else : for icoup in range(0,len(new_couplings)) : if(new_couplings[icoup] and all_couplings[icoup]) : all_couplings[icoup] = '(%s) + (%s) ' % (new_couplings[icoup],all_couplings[icoup]) elif(new_couplings[icoup]) : all_couplings[icoup] = new_couplings[icoup] # return the results return (ordering,all_couplings) def processTensorCouplings(lorentztag,vertex,model,parmsubs,all_couplings,order) : # check for fermion vertices (i.e. has L/R couplings) fermions = "FF" in lorentztag # test and process the values of the couplings tval = ["Unknown"]*3 value = ["Unknown"]*3 # loop over the colours for icolor in range(0,len(all_couplings)) : lmax = len(all_couplings[icolor]) - if(fermions) : lmax /=3 + if(fermions) : lmax //=3 # loop over the different terms for ix in range(0,lmax) : test = [False]*3 imax=3 # normal case if( not fermions ) : test[0] = all_couplings[icolor][ix] imax=1 else : # first case vector but no L/R couplings if( not all_couplings[icolor][lmax+ix] and not all_couplings[icolor][2*lmax+ix] ) : test[0] = all_couplings[icolor][ix] imax=1 # special for mass terms and massless particles if(not all_couplings[icolor][ix]) : code = abs(vertex.particles[order[0]-1].pdg_code) if(ix==6 and (code ==12 or code ==14 or code==16) ) : continue else : raise SkipThisVertex() # second case L/R couplings elif( not all_couplings[icolor][ix] ) : # L=R, replace with vector if(all_couplings[icolor][lmax+ix] == all_couplings[icolor][2*lmax+ix]) : test[0] = all_couplings[icolor][lmax+ix] else : test[1] = all_couplings[icolor][lmax+ix] test[2] = all_couplings[icolor][2*lmax+ix] else : raise SkipThisVertex() # special handling of mass terms # scalar divide by m**2 if((ix==3 and lorentztag=="SST") or ( ix>=10 and ix<=12 and lorentztag=="VVT" )) : for i in range(0,len(test)) : if(test[i]) : test[i] = '(%s)/%s**2' % (test[i],vertex.particles[order[0]-1].mass.value) # fermions divide by 4*m elif(ix==6 and lorentztag=="FFT" and float(vertex.particles[order[0]-1].mass.value) != 0. ) : for i in range(0,len(test)) : if(test[i]) : test[i] = '-(%s)/%s/4' % (test[i],vertex.particles[order[0]-1].mass.value) # set values on first pass if((tval[0]=="Unknown" and fermions ) or (not fermions and tval[0]=="Unknown" and tval[1]=="Unknown" and tval[2]=="Unknown")) : value = test for i in range(0,len(test)) : if(test[i]) : tval[i] = evaluate(test[i],model,parmsubs) else : for i in range(0,imax) : if(not test[i] and not tval[i]) : continue # special for vector gauge terms if(lorentztag=="VVT" and ix>=13) : continue if(not test[i] or tval[i]=="Unknown") : # special for mass terms and vectors if(lorentztag=="VVT" and ix >=10 and ix <=12 and float(vertex.particles[order[0]-1].mass.value) == 0. ) : continue raise SkipThisVertex() tval2 = evaluate(test[i],model,parmsubs) if(abs(tval[i]-tval2)>1e-6) : # special for fermion mass term if fermions massless if(lorentztag=="FFT" and ix ==6 and tval2 == 0. and float(vertex.particles[order[0]-1].mass.value) == 0. ) : continue raise SkipThisVertex() # simple clean up for i in range(0,len(value)): if(value[i]) : value[i] = value[i].replace("(1.0) * ","").replace(" * (1)","") # put everything together coup_left = 0. coup_right = 0. coup_norm = 0. if(lorentztag == "SST" or lorentztag == "VVT" or lorentztag == "VVVT" or lorentztag == "FFT" ) : coup_norm = value[0] if(value[1] or value[2]) : raise SkipThisVertex() elif(lorentztag=="FFVT") : if(not value[1] and not value[2]) : coup_norm = value[0] coup_left = "1." coup_right = "1." elif(not value[0]) : coup_norm = "1." if(value[1] and value[2]) : coup_left = value[1] coup_right = value[2] elif(value[1]) : coup_left = value[1] coup_right = "0." elif(value[2]) : coup_left = "0." coup_right = value[2] else : raise SkipThisVertex() else : raise SkipThisVertex() else : raise SkipThisVertex() # return the answer return (coup_left,coup_right,coup_norm) def extractStructures(L) : structure1 = L.structure.replace(")-P",") - P").replace(")+P",") + P") structure1 = structure1.split() structures =[] sign='' for struct in structure1 : if(struct=='+') : continue elif(struct=='-') : sign='-' else : structures.append(sign+struct.strip()) sign='' return structures def changeSign(sign1,sign2) : if((sign1=="+" and sign2=="+") or\ (sign1=="-" and sign2=="-")) : return "+" else : return "-" def epsilonOrder(eps) : terms,sign = extractAntiSymmetricIndices(eps,"Epsilon(") return (sign,"Epsilon(%s,%s,%s,%s)" % (terms[0],terms[1],terms[2],terms[3])) def VVSEpsilon(couplings,struct) : if(struct.find("Epsilon")<0) : return fact="" sign="+" if(struct[-1]==")") : fact=struct.split("(")[0] if(fact.find("Epsilon")>=0) : fact="" else : struct=struct.split("(",1)[1][:-1] if(fact[0]=="-") : sign="-" fact=fact[1:] split = struct.split("*") # find the epsilon eps="" for piece in split : if(piece.find("Epsilon")>=0) : eps=piece split.remove(piece) break # and any prefactors for piece in split : if(piece.find("P(")<0) : split.remove(piece) if(piece[0]=="+" or piece[0]=="-") : sign=changeSign(sign,piece[0]) piece=piece[1:] if(fact=="") : fact=piece else : fact = "( %s ) * ( %s )" % (fact , piece) # now sort out the momenta for piece in split : terms=piece.split(",") terms[0]=terms[0].strip("P(") terms[1]=terms[1].strip(")") eps=eps.replace(terms[0],"P%s"%terms[1]) (nsign,eps)=epsilonOrder(eps) if(nsign>0) : sign=changeSign(sign,"-") if(fact=="") : fact="1." if(eps!="Epsilon(1,2,P1,P2)") : return if(couplings[6]==0.) : couplings[6] = "( %s%s )" % (sign,fact) else : couplings[6] = "( %s ) + ( %s%s )" % (couplings[6],sign,fact) def scalarVectorCouplings(value,prefactors,L,lorentztag,all_couplings,order) : # set up the types of term we are looking for if(lorentztag=="VVS") : couplings=[0.,0.,0.,0.,0.,0.,0.] terms=[['P(-1,%s)' % order[0], 'P(-1,%s)' % order[1], 'Metric(%s,%s)' %(order[0],order[1])], ['P(1,%s)' % order[0], 'P(2,%s)' % order[0]], ['P(1,%s)' % order[0], 'P(2,%s)' % order[1]], ['P(1,%s)' % order[1], 'P(2,%s)' % order[0]], ['P(1,%s)' % order[1], 'P(2,%s)' % order[1]], ['Metric(%s,%s)'%(order[0],order[1])]] elif(lorentztag=="VVSS") : couplings=[0.] terms=[['Metric(%s,%s)' % (order[0],order[1])]] elif(lorentztag=="VSS"): couplings=[0.,0.] terms=[['P(%s,%s)' % (order[0],order[2])], ['P(%s,%s)' % (order[0],order[1])]] # extract the lorentz structures structures = extractStructures(L) # handle the scalar couplings itype=-1 try : for term in terms: itype+=1 for perm in itertools.permutations(term): label = '*'.join(perm) for istruct in range(0,len(structures)) : if label in structures[istruct] : reminder = structures[istruct].replace(label,'1.',1) couplings[itype]+=eval(reminder, {'cmath':cmath} ) structures[istruct]='Done' except : raise SkipThisVertex() # special for VVS and epsilon # handle the pseudoscalar couplings for struct in structures : if(struct != "Done" ) : if(lorentztag=="VVS") : VVSEpsilon(couplings,struct) else : raise SkipThisVertex() # put it all together if(len(all_couplings)==0) : for ic in range(0,len(couplings)) : if(couplings[ic]!=0.) : all_couplings.append('(%s) * (%s) * (%s)' % (prefactors,value,couplings[ic])) else : all_couplings.append(False) else : for ic in range(0,len(couplings)) : if(couplings[ic]!=0. and all_couplings[ic]) : all_couplings[ic] = '(%s) * (%s) * (%s) + (%s) ' % (prefactors,value, couplings[ic],all_couplings[ic]) elif(couplings[ic]!=0) : all_couplings[ic] = '(%s) * (%s) * (%s) ' % (prefactors,value,couplings[ic]) return all_couplings def processScalarVectorCouplings(lorentztag,vertex,model,parmsubs,all_couplings,header,order) : # check the values tval = [False]*len(all_couplings[0]) value =[False]*len(all_couplings[0]) for icolor in range(0,len(all_couplings)) : for ix in range(0,len(all_couplings[icolor])) : if(not value[ix]) : value[ix] = all_couplings[icolor][ix] if(value[ix] and not tval[ix]) : tval[ix] = evaluate(value[ix],model,parmsubs) elif(value[ix]) : tval2 = evaluate(all_couplings[icolor][0],model,parmsubs) if(abs(tval[ix]-tval2)>1e-6) : raise SkipThisVertex() append = "" symbols = set() coup_norm=0. if(lorentztag=="VVS") : if(not value[0] and not value[1] and not value[2] and \ not value[3] and not value[4] and not value[6] and value[5]) : coup_norm=value[5] else : for ix in range(0,len(value)) : if(value[ix]) : value[ix], sym = py2cpp(value[ix]) symbols |= sym else : value[ix]=0. lorentztag = 'GeneralVVS' header="kinematics(true);" # g_mu,nv piece of coupling if(value[5]!=0.) : append +='a00( %s + Complex(( %s )* GeV2/invariant(1,2)));\n' % ( value[0],value[5]) else : append +='a00( %s );\n' % value[0] # other couplings append += 'a11( %s );\n a12( %s );\n a21( %s );\n a22( %s );\n aEp( %s );\n' % \ ( value[1],value[2],value[3],value[4],value[6] ) coup_norm="1." elif(lorentztag=="VVSS") : coup_norm = value[0] elif(lorentztag=="VSS") : if(abs(tval[0]+tval[1])>1e-6) : for ix in range(0,len(value)) : if(value[ix]) : value[ix], sym = py2cpp(value[ix]) symbols |= sym else : value[ix]=0. coup_norm = "1." append = 'if(p2->id()==%s) { a( %s ); b( %s);}\n else { a( %s ); b( %s);}' \ % (vertex.particles[order[1]-1].pdg_code, value[0],value[1],value[1],value[0]) else : coup_norm = value[1] append = 'if(p2->id()!=%s){norm(-norm());}' \ % vertex.particles[order[1]-1].pdg_code # return the answer return (coup_norm,append,lorentztag,header,symbols) def getIndices(term) : if(term[0:2]=="P(") : indices = term.strip(")").strip("P(").split(",") mom = int(indices[1]) index = int(indices[0]) return (True,mom,index) else : return (False,0,0) def lorentzScalar(vertex,L) : dotProduct = """(invariant( i[{i1}], i[{i2}] )/GeV2)""" structures=L.structure.split() output="(" for struct in structures: if(struct=="+" or struct=="-") : output+=struct continue structure = struct.split("*") worked = False mom=-1 newTerm="" while True : term = structure[-1] structure.pop() (momentum,mom,index) = getIndices(term) if( not momentum) : break # look for the matching momenta for term in structure : (momentum,mom2,index2) = getIndices(term) if(index2==index) : structure.remove(term) dot = dotProduct.format(i1=mom-1,i2=mom2-1) if(newTerm=="") : newTerm = dot else : newTerm = " ( %s) * ( %s ) " % (newTerm,dot) if(len(structure)==0) : worked = True break if(not worked) : return False else : output+=newTerm output+=")" return output kinematicsline = """\ long id [3]={{{id1},{id2},{id3}}}; long id2[3]={{p1->id(),p2->id(),p3->id()}}; unsigned int i[3]; for(unsigned int ix=0;ix<3;++ix) {{ for(unsigned int iy=0;iy<3;++iy) {{ if(id[ix]==id2[iy]) {{ i[ix] = iy; id2[iy]=0; break; }} }} }} double hw_kine1 = {kine}; """ kinematicsline2 = """\ long id [4]={{{id1},{id2},{id3},{id4}}}; long id2[4]={{p1->id(),p2->id(),p3->id(),p4->id()}}; unsigned int i[4]; for(unsigned int ix=0;ix<4;++ix) {{ for(unsigned int iy=0;iy<4;++iy) {{ if(id[ix]==id2[iy]) {{ i[ix] = iy; id2[iy]=0; break; }} }} }} double hw_kine1 = {kine}; """ kinematicsline3 ="""\ double hw_kine{i} = {kine}; """ def scalarCouplings(vertex,value,prefactors,L,lorentztag, all_couplings,prepend,header) : try : val = int(L.structure) except : output = lorentzScalar(vertex,L) if( not output ) : raise SkipThisVertex() else : if(prepend=="") : if(lorentztag=="SSS") : # order doesn't matter here, all same spin prepend = kinematicsline.format(id1=vertex.particles[0].pdg_code, id2=vertex.particles[1].pdg_code, id3=vertex.particles[2].pdg_code, kine=output) else : # order doesn't matter here, all same spin prepend = kinematicsline2.format(id1=vertex.particles[0].pdg_code, id2=vertex.particles[1].pdg_code, id3=vertex.particles[2].pdg_code, id4=vertex.particles[3].pdg_code, kine=output) value = "(%s) *(hw_kine1)" % value else : osplit=prepend.split("\n") i=-1 while osplit[i]=="": i=i-1 ikin=int(osplit[i].split("=")[0].replace("double hw_kine",""))+1 prepend +=kinematicsline3.format(kine=output,i=ikin) value = "(%s) *(hw_kine%s)" % (value,ikin) header="kinematics(true);" if(len(all_couplings)==0) : all_couplings.append('(%s) * (%s)' % (prefactors,value)) else : all_couplings[0] = '(%s) * (%s) + (%s)' % (prefactors,value,all_couplings[0]) return (prepend, header,all_couplings) def processScalarCouplings(model,parmsubs,all_couplings) : tval = False value = False for icolor in range(0,len(all_couplings)) : if(len(all_couplings[icolor])!=1) : raise SkipThisVertex() if(not value) : value = all_couplings[icolor][0] m = re.findall('hw_kine[0-9]*', all_couplings[icolor][0]) if m: for kine in m: # bizarre number for checks, must be a better option parmsubs[kine] = 987654321. if(not tval) : tval = evaluate(value,model,parmsubs) else : tval2 = evaluate(all_couplings[icolor][0],model,parmsubs) if(abs(tval[i]-tval2)>1e-6) : raise SkipThisVertex() # cleanup and return the answer return value.replace("(1.0) * ","").replace(" * (1)","") def vectorCouplings(vertex,value,prefactors,L,lorentztag,pos, all_couplings,append,qcd,order) : structures=extractStructures(L) terms=[] signs=[] if(lorentztag=="VVV") : terms=[['P(%s,%s)' % (order[2],order[0]),'Metric(%s,%s)' % (order[0],order[1])], ['P(%s,%s)' % (order[2],order[1]),'Metric(%s,%s)' % (order[0],order[1])], ['P(%s,%s)' % (order[1],order[0]),'Metric(%s,%s)' % (order[0],order[2])], ['P(%s,%s)' % (order[1],order[2]),'Metric(%s,%s)' % (order[0],order[2])], ['P(%s,%s)' % (order[0],order[1]),'Metric(%s,%s)' % (order[1],order[2])], ['P(%s,%s)' % (order[0],order[2]),'Metric(%s,%s)' % (order[1],order[2])]] signs=[1.,-1.,-1.,1.,1.,-1.] elif(lorentztag=="VVVV") : terms=[['Metric(%s,%s)' % (order[0],order[3]),'Metric(%s,%s)' % (order[1],order[2])], ['Metric(%s,%s)' % (order[0],order[2]),'Metric(%s,%s)' % (order[1],order[3])], ['Metric(%s,%s)' % (order[0],order[1]),'Metric(%s,%s)' % (order[2],order[3])]] signs=[1.,1.,1.] elif(lorentztag=="VVVS") : terms=[['P(%s,%s)' % (order[2],order[0]),'Metric(%s,%s)' % (order[0],order[1])], ['P(%s,%s)' % (order[2],order[1]),'Metric(%s,%s)' % (order[0],order[1])], ['P(%s,%s)' % (order[1],order[0]),'Metric(%s,%s)' % (order[0],order[2])], ['P(%s,%s)' % (order[1],order[2]),'Metric(%s,%s)' % (order[0],order[2])], ['P(%s,%s)' % (order[0],order[1]),'Metric(%s,%s)' % (order[1],order[2])], ['P(%s,%s)' % (order[0],order[2]),'Metric(%s,%s)' % (order[1],order[2])], ['Epsilon(1,2,3,-1)','P(-1,1)'],['Epsilon(1,2,3,-1)','P(-1,2)'], ['Epsilon(1,2,3,-1)','P(-1,3)']] signs=[1.,-1.,-1.,1.,1.,-1.,1.,1.,1.] # extract the couplings new_couplings = [False]*len(terms) iterm=0 try : for term in terms: for perm in itertools.permutations(term): label = '*'.join(perm) for istruct in range(0,len(structures)) : if label in structures[istruct] : reminder = structures[istruct].replace(label,'1.',1) structures[istruct] = "Done" val = eval(reminder, {'cmath':cmath} )*signs[iterm] if(new_couplings[iterm]) : new_couplings[iterm] += val else : new_couplings[iterm] = val iterm += 1 except : raise SkipThisVertex() # check we've handled all the terms for val in structures: if(val!="Done") : raise SkipThisVertex() # set the couplings for icoup in range(0,len(new_couplings)) : if(new_couplings[icoup]) : new_couplings[icoup] = '(%s) * (%s) *(%s)' % (new_couplings[icoup],prefactors,value) if(len(all_couplings)==0) : all_couplings=new_couplings else : for icoup in range(0,len(new_couplings)) : if(new_couplings[icoup] and all_couplings[icoup]) : all_couplings[icoup] = '(%s) * (%s) *(%s) + (%s) ' % (new_couplings[icoup],prefactors,value,all_couplings[icoup]) elif(new_couplings[icoup]) : all_couplings[icoup] = new_couplings[icoup] # ordering for VVV type vertices if(len(pos[8]) != 3 and (lorentztag=="VVV" or lorentztag=="VVVS")) : append = VVVordering(vertex) return all_couplings,append def processVectorCouplings(lorentztag,vertex,model,parmsubs,all_couplings,append,header) : value = False tval = False if(lorentztag=="VVV") : for icolor in range(0,len(all_couplings)) : # loop over the different terms for ix in range(0,len(all_couplings[icolor])) : if(not value) : value = all_couplings[icolor][ix] tval = evaluate(value,model,parmsubs) else : tval2 = evaluate(all_couplings[icolor][ix],model,parmsubs) if(abs(tval-tval2)>1e-6) : raise SkipThisVertex() elif(lorentztag=="VVVV") : order=[] colours = vertex.color if(len(colours)==1) : tval=[] for i in range(0,3) : tval.append(evaluate(all_couplings[0][i],model,parmsubs)) if(compare(tval[2],-2.*tval[1]) and compare(tval[2],-2.*tval[0]) ) : order=[0,1,2,3] value = "0.5*(%s)" % all_couplings[0][2] elif(compare(tval[1],-2.*tval[2]) and compare(tval[1],-2.*tval[0]) ) : order=[0,2,1,3] value = "0.5*(%s)" % all_couplings[0][1] elif(compare(tval[0],-2.*tval[2]) and compare(tval[0],-2.*tval[1]) ) : order=[0,3,1,2] value = "0.5*(%s)" % all_couplings[0][0] else: 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]);" # order doesn't matter here same spin append = 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] ) else : for icolor in range(0,len(all_couplings)) : col=colours[icolor].split("*") if(len(col)==2 and "f(" in col[0] and "f(" in col[1]) : sign = 1 for i in range(0,2) : col[i],stemp = extractAntiSymmetricIndices(col[i],"f(") for ix in range(0,len(col[i])): col[i][ix]=int(col[i][ix]) sign *=stemp if(col[0][0]>col[1][0]) : col[0],col[1] = col[1],col[0] # first flow if(col[0][0]==1 and col[0][1]==2 and col[1][0] ==3 and col[1][1] == 4) : if(all_couplings[icolor][2] or not all_couplings[icolor][0] or not all_couplings[icolor][1]) : raise SkipThisVertex() if(not value) : value = all_couplings[icolor][1] tval = evaluate(value,model,parmsubs) tval2 = -evaluate(all_couplings[icolor][0],model,parmsubs) tval3 = evaluate(all_couplings[icolor][1],model,parmsubs) elif(col[0][0]==1 and col[0][1]==3 and col[1][0] ==2 and col[1][1] == 4) : if(all_couplings[icolor][1] or not all_couplings[icolor][0] or not all_couplings[icolor][2]) : raise SkipThisVertex() if(not value) : value = all_couplings[icolor][2] tval = evaluate(value,model,parmsubs) tval2 = -evaluate(all_couplings[icolor][0],model,parmsubs) tval3 = evaluate(all_couplings[icolor][2],model,parmsubs) elif(col[0][0]==1 and col[0][1]==4 and col[1][0] ==2 and col[1][1] == 3) : if(all_couplings[icolor][0] or not all_couplings[icolor][1] or not all_couplings[icolor][2]) : raise SkipThisVertex() if(not value) : value = all_couplings[icolor][2] tval = evaluate(value,model,parmsubs) tval2 = -evaluate(all_couplings[icolor][1],model,parmsubs) tval3 = evaluate(all_couplings[icolor][2],model,parmsubs) else : raise SkipThisVertex() if(abs(tval-tval2)>1e-6 or abs(tval-tval3)>1e-6 ) : raise SkipThisVertex() append = 'setType(1);\nsetOrder(0,1,2,3);' else : - print 'unknown colour structure for VVVV vertex' + print('unknown colour structure for VVVV vertex') raise SkipThisVertex() elif(lorentztag=="VVVS") : try : # two distinct cases 0-5 = , 6-8= if(all_couplings[0][0]) : imin=0 imax=6 header="scalar(true);" else : imin=6 imax=9 header="scalar(false);" for icolor in range(0,len(all_couplings)) : # loop over the different terms for ix in range(imin,imax) : if(not value) : value = all_couplings[icolor][ix] tval = evaluate(value,model,parmsubs) else : tval2 = evaluate(value,model,parmsubs) if(abs(tval-tval2)>1e-6) : raise SkipThisVertex() except : SkipThisVertex() # cleanup and return the answer value = value.replace("(1.0) * ","").replace(" * (1)","") return (value,append,header) def fermionCouplings(value,prefactors,L,all_couplings,order) : new_couplings=[False,False] try : new_couplings[0],new_couplings[1] = parse_lorentz(L.structure) except : raise SkipThisVertex() for i in range(0,2) : if new_couplings[i]: new_couplings[i] = '(%s) * (%s) * (%s)' % (prefactors,new_couplings[i],value) if(len(all_couplings)==0) : all_couplings=new_couplings else : for i in range(0,len(new_couplings)) : if(new_couplings[i] and all_couplings[i]) : all_couplings[i] = '(%s) + (%s) ' % (new_couplings[i],all_couplings[i]) elif(new_couplings[i]) : all_couplings[i] = new_couplings[i] return all_couplings def processFermionCouplings(lorentztag,vertex,model,parmsubs,all_couplings,order) : leftcontent = all_couplings[0][0] if all_couplings[0][0] else "0." rightcontent = all_couplings[0][1] if all_couplings[0][1] else "0." tval=[evaluate( leftcontent,model,parmsubs), evaluate(rightcontent,model,parmsubs)] for icolor in range(0,len(all_couplings)) : # loop over the different terms for ix in range(0,len(all_couplings[icolor])) : tval2 = evaluate(all_couplings[icolor][ix],model,parmsubs) if all_couplings[icolor][ix] else 0. if(abs(tval[ix]-tval2)>1e-6) : raise SkipThisVertex() normcontent = "1." append="" if lorentztag == 'FFV': append = ('if(p1->id()!=%s) {Complex ltemp=left(), rtemp=right(); left(-rtemp); right(-ltemp);}' % vertex.particles[order[0]-1].pdg_code) return normcontent,leftcontent,rightcontent,append def RSCouplings(value,prefactors,L,all_couplings,order) : raise SkipThisVertex() diff --git a/Models/Feynrules/python/ufo2peg/collapse_vertices.py b/Models/Feynrules/python/ufo2peg/collapse_vertices.py --- a/Models/Feynrules/python/ufo2peg/collapse_vertices.py +++ b/Models/Feynrules/python/ufo2peg/collapse_vertices.py @@ -1,29 +1,29 @@ class CollapsedVertex(object): def __init__(self,v): self.name = v.name[2:] self.particles_list = [v.particles] self.color = tuple(v.color) self.lorentz = tuple(v.lorentz) - self.couplings = tuple(v.couplings.iteritems()) + self.couplings = tuple(v.couplings.items()) def __hash__(self): return hash((self.color, self.lorentz, self.couplings)) def __eq__(self,other): return ( (self.color, self.lorentz, self.couplings) == (other.color, other.lorentz, other.couplings) ) def collapse_vertices(vs): """Collect particle lists with common structure from the Feynrules vertices.""" vertices = {} for v in vs: cv = CollapsedVertex(v) if cv in vertices: vertices[cv].particles_list += cv.particles_list vertices[cv].name += '_%s' % cv.name else: vertices[cv] = cv return vertices.keys() diff --git a/Models/Feynrules/python/ufo2peg/general_lorentz.py b/Models/Feynrules/python/ufo2peg/general_lorentz.py --- a/Models/Feynrules/python/ufo2peg/general_lorentz.py +++ b/Models/Feynrules/python/ufo2peg/general_lorentz.py @@ -1,3029 +1,3030 @@ +from __future__ import print_function import copy from .helpers import SkipThisVertex,def_from_model from .converter import py2cpp import string,re from string import Template epsValue=[[[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]], [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]], [[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]], [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]], [[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]], [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]], [[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]], [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]]] epsValue[0][1][2][3] = -1. epsValue[0][1][3][2] = 1. epsValue[0][2][1][3] = 1. epsValue[0][2][3][1] = -1. epsValue[0][3][1][2] = -1. epsValue[0][3][2][1] = 1. epsValue[1][0][2][3] = 1. epsValue[1][0][3][2] = -1. epsValue[1][2][0][3] = -1. epsValue[1][2][3][0] = 1. epsValue[1][3][0][2] = 1. epsValue[1][3][2][0] = -1. epsValue[2][0][1][3] = -1. epsValue[2][0][3][1] = 1. epsValue[2][1][0][3] = 1. epsValue[2][1][3][0] = -1. epsValue[2][3][0][1] = -1. epsValue[2][3][1][0] = 1. epsValue[3][0][1][2] = 1. epsValue[3][0][2][1] = -1. epsValue[3][1][0][2] = -1. epsValue[3][1][2][0] = 1. epsValue[3][2][0][1] = 1. epsValue[3][2][1][0] = -1. # self contracted tensor propagator tPropA=[[],[],[],[]] tPropA[0].append(Template("-2. / 3. * (M${iloc}2 + 2 * P${iloc}t ** 2) * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[0].append(Template("-4. / 3. * P${iloc}t * P${iloc}x * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[0].append(Template("-4. / 3. * P${iloc}t * P${iloc}y * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[0].append(Template("-4. / 3. * P${iloc}t * P${iloc}z * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[1].append(Template("-4. / 3. * P${iloc}t * P${iloc}x * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[1].append(Template(" 2. / 3. * (M${iloc}2 - 2 * P${iloc}x ** 2) * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[1].append(Template("-4. / 3. * P${iloc}x * P${iloc}y * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[1].append(Template("-4. / 3. * P${iloc}x * P${iloc}z * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[2].append(Template("-4. / 3. * P${iloc}t * P${iloc}y * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[2].append(Template("-4. / 3. * P${iloc}x * P${iloc}y * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[2].append(Template(" 2. / 3. * (M${iloc}2 - 2 * P${iloc}y ** 2) * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[2].append(Template("-4. / 3. * P${iloc}y * P${iloc}z * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[3].append(Template("-4. / 3. * P${iloc}t * P${iloc}z * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[3].append(Template("-4. / 3. * P${iloc}x * P${iloc}z * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[3].append(Template("-4. / 3. * P${iloc}y * P${iloc}z * (M${iloc}2 -p2)*OM${iloc}**2")) tPropA[3].append(Template(" 2. / 3. * (M${iloc}2 - 2 * P${iloc}z ** 2) * (M${iloc}2 -p2)*OM${iloc}**2")) # tensor propagator 1 index contracted tPropB=[[[],[],[],[]],[[],[],[],[]],[[],[],[],[]],[[],[],[],[]]] tPropB[0][0].append(Template("4. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (1. - OM${iloc} * P${iloc}t ** 2)")) tPropB[0][0].append(Template("-2 * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}x - 2. / 3. * (1. - OM${iloc} * P${iloc}t ** 2) * (${V}x - ${dot}*OM${iloc} * P${iloc}x)")) tPropB[0][0].append(Template(" -2 * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}y - 2. / 3. * (1. - OM${iloc} * P${iloc}t ** 2) * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[0][0].append(Template(" -2 * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}z - 2. / 3. * (1. - OM${iloc} * P${iloc}t ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[0][1].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}x / 3. + (1. - OM${iloc} * P${iloc}t ** 2) * (${V}x - ${dot}*OM${iloc} * P${iloc}x)")) tPropB[0][1].append(Template(" (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1. - OM${iloc} * P${iloc}x ** 2) - OM${iloc} * P${iloc}t * P${iloc}x * (${V}x - ${dot}*OM${iloc} * P${iloc}x) / 3.")) tPropB[0][1].append(Template("-(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}y - OM${iloc} * P${iloc}t * P${iloc}y * (${V}x - ${dot}*OM${iloc} * P${iloc}x) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}x * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[0][1].append(Template("-(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}z - OM${iloc} * P${iloc}t * P${iloc}z * (${V}x - ${dot}*OM${iloc} * P${iloc}x) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}x * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[0][2].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}y / 3. + (1. - OM${iloc} * P${iloc}t ** 2) * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[0][2].append(Template("-(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}y - OM${iloc} * P${iloc}t * P${iloc}x * (${V}y - ${dot}*OM${iloc} * P${iloc}y) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}y * (${V}x - ${dot}*OM${iloc} * P${iloc}x)")) tPropB[0][2].append(Template(" (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1. - OM${iloc} * P${iloc}y ** 2) - OM${iloc} * P${iloc}t * P${iloc}y * (${V}y - ${dot}*OM${iloc} * P${iloc}y) / 3.")) tPropB[0][2].append(Template("-(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}y * P${iloc}z - OM${iloc} * P${iloc}t * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[0][3].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}z / 3. + (1. - OM${iloc} * P${iloc}t ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[0][3].append(Template("-(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}z - OM${iloc} * P${iloc}t * P${iloc}x * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}z * (${V}x - ${dot}*OM${iloc} * P${iloc}x)")) tPropB[0][3].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}y * P${iloc}z - OM${iloc} * P${iloc}t * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[0][3].append(Template("(${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1. - OM${iloc} * P${iloc}z ** 2) - OM${iloc} * P${iloc}t * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) / 3.")) tPropB[1][0].append(Template("-(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}x / 3. + (1 - OM${iloc} * P${iloc}t ** 2) * (${V}x - ${dot}*OM${iloc} * P${iloc}x)")) tPropB[1][0].append(Template(" (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1 - OM${iloc} * P${iloc}x ** 2) - OM${iloc} * P${iloc}t * P${iloc}x * (${V}x - ${dot}*OM${iloc} * P${iloc}x) / 3.")) tPropB[1][0].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}y - OM${iloc} * P${iloc}t * P${iloc}y * (${V}x - ${dot}*OM${iloc} * P${iloc}x) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}x * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[1][0].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}z - OM${iloc} * P${iloc}t * P${iloc}z * (${V}x - ${dot}*OM${iloc} * P${iloc}x) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}x * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[1][1].append(Template(" -2*OM${iloc} * P${iloc}t * P${iloc}x * (${V}x - ${dot}*OM${iloc} * P${iloc}x) - 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropB[1][1].append(Template(" 4. / 3. * (${V}x - ${dot}*OM${iloc} * P${iloc}x) * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropB[1][1].append(Template(" -2 * (${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}x * P${iloc}y - 2. / 3. * (-1 - OM${iloc} * P${iloc}x ** 2) * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[1][1].append(Template(" -2 * (${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}x * P${iloc}z - 2. / 3. * (-1 - OM${iloc} * P${iloc}x ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[1][2].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}y * (${V}x - ${dot}*OM${iloc} * P${iloc}x) - OM${iloc} * P${iloc}t * P${iloc}x * (${V}y - ${dot}*OM${iloc} * P${iloc}y) + 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}y")) tPropB[1][2].append(Template(" -(${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}x * P${iloc}y / 3. + (-1 - OM${iloc} * P${iloc}x ** 2) * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[1][2].append(Template(" (${V}x - ${dot}*OM${iloc} * P${iloc}x) * (-1 - OM${iloc} * P${iloc}y ** 2) - OM${iloc} * P${iloc}x * P${iloc}y * (${V}y - ${dot}*OM${iloc} * P${iloc}y) / 3.")) tPropB[1][2].append(Template("-(${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}y * P${iloc}z - OM${iloc} * P${iloc}x * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y) + 2. / 3.*OM${iloc} * P${iloc}x * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[1][3].append(Template("-OM${iloc} * P${iloc}t * P${iloc}z * (${V}x - ${dot}*OM${iloc} * P${iloc}x) - OM${iloc} * P${iloc}t * P${iloc}x * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}z")) tPropB[1][3].append(Template(" -(${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}x * P${iloc}z / 3. + (-1 - OM${iloc} * P${iloc}x ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[1][3].append(Template(" -(${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}y * P${iloc}z - OM${iloc} * P${iloc}x * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3.*OM${iloc} * P${iloc}x * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[1][3].append(Template("(${V}x - ${dot}*OM${iloc} * P${iloc}x) * (-1 - OM${iloc} * P${iloc}z ** 2) - OM${iloc} * P${iloc}x * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) / 3.")) tPropB[2][0].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}y / 3. + (1 - OM${iloc} * P${iloc}t ** 2) * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[2][0].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}y - OM${iloc} * P${iloc}t * P${iloc}x * (${V}y - ${dot}*OM${iloc} * P${iloc}y) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}y * (${V}x - ${dot}*OM${iloc} * P${iloc}x)")) tPropB[2][0].append(Template(" (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1 - OM${iloc} * P${iloc}y ** 2) - OM${iloc} * P${iloc}t * P${iloc}y * (${V}y - ${dot}*OM${iloc} * P${iloc}y) / 3.")) tPropB[2][0].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}y * P${iloc}z - OM${iloc} * P${iloc}t * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[2][1].append(Template("-OM${iloc} * P${iloc}t * P${iloc}y * (${V}x - ${dot}*OM${iloc} * P${iloc}x) - OM${iloc} * P${iloc}t * P${iloc}x * (${V}y - ${dot}*OM${iloc} * P${iloc}y) + 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}y")) tPropB[2][1].append(Template("-(${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}x * P${iloc}y / 3. + (-1 - OM${iloc} * P${iloc}x ** 2) * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[2][1].append(Template(" (${V}x - ${dot}*OM${iloc} * P${iloc}x) * (-1 - OM${iloc} * P${iloc}y ** 2) - OM${iloc} * P${iloc}x * P${iloc}y * (${V}y - ${dot}*OM${iloc} * P${iloc}y) / 3.")) tPropB[2][1].append(Template("-(${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}y * P${iloc}z - OM${iloc} * P${iloc}x * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y) + 2. / 3.*OM${iloc} * P${iloc}x * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[2][2].append(Template(" -2*OM${iloc} * P${iloc}t * P${iloc}y * (${V}y - ${dot}*OM${iloc} * P${iloc}y) - 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropB[2][2].append(Template(" -2*OM${iloc} * P${iloc}x * P${iloc}y * (${V}y - ${dot}*OM${iloc} * P${iloc}y) - 2. / 3. * (${V}x - ${dot}*OM${iloc} * P${iloc}x) * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropB[2][2].append(Template("4. / 3. * (${V}y - ${dot}*OM${iloc} * P${iloc}y) * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropB[2][2].append(Template(" -2 * (${V}y - ${dot}*OM${iloc} * P${iloc}y)*OM${iloc} * P${iloc}y * P${iloc}z - 2. / 3. * (-1 - OM${iloc} * P${iloc}y ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[2][3].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y) - OM${iloc} * P${iloc}t * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropB[2][3].append(Template(" -OM${iloc} * P${iloc}x * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y) - OM${iloc} * P${iloc}x * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3. * (${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropB[2][3].append(Template(" -(${V}y - ${dot}*OM${iloc} * P${iloc}y)*OM${iloc} * P${iloc}y * P${iloc}z / 3. + (-1 - OM${iloc} * P${iloc}y ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[2][3].append(Template(" (${V}y - ${dot}*OM${iloc} * P${iloc}y) * (-1 - OM${iloc} * P${iloc}z ** 2) - OM${iloc} * P${iloc}y * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) / 3.")) tPropB[3][0].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}t * P${iloc}z / 3. + (1 - OM${iloc} * P${iloc}t ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[3][0].append(Template(" -(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}z - OM${iloc} * P${iloc}t * P${iloc}x * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}z * (${V}x - ${dot}*OM${iloc} * P${iloc}x)")) tPropB[3][0].append(Template("-(${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}y * P${iloc}z - OM${iloc} * P${iloc}t * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3.*OM${iloc} * P${iloc}t * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[3][0].append(Template(" (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1 - OM${iloc} * P${iloc}z ** 2) - OM${iloc} * P${iloc}t * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) / 3.")) tPropB[3][1].append(Template("-OM${iloc} * P${iloc}t * P${iloc}z * (${V}x - ${dot}*OM${iloc} * P${iloc}x) - OM${iloc} * P${iloc}t * P${iloc}x * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}x * P${iloc}z")) tPropB[3][1].append(Template(" -(${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}x * P${iloc}z / 3. + (-1 - OM${iloc} * P${iloc}x ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[3][1].append(Template(" -(${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}y * P${iloc}z - OM${iloc} * P${iloc}x * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3.*OM${iloc} * P${iloc}x * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y)")) tPropB[3][1].append(Template(" (${V}x - ${dot}*OM${iloc} * P${iloc}x) * (-1 - OM${iloc} * P${iloc}z ** 2) - OM${iloc} * P${iloc}x * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) / 3.")) tPropB[3][2].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y) - OM${iloc} * P${iloc}t * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropB[3][2].append(Template("-OM${iloc} * P${iloc}x * P${iloc}z * (${V}y - ${dot}*OM${iloc} * P${iloc}y) - OM${iloc} * P${iloc}x * P${iloc}y * (${V}z - ${dot}*OM${iloc} * P${iloc}z) + 2. / 3. * (${V}x - ${dot}*OM${iloc} * P${iloc}x)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropB[3][2].append(Template(" -(${V}y - ${dot}*OM${iloc} * P${iloc}y)*OM${iloc} * P${iloc}y * P${iloc}z / 3. + (-1 - OM${iloc} * P${iloc}y ** 2) * (${V}z - ${dot}*OM${iloc} * P${iloc}z)")) tPropB[3][2].append(Template(" (${V}y - ${dot}*OM${iloc} * P${iloc}y) * (-1 - OM${iloc} * P${iloc}z ** 2) - OM${iloc} * P${iloc}y * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) / 3.")) tPropB[3][3].append(Template(" -2*OM${iloc} * P${iloc}t * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) - 2. / 3. * (${V}t - ${dot}*OM${iloc} * P${iloc}t) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropB[3][3].append(Template("-2*OM${iloc} * P${iloc}x * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) - 2. / 3. * (${V}x - ${dot}*OM${iloc} * P${iloc}x) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropB[3][3].append(Template("-2*OM${iloc} * P${iloc}y * P${iloc}z * (${V}z - ${dot}*OM${iloc} * P${iloc}z) - 2. / 3. * (${V}y - ${dot}*OM${iloc} * P${iloc}y) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropB[3][3].append(Template("4. / 3. * (${V}z - ${dot}*OM${iloc} * P${iloc}z) * (-1 - OM${iloc} * P${iloc}z ** 2)")) # tensor propagator, no contracted indices tPropC=[[[[],[],[],[]],[[],[],[],[]],[[],[],[],[]],[[],[],[],[]]], [[[],[],[],[]],[[],[],[],[]],[[],[],[],[]],[[],[],[],[]]], [[[],[],[],[]],[[],[],[],[]],[[],[],[],[]],[[],[],[],[]]], [[[],[],[],[]],[[],[],[],[]],[[],[],[],[]],[[],[],[],[]]]] tPropC[0][0][0].append(Template("4./3. * (1 - OM${iloc} * P${iloc}t ** 2) ** 2")) tPropC[0][0][0].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}x")) tPropC[0][0][0].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}y")) tPropC[0][0][0].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}z")) tPropC[0][0][1].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}x")) tPropC[0][0][1].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x ** 2 - 2./3. * (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[0][0][1].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y")) tPropC[0][0][1].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z")) tPropC[0][0][2].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}y")) tPropC[0][0][2].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y")) tPropC[0][0][2].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y ** 2 - 2./3. * (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[0][0][2].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[0][0][3].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}z")) tPropC[0][0][3].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z")) tPropC[0][0][3].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[0][0][3].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}z ** 2 - 2./3. * (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[0][1][0].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}x")) tPropC[0][1][0].append(Template("(1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}x ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x ** 2 /3.")) tPropC[0][1][0].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y /3.")) tPropC[0][1][0].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z /3.")) tPropC[0][1][1].append(Template("(1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}x ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x ** 2 /3.")) tPropC[0][1][1].append(Template("-4./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[0][1][1].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y /3. - OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[0][1][1].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[0][1][2].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y /3.")) tPropC[0][1][2].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y /3. - OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[0][1][2].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x + 2./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[0][1][2].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[0][1][3].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z /3.")) tPropC[0][1][3].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[0][1][3].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[0][1][3].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x + 2./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[0][2][0].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}y")) tPropC[0][2][0].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y /3.")) tPropC[0][2][0].append(Template("(1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y ** 2 /3.")) tPropC[0][2][0].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[0][2][1].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y /3.")) tPropC[0][2][1].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[0][2][1].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x /3.")) tPropC[0][2][1].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[0][2][2].append(Template("(1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y ** 2 /3.")) tPropC[0][2][2].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x /3.")) tPropC[0][2][2].append(Template("-4./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[0][2][2].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[0][2][3].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[0][2][3].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[0][2][3].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[0][2][3].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[0][3][0].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}z")) tPropC[0][3][0].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z /3.")) tPropC[0][3][0].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[0][3][0].append(Template("(1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}z ** 2 /3.")) tPropC[0][3][1].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z /3.")) tPropC[0][3][1].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[0][3][1].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[0][3][1].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x /3.")) tPropC[0][3][2].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[0][3][2].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[0][3][2].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[0][3][2].append(Template("-OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[0][3][3].append(Template("(1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}z ** 2 /3.")) tPropC[0][3][3].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x /3.")) tPropC[0][3][3].append(Template("-OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[0][3][3].append(Template("-4./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[1][0][0].append(Template(" -4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}x")) tPropC[1][0][0].append(Template(" (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}x ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x ** 2 /3.")) tPropC[1][0][0].append(Template(" -(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y /3.")) tPropC[1][0][0].append(Template(" -(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z /3.")) tPropC[1][0][1].append(Template(" (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}x ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x ** 2 /3.")) tPropC[1][0][1].append(Template(" -4./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][0][1].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y /3. - OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][0][1].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][0][2].append(Template(" -(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y /3.")) tPropC[1][0][2].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y /3. - OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][0][2].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x + 2./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[1][0][2].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[1][0][3].append(Template(" -(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z /3.")) tPropC[1][0][3].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][0][3].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[1][0][3].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x + 2./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[1][1][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x ** 2 - 2./3. * (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][1][0].append(Template(" -4./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][1][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][1][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][1][1].append(Template(" -4./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][1][1].append(Template(" 4./3. * (-1 - OM${iloc} * P${iloc}x ** 2) ** 2")) tPropC[1][1][1].append(Template(" -4./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}x * P${iloc}y")) tPropC[1][1][1].append(Template(" -4./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}x * P${iloc}z")) tPropC[1][1][2].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][1][2].append(Template(" -4./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}x * P${iloc}y")) tPropC[1][1][2].append(Template(" 2*OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y ** 2 - 2./3. * (-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[1][1][2].append(Template(" 2*OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z + 2./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[1][1][3].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][1][3].append(Template(" -4./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}x * P${iloc}z")) tPropC[1][1][3].append(Template(" 2*OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z + 2./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[1][1][3].append(Template(" 2*OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}z ** 2 - 2./3. * (-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[1][2][0].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y")) tPropC[1][2][0].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y /3. - OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][2][0].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x /3.")) tPropC[1][2][0].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[1][2][1].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y /3. - OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][2][1].append(Template(" -4./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}x * P${iloc}y ")) tPropC[1][2][1].append(Template(" (-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y ** 2 /3. ")) tPropC[1][2][1].append(Template(" -(-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[1][2][2].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x /3.")) tPropC[1][2][2].append(Template(" (-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y ** 2 /3. ")) tPropC[1][2][2].append(Template(" -4./3.*OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[1][2][2].append(Template(" OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[1][2][3].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[1][2][3].append(Template(" -(-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[1][2][3].append(Template(" OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[1][2][3].append(Template(" 2*OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[1][3][0].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z")) tPropC[1][3][0].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][3][0].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[1][3][0].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x /3.")) tPropC[1][3][1].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[1][3][1].append(Template("-4./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}x * P${iloc}z ")) tPropC[1][3][1].append(Template("-(-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[1][3][1].append(Template("(-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}z ** 2 /3. ")) tPropC[1][3][2].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[1][3][2].append(Template("-(-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[1][3][2].append(Template("2*OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[1][3][2].append(Template("-OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[1][3][3].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x /3.")) tPropC[1][3][3].append(Template("(-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}z ** 2 /3. ")) tPropC[1][3][3].append(Template("-OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[1][3][3].append(Template("-4./3.*OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[2][0][0].append(Template("-4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}y ")) tPropC[2][0][0].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y /3. ")) tPropC[2][0][0].append(Template("(1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y ** 2 /3. ")) tPropC[2][0][0].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z /3. ")) tPropC[2][0][1].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y /3. ")) tPropC[2][0][1].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[2][0][1].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x /3.")) tPropC[2][0][1].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[2][0][2].append(Template("(1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y ** 2 /3. ")) tPropC[2][0][2].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x /3.")) tPropC[2][0][2].append(Template("-4./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][0][2].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[2][0][3].append(Template("-(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z /3. ")) tPropC[2][0][3].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[2][0][3].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[2][0][3].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[2][1][0].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}y + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}y")) tPropC[2][1][0].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y /3. - OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[2][1][0].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x /3.")) tPropC[2][1][0].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[2][1][1].append(Template("OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}y /3. - OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[2][1][1].append(Template("-4./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}x * P${iloc}y ")) tPropC[2][1][1].append(Template("(-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y ** 2 /3. ")) tPropC[2][1][1].append(Template("-(-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[2][1][2].append(Template("-OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x /3.")) tPropC[2][1][2].append(Template("(-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y ** 2 /3. ")) tPropC[2][1][2].append(Template("-4./3.*OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][1][2].append(Template("OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[2][1][3].append(Template("4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[2][1][3].append(Template("-(-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[2][1][3].append(Template("OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[2][1][3].append(Template("2*OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[2][2][0].append(Template("2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y ** 2 - 2./3. * (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][0].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x + 2./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][0].append(Template("-4./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][0].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][1].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}x + 2./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][1].append(Template("2*OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y ** 2 - 2./3. * (-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][1].append(Template("-4./3.*OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][1].append(Template("2*OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][2].append(Template("-4./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][2].append(Template("-4./3.*OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][2].append(Template("4./3. * (-1 - OM${iloc} * P${iloc}y ** 2) ** 2 ")) tPropC[2][2][2].append(Template("-4./3. * (-1 - OM${iloc} * P${iloc}y ** 2)*OM${iloc} * P${iloc}y * P${iloc}z ")) tPropC[2][2][3].append(Template("2*OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][3].append(Template("2*OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2) ")) tPropC[2][2][3].append(Template("-4./3. * (-1 - OM${iloc} * P${iloc}y ** 2)*OM${iloc} * P${iloc}y * P${iloc}z ")) tPropC[2][2][3].append(Template("2*OM${iloc}**2 * P${iloc}y ** 2 * P${iloc}z ** 2 - 2./3. * (-1 - OM${iloc} * P${iloc}y ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[2][3][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[2][3][0].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[2][3][0].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[2][3][0].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[2][3][1].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[2][3][1].append(Template(" 2*OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z + 2./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[2][3][1].append(Template(" OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[2][3][1].append(Template(" -OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[2][3][2].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[2][3][2].append(Template(" OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[2][3][2].append(Template(" -4./3. * (-1 - OM${iloc} * P${iloc}y ** 2)*OM${iloc} * P${iloc}y * P${iloc}z ")) tPropC[2][3][2].append(Template(" (-1 - OM${iloc} * P${iloc}y ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}y ** 2 * P${iloc}z ** 2 /3. ")) tPropC[2][3][3].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[2][3][3].append(Template(" -OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[2][3][3].append(Template(" (-1 - OM${iloc} * P${iloc}y ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}y ** 2 * P${iloc}z ** 2 /3. ")) tPropC[2][3][3].append(Template(" -4./3.*OM${iloc} * P${iloc}y * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2) ")) tPropC[3][0][0].append(Template(" -4./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}t * P${iloc}z ")) tPropC[3][0][0].append(Template(" -(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z /3. ")) tPropC[3][0][0].append(Template(" -(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z /3. ")) tPropC[3][0][0].append(Template(" (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}z ** 2 /3. ")) tPropC[3][0][1].append(Template(" -(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z /3. ")) tPropC[3][0][1].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[3][0][1].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[3][0][1].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x /3.")) tPropC[3][0][2].append(Template(" -(1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z /3. ")) tPropC[3][0][2].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[3][0][2].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[3][0][2].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[3][0][3].append(Template(" (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}z ** 2 /3. ")) tPropC[3][0][3].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x /3.")) tPropC[3][0][3].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[3][0][3].append(Template(" -4./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][1][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}x * P${iloc}z + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}x * P${iloc}z")) tPropC[3][1][0].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[3][1][0].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[3][1][0].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x /3.")) tPropC[3][1][1].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}x ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}x ** 2)")) tPropC[3][1][1].append(Template(" -4./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}x * P${iloc}z ")) tPropC[3][1][1].append(Template(" -(-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[3][1][1].append(Template(" (-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}z ** 2 /3. ")) tPropC[3][1][2].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z ")) tPropC[3][1][2].append(Template(" -(-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z /3.")) tPropC[3][1][2].append(Template(" 2*OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z + 2./3.*OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[3][1][2].append(Template(" -OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[3][1][3].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x /3.")) tPropC[3][1][3].append(Template(" (-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}z ** 2 /3. ")) tPropC[3][1][3].append(Template(" -OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[3][1][3].append(Template(" -4./3.*OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2) ")) tPropC[3][2][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}y * P${iloc}z + 2./3. * (1 - OM${iloc} * P${iloc}t ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[3][2][0].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[3][2][0].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[3][2][0].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[3][2][1].append(Template(" 4./3.*OM${iloc}**2 * P${iloc}t * P${iloc}y * P${iloc}x * P${iloc}z")) tPropC[3][2][1].append(Template(" 2*OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}y * P${iloc}z + 2./3. * (-1 - OM${iloc} * P${iloc}x ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[3][2][1].append(Template(" OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[3][2][1].append(Template(" -OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[3][2][2].append(Template(" OM${iloc}**2 * P${iloc}t * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[3][2][2].append(Template(" OM${iloc}**2 * P${iloc}x * P${iloc}y ** 2 * P${iloc}z /3. - OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}y ** 2)")) tPropC[3][2][2].append(Template(" -4./3. * (-1 - OM${iloc} * P${iloc}y ** 2)*OM${iloc} * P${iloc}y * P${iloc}z")) tPropC[3][2][2].append(Template(" (-1 - OM${iloc} * P${iloc}y ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}y ** 2 * P${iloc}z ** 2 /3. ")) tPropC[3][2][3].append(Template(" -OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[3][2][3].append(Template(" -OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y /3.")) tPropC[3][2][3].append(Template(" (-1 - OM${iloc} * P${iloc}y ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2) + OM${iloc}**2 * P${iloc}y ** 2 * P${iloc}z ** 2 /3. ")) tPropC[3][2][3].append(Template(" -4./3.*OM${iloc} * P${iloc}y * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t ** 2 * P${iloc}z ** 2 - 2./3. * (1 - OM${iloc} * P${iloc}t ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x + 2./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][0].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][0].append(Template(" -4./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2) ")) tPropC[3][3][1].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}x + 2./3.*OM${iloc} * P${iloc}t * P${iloc}x * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][1].append(Template(" 2*OM${iloc}**2 * P${iloc}x ** 2 * P${iloc}z ** 2 - 2./3. * (-1 - OM${iloc} * P${iloc}x ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][1].append(Template(" 2*OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][1].append(Template(" -4./3.*OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2) ")) tPropC[3][3][2].append(Template(" 2*OM${iloc}**2 * P${iloc}t * P${iloc}z ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}t * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][2].append(Template(" 2*OM${iloc}**2 * P${iloc}x * P${iloc}z ** 2 * P${iloc}y + 2./3.*OM${iloc} * P${iloc}x * P${iloc}y * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][2].append(Template(" 2*OM${iloc}**2 * P${iloc}y ** 2 * P${iloc}z ** 2 - 2./3. * (-1 - OM${iloc} * P${iloc}y ** 2) * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][2].append(Template(" -4./3.*OM${iloc} * P${iloc}y * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][3].append(Template(" -4./3.*OM${iloc} * P${iloc}t * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][3].append(Template(" -4./3.*OM${iloc} * P${iloc}x * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][3].append(Template(" -4./3.*OM${iloc} * P${iloc}y * P${iloc}z * (-1 - OM${iloc} * P${iloc}z ** 2)")) tPropC[3][3][3].append(Template(" 4./3. * (-1 - OM${iloc} * P${iloc}z ** 2) ** 2")) imap=["t","x","y","z"] RSDotProduct = Template("${s}ts${si}*${v}t-${s}xs${si}*${v}x-${s}ys${si}*${v}y-${s}zs${si}*${v}z") vTemplateT="""\ {header} {{ if({type}W{iloc}.id()=={id}) {{ return {normal}; }} else {{ return {transpose}; }} }}; """ vTemplate4="""\ {header} {{ {swap} if(id{iloc1}=={id1}) {{ if(id{iloc2}=={id2}) {{ return {res1}; }} else {{ return {res2}; }} }} else {{ if(id{iloc2}=={id2}) {{ return {res3}; }} else {{ return {res4}; }} }} }}; """ vecTemplate="""\ Energy2 p2 = P{iloc}.m2(); LorentzPolarizationVector vtemp = {res}; Complex fact = -Complex(0.,1.)*({cf})*propagator(iopt,p2,out,mass,width); if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); complex mass2 = sqr(mass); if(mass.real()==ZERO) {{ vtemp =fact*vtemp; }} else {{ complex dot = P{iloc}*vtemp; vtemp = fact*(vtemp-dot/mass2*P{iloc}); }} return VectorWaveFunction(P{iloc},out,vtemp.x(),vtemp.y(),vtemp.z(),vtemp.t()); """ sTemplate="""\ Energy2 p2 = P{iloc}.m2(); if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); Complex fact = Complex(0.,1.)*({cf})*propagator(iopt,p2,out,mass,width); Lorentz{offTypeA} newSpin = fact*({res}); return {offTypeB}(P{iloc},out,newSpin); """ RSTemplate="""\ if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); Energy2 p2 = P{iloc}.m2(); Complex fact = Complex(0.,-1.)*({cf})*propagator(iopt,p2,out,mass,width); complex Omass = mass.real()==ZERO ? InvEnergy(ZERO) : 1./mass; Lorentz{offTypeA} newSpin = fact*({res}); return {offTypeB}(P{iloc},out,newSpin); """ scaTemplate="""\ if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); Energy2 p2 = P{iloc}.m2(); Complex fact = Complex(0.,1.)*({cf})*propagator(iopt,p2,out,mass,width); complex output = fact*({res}); return ScalarWaveFunction(P{iloc},out,output); """ tenTemplate="""\ if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); InvEnergy2 OM{iloc} = mass.real()==ZERO ? InvEnergy2(ZERO) : 1./sqr(mass.real()); Energy2 M{iloc}2 = sqr(mass.real()); Energy2 p2 = P{iloc}.m2(); Complex fact = Complex(0.,1.)*({cf})*propagator(iopt,p2,out,mass,width); LorentzTensor output = fact*({res}); return TensorWaveFunction(P{iloc},out,output); """ # various strings for matrixes I4 = "Matrix([[1.,0,0,0],[0,1.,0,0],[0,0,1.,0],[0,0,0,1.]])" G5 = "Matrix([[-1.,0,0,0],[0,-1.,0,0],[0,0,1.,0],[0,0,0,1.]])" PM = "Matrix([[1.,0,0,0],[0,1.,0,0],[0,0,0,0],[0,0,0,0]])" PP = "Matrix([[0,0,0,0],[0,0,0,0],[0,0,1.,0],[0,0,0,1.]])" vslash = Template("Matrix([[0,0,${v}TMZ,-${v}XMY],[0,0,-${v}XPY,${v}TPZ],[${v}TPZ,${v}XMY,0,0],[${v}XPY,${v}TMZ,0,0]])") vslashS = Template("${v}TPZ=Symbol(\"${v}TPZ\")\n${v}TMZ=Symbol(\"${v}TMZ\")\n${v}XPY=Symbol(\"${v}XPY\")\n${v}XMY=Symbol(\"${v}XMY\")\n") momCom = Template("${v}t = Symbol(\"${v}t\")\n${v}x = Symbol(\"${v}x\")\n${v}y = Symbol(\"${v}y\")\n${v}z = Symbol(\"${v}z\")\n") vslashD = Template("complex<${var}> ${v}TPZ = ${v}.t()+${v}.z();\n complex<${var}> ${v}TMZ = ${v}.t()-${v}.z();\n complex<${var}> ${v}XPY = ${v}.x()+Complex(0.,1.)*${v}.y();\n complex<${var}> ${v}XMY = ${v}.x()-Complex(0.,1.)*${v}.y();") vslashM = Template("Matrix([[$m,0,${v}TMZ,-${v}XMY],[0,$m,-${v}XPY,${v}TPZ],[${v}TPZ,${v}XMY,$m,0],[${v}XPY,${v}TMZ,0,$m]])") vslashM2 = Template("Matrix([[$m,0,-${v}TMZ,${v}XMY],[0,$m,${v}XPY,-${v}TPZ],[-${v}TPZ,-${v}XMY,$m,0],[-${v}XPY,-${v}TMZ,0,$m]])") vslashMS = Template("${v}TPZ=Symbol(\"${v}TPZ\")\n${v}TMZ=Symbol(\"${v}TMZ\")\n${v}XPY=Symbol(\"${v}XPY\")\n${v}XMY=Symbol(\"${v}XMY\")\n${m}=Symbol(\"${m}\")\nO${m}=Symbol(\"O${m}\")\n") rslash = Template("Matrix([[$m,0,${v}TMZ,-${v}XMY],[0,$m,-${v}XPY,${v}TPZ],[${v}TPZ,${v}XMY,$m,0],[${v}XPY,${v}TMZ,0,$m]])*( ($${eta}-2./3.*O${m}**2*${v}$${A}*${v}$${B})*Matrix([[1.,0,0,0],[0,1.,0,0],[0,0,1.,0],[0,0,0,1.]]) -1./3.*$${DA}*$${DB} -1./3.*O${m}*(${v}$${B}*$${DA}-${v}$${A}*$${DB}))") rslashB = Template("Matrix([[$m,0,${v}TMZ,-${v}XMY],[0,$m,-${v}XPY,${v}TPZ],[${v}TPZ,${v}XMY,$m,0],[${v}XPY,${v}TMZ,0,$m]])*( (${v2}$${A}-2./3.*O${m}**2*${v}$${A}*${dot})*Matrix([[1.,0,0,0],[0,1.,0,0],[0,0,1.,0],[0,0,0,1.]]) -1./3.*$${DA}*${DB} -1./3.*O${m}*(${dot}*$${DA}-${v}$${A}*${DB}))") rslash2 = Template("Matrix([[$m,0,-${v}TMZ,${v}XMY],[0,$m,${v}XPY,-${v}TPZ],[-${v}TPZ,-${v}XMY,$m,0],[-${v}XPY,-${v}TMZ,0,$m]])*( ($${eta}-2./3.*O${m}**2*${v}$${A}*${v}$${B})*Matrix([[1.,0,0,0],[0,1.,0,0],[0,0,1.,0],[0,0,0,1.]]) -1./3.*$${DA}*$${DB} +1./3.*O${m}*(${v}$${B}*$${DA}-${v}$${A}*$${DB}))") rslash2B = Template("Matrix([[$m,0,-${v}TMZ,${v}XMY],[0,$m,${v}XPY,-${v}TPZ],[-${v}TPZ,-${v}XMY,$m,0],[-${v}XPY,-${v}TMZ,0,$m]])*( (${v2}$${B}-2./3.*O${m}**2*${dot}*${v}$${B})*Matrix([[1.,0,0,0],[0,1.,0,0],[0,0,1.,0],[0,0,0,1.]]) -1./3.*${DA}*$${DB} +1./3.*O${m}*(${v}$${B}*${DA}-${dot}*$${DB}))") dirac=["Matrix([[0,0,1.,0],[0,0,0,1.],[1.,0,0,0],[0,1.,0,0]])","Matrix([[0,0,0,1.],[0,0,1.,0],[0,-1.,0,0],[-1.,0,0,0]])", "Matrix([[0,0,0,complex(0, -1.)],[0,0,complex(0, 1.),0],[0,complex(0, 1.),0,0],[complex(0, -1.),0,0,0]])", "Matrix([[0,0,1.,0],[0,0,0,-1.],[-1.,0,0,0],[0,1.,0,0]])"] CC = "Matrix([[0,1.,0,0],[-1.,0,0,0],[0,0,0,-1.],[0,0,1.,0]])" CD = "Matrix([[0,-1.,0,0],[1.,0,0,0],[0,0,0,1.],[0,0,-1.,0]])" evaluateTemplate = """\ {decl} {{ {momenta} {waves} {swap} {symbols} {couplings} {defns} {result} }} """ spinor = Template("Matrix([[${s}s1],[${s}s2],[${s}s3],[${s}s4]])") sbar = Template("Matrix([[${s}s1,${s}s2,${s}s3,${s}s4]])") sline = Template("${s}s1=Symbol(\"${s}s1\")\n${s}s2=Symbol(\"${s}s2\")\n${s}s3=Symbol(\"${s}s3\")\n${s}s4=Symbol(\"${s}s4\")\n") RSSpinorTemplate = Template("${type}(${outxs1},\n${outxs2},\n${outxs3},\n${outxs4},\n${outys1},\n${outys2},\n${outys3},\n${outys4},\n${outzs1},\n${outzs2},\n${outzs3},\n${outzs4},\n${outts1},\n${outts2},\n${outts3},\n${outts4})") SpinorTemplate = Template("${type}(${outs1},\n${outs2},\n${outs3},\n${outs4})") class LorentzIndex : """ A simple classs to store a Lorentz index """ type="" value=0 dimension=0 def __repr__(self): if(self.type=="V" and not isinstance(self.value,int)) : return self.value else : return "%s%s" % (self.type,self.value) def __init__(self,val) : if(isinstance(val,int)) : self.dimension=0 if(val<0) : self.type="D" self.value = val elif(val>0 and val/1000==0) : self.type="E" self.value = val elif(val>0 and val/1000==1) : self.type="T1" self.value = val%1000 elif(val>0 and val/1000==2) : self.type="T2" self.value = val%1000 else : - print "Unknown value in Lorentz index:",val + print("Unknown value in Lorentz index:",val) raise SkipThisVertex() else : - print "Unknown value in Lorentz index:",val + print("Unknown value in Lorentz index:",val) raise SkipThisVertex() def __eq__(self,other): if(not isinstance(other,LorentzIndex)) : return False return ( (self.type, self.value) == (other.type, other.value) ) def __hash__(self) : return hash((self.type,self.value)) class DiracMatrix: """A simple class to store Dirac matrices""" name ="" value="" index=0 def __init(self) : self.name="" self.value="" self.index=0 def __repr__(self) : if(self.value==0) : return "%s" % self.index else : return "%s" % self.value class LorentzStructure: """A simple class to store a Lorentz structures""" name="" value=0 lorentz=[] spin=[] def __init(self) : self.name="" self.value=0 self.lorentz=[] self.spin=[] def __repr__(self): output = self.name if((self.name=="P" or self.name=="Tensor") and self.value!=0) : output += "%s" % self.value if(self.name=="int" or self.name=="sign") : output += "=%s" % self.value elif(len(self.spin)==0) : output += "(" for val in self.lorentz : output += "%s," % val output=output.rstrip(",") output+=")" elif(len(self.lorentz)==0) : output += "(" for val in self.spin : output += "%s," % val output=output.rstrip(",") output+=")" else : output += "(" for val in self.lorentz : output += "%s," % val for val in self.spin : output += "%s," % val output=output.rstrip(",") output+=")" return output def LorentzCompare(a,b) : if(a.name=="int" and b.name=="int") : return int(abs(b.value)-abs(a.value)) elif(a.name=="int") : return -1 elif(b.name=="int") : return 1 elif(len(a.spin)==0) : if(len(b.spin)==0) : return len(b.lorentz)-len(a.lorentz) else : return -1 elif(len(b.spin)==0) : return 1 else : if(len(a.spin)==0 or len(b.spin)==0) : - print 'Index problem in lorentz compare',\ - a.name,b.name,a.spin,b.spin + print('Index problem in lorentz compare', + a.name,b.name,a.spin,b.spin) raise SkipThisVertex() if(a.spin[0]>0 or b.spin[1]>0 ) : return -1 if(a.spin[1]>0 or b.spin[0]>0 ) : return 1 if(a.spin[1]==b.spin[0]) : return -1 if(b.spin[1]==a.spin[0]) : return 1 return 0 def extractIndices(struct) : if(struct.find("(")<0) : return [] temp=struct.split("(")[1].split(")")[0].split(",") output=[] for val in temp : output.append(int(val)) return output def parse_structure(structure,spins) : output=[] found = True while(found) : found = False # signs between terms if(structure=="+" or structure=="-") : output.append(LorentzStructure()) output[0].name="sign" output[0].value=structure[0]+"1." output[0].value=float(output[0].value) output[0].lorentz=[] output[0].spin=[] return output # simple numeric pre/post factors elif((structure[0]=="-" or structure[0]=="+") and structure[-1]==")" and structure[1]=="(") : output.append(LorentzStructure()) output[-1].name="int" output[-1].value=structure[0]+"1." output[-1].value=float(output[-1].value) output[-1].lorentz=[] output[-1].spin=[] structure=structure[2:-1] found=True elif(structure[0]=="(") : temp=structure.rsplit(")",1) structure=temp[0][1:] output.append(LorentzStructure()) output[-1].name="int" output[-1].value="1."+temp[1] output[-1].value=float(eval(output[-1].value)) output[-1].lorentz=[] output[-1].spin=[] found=True elif(structure[0:2]=="-(") : temp=structure.rsplit(")",1) structure=temp[0][2:] output.append(LorentzStructure()) output[-1].name="int" output[-1].value="-1."+temp[1] output[-1].value=float(eval(output[-1].value)) output[-1].lorentz=[] output[-1].spin=[] found=True # special handling for powers power = False if("**" in structure ) : power = True structure = structure.replace("**","^") structures = structure.split("*") if(power) : for j in range(0,len(structures)): if(structures[j].find("^")>=0) : temp = structures[j].split("^") structures[j] = temp[0] for i in range(0,int(temp[1])-1) : structures.append(temp[0]) # split up the structure for struct in structures: ind = extractIndices(struct) # different types of object # object with only spin indices if(struct.find("Identity")==0 or struct.find("Proj")==0 or struct.find("Gamma5")==0) : output.append(LorentzStructure()) output[-1].spin=ind output[-1].lorentz=[] output[-1].name=struct.split("(")[0] output[-1].value=0 if(len(struct.replace("%s(%s,%s)" % (output[-1].name,ind[0],ind[1]),""))!=0) : - print "Problem handling %s structure " % output[-1].name + print("Problem handling %s structure " % output[-1].name) raise SkipThisVertex() # objects with 2 lorentz indices elif(struct.find("Metric")==0) : output.append(LorentzStructure()) output[-1].lorentz=[LorentzIndex(ind[0]),LorentzIndex(ind[1])] output[-1].name=struct.split("(")[0] output[-1].value=0 output[-1].spin=[] if(len(struct.replace("%s(%s,%s)" % (output[-1].name,ind[0],ind[1]),""))!=0) : - print "Problem handling %s structure " % output[-1].name + print("Problem handling %s structure " % output[-1].name) raise SkipThisVertex() elif(struct.find("P(")==0) : output.append(LorentzStructure()) output[-1].lorentz=[LorentzIndex(ind[0])] output[-1].name=struct.split("(")[0] output[-1].value=ind[1] output[-1].spin=[] if(len(struct.replace("%s(%s,%s)" % (output[-1].name,ind[0],ind[1]),""))!=0) : - print "Problem handling %s structure " % output[-1].name + print("Problem handling %s structure " % output[-1].name) raise SkipThisVertex() # 1 lorentz and 1 spin index elif(struct.find("Gamma")==0) : output.append(LorentzStructure()) output[-1].lorentz=[LorentzIndex(ind[0])] output[-1].spin=[ind[1],ind[2]] output[-1].name=struct.split("(")[0] output[-1].value=1 if(len(struct.replace("%s(%s,%s,%s)" % (output[-1].name,ind[0],ind[1],ind[2]),""))!=0) : - print "problem parsing gamma matrix",struct + print("problem parsing gamma matrix",struct) raise SkipThisVertex() # objects with 4 lorentz indices elif(struct.find("Epsilon")==0) : output.append(LorentzStructure()) output[-1].lorentz=[] for i in range(0,len(ind)) : output[-1].lorentz.append(LorentzIndex(ind[i])) output[-1].spin=[] output[-1].name=struct.split("(")[0] output[-1].value=1 if(len(struct.replace("%s(%s,%s,%s,%s)" % (output[-1].name,ind[0],ind[1],ind[2],ind[3]),""))!=0) : - print 'Problem parsing epsilon',struct + print('Problem parsing epsilon',struct) raise SkipThisVertex() # scalars else : try : output.append(LorentzStructure()) output[-1].value=float(struct) output[-1].name="int" output[-1].lorentz=[] output[-1].spin=[] except : if(struct.find("complex")==0) : vals = struct[0:-1].replace("complex(","").split(",") output[-1].value=complex(float(vals[0]),float(vals[1])) output[-1].name="int" output[-1].lorentz=[] output[-1].spin=[] else : - print 'Problem parsing scalar',struct + print('Problem parsing scalar',struct) raise SkipThisVertex() # now do the sorting if(len(output)==1) : return output output = sorted(output,cmp=LorentzCompare) # fix indices in the RS case if(4 in spins) : for i in range(0,len(output)) : for ll in range(0,len(output[i].lorentz)) : if(spins[output[i].lorentz[ll].value-1]==4 and output[i].lorentz[ll].type=="E") : output[i].lorentz[ll].type="R" # return the answer return output def constructDotProduct(ind1,ind2,defns) : (ind1,ind2) = sorted((ind1,ind2),cmp=indSort) dimension=ind1.dimension+ind2.dimension # this product already dealt with ? if((ind1,ind2) in defns) : name = defns[(ind1,ind2)][0] # handle the product else : name = "dot%s" % (len(defns)+1) unit = computeUnit(dimension) defns[(ind1,ind2)] = [name,"complex<%s> %s = %s*%s;" % (unit,name,ind1,ind2)] return (name,dimension) def contract(parsed) : for j in range(0,len(parsed)) : if(parsed[j]=="") : continue if(parsed[j].name=="P") : # simplest case if(parsed[j].lorentz[0].type=="E" or parsed[j].lorentz[0].type=="P") : newIndex = LorentzIndex(parsed[j].value) newIndex.type="P" newIndex.dimension=1 parsed[j].name="Metric" parsed[j].lorentz.append(newIndex) parsed[j].lorentz = sorted(parsed[j].lorentz,cmp=indSort) continue ll=1 found=False for k in range(0,len(parsed)) : if(j==k or parsed[k]=="" ) : continue for i in range(0,len(parsed[k].lorentz)) : if(parsed[k].lorentz[i] == parsed[j].lorentz[0]) : parsed[k].lorentz[i].type="P" parsed[k].lorentz[i].value = parsed[j].value parsed[k].lorentz[i].dimension=1 if(parsed[k].name=="P") : parsed[k].lorentz.append(LorentzIndex(parsed[k].value)) parsed[k].lorentz[1].type="P" parsed[k].lorentz[1].dimension=1 parsed[k].name="Metric" parsed[k].value = 0 found=True break if(found) : parsed[j]="" break return [x for x in parsed if x != ""] def computeUnit(dimension) : if(isinstance(dimension,int)) : dtemp = dimension else : dtemp=dimension[1]+dimension[2] if(dtemp==0) : unit="double" elif(dtemp==1) : unit="Energy" elif(dtemp==-1) : unit="InvEnergy" elif(dtemp>0) : unit="Energy%s" % (dtemp) elif(dtemp<0) : unit="InvEnergy%s" % (dtemp) return unit def computeUnit2(dimension,vDim) : # first correct for any coupling power in vertex totalDim = int(dimension[0])+dimension[2]+vDim-4 output="" if(totalDim!=0) : if(totalDim>0) : if(totalDim==1) : output = "1./GeV" elif(totalDim==2) : output = "1./GeV2" else : output="1." for i in range(0,totalDim) : output +="/GeV" else : if(totalDim==-1) : output = "GeV" elif(totalDim==-2) : output = "GeV2" else : output="1." for i in range(0,-totalDim) : output +="*GeV" expr="" # now remove the remaining dimensionality removal=dimension[1]-int(dimension[0])-vDim+4 if(removal!=0) : if(removal>0) : if(removal==1) : expr = "UnitRemovalInvE" else : expr = "UnitRemovalInvE%s" % removal else : if(removal==-1) : expr = "UnitRemovalE" else : expr = "UnitRemovalE%s" % (-removal) if(output=="") : return expr elif(expr=="") : return output else : return "%s*%s" %(output,expr) # order the indices of a dot product def indSort(a,b) : if(not isinstance(a,LorentzIndex) or not isinstance(b,LorentzIndex)) : - print "Trying to sort something that's not a Lorentz index",a,b + print("Trying to sort something that's not a Lorentz index",a,b) raise SkipThisVertex() if(a.type==b.type) : i1=a.value i2=b.value if(i1>i2) : return 1 elif(i10) : output="+" else : output="-" parsed=[] return (output,parsed,dimension,eps) # replace integers (really lorentz scalars) for j in range(0,len(parsed)) : if(parsed[j]!="" and parsed[j].name=="int") : output *= parsed[j].value parsed[j]="" # bracket this for safety if(output!="") : output = "(%s)" % output # special for tensor indices if("T" in lorentztag) : for j in range(0,len(parsed)) : if(parsed[j]=="") :continue # check for tensor index found=False for li in parsed[j].lorentz : if(li.type[0]=="T") : index = li found=True break if(not found) : continue # workout the other index for the tensor index2 = LorentzIndex(li.value) if(index.type=="T1") : index2.type="T2" else : index2.type="T1" # special is tensor contracted with itself if(parsed[j].name=="Metric" and index2 == parsed[j].lorentz[1]) : parsed[j].name = "Tensor" parsed[j].value = index.value parsed[j].lorentz = [] if(iloc!=index.value) : name= "traceT%s" % parsed[j].value if( name in defns ) : output += "*(%s)" % defns[name][0] else : defns[name] = [name,"Complex %s = T%s.trace();" % (name,parsed[j].value)] output += "*(%s)" % defns[name][0] parsed[j]="" continue # otherwise search for the match for k in range(j+1,len(parsed)) : found = False for li in parsed[k].lorentz : if(li == index2) : found=True break if(not found) : continue if(parsed[j].name=="P") : newIndex1 = LorentzIndex(parsed[j].value) newIndex1.type="P" newIndex1.dimension=1 elif(parsed[j].name=="Metric") : for li in parsed[j].lorentz : if(li != index) : newIndex1=li break else : - print 'Unknown object with tensor index, first object',parsed[j] + print('Unknown object with tensor index, first object',parsed[j]) raise SkipThisVertex() if(parsed[k].name=="P") : newIndex2 = LorentzIndex(parsed[k].value) newIndex2.type="P" newIndex2.dimension=1 elif(parsed[k].name=="Metric") : for li in parsed[k].lorentz : if(li != index) : newIndex2=li break elif(parsed[k].name=="Gamma") : # if we can't contract if(index.value==iloc or (newIndex1.type=="E" and newIndex1.value==iloc)) : newIndex2=index2 # otherwise contract else : unit=computeUnit(newIndex1.dimension) if(index.type=="T1") : name="T%s%sF" % (index.value,newIndex1) defns[name] = [name,"LorentzVector > %s = T%s.preDot(%s);" % (unit,name,index.value,newIndex1)] else : name="T%s%sS" % (index.value,newIndex1) defns[name] = [name,"LorentzVector > %s = T%s.postDot(%s);" % (unit,name,index.value,newIndex1)] parsed[j]="" gIndex=LorentzIndex(-1) gIndex.type="V" gIndex.value=name gIndex.dimension=newIndex1.dimension parsed[k].lorentz[0] = gIndex break else : - print 'Unknown object with tensor index, second object',parsed[j],parsed[k] + print('Unknown object with tensor index, second object',parsed[j],parsed[k]) raise SkipThisVertex() if(index2.type=="T1") : newIndex1,newIndex2=newIndex2,newIndex1 parsed[j].name = "Tensor" parsed[j].value= int(index.value) parsed[j].lorentz= [newIndex1,newIndex2] if(parsed[k].name!="Gamma") : parsed[k]="" break # main handling of lorentz structures for j in range(0,len(parsed)) : if(parsed[j]=="") : continue if(parsed[j].name=="Metric") : # check whether or not we can contract canContract=True for ll in parsed[j].lorentz : if((ll.type=="E" and ll.value==iloc) or ll.type=="R") : canContract = False break if(not canContract) : continue # if we can do it (name,dTemp) = constructDotProduct(parsed[j].lorentz[0],parsed[j].lorentz[1],defns) output += "*(%s)" % name dimension[2] += dTemp parsed[j]="" elif(parsed[j].name=="Epsilon") : if(not eps) : eps = True # work out which, if any of the indices can be summed over summable=[] for ix in range(0,len(parsed[j].lorentz)) : if(parsed[j].lorentz[ix].type=="P" or (parsed[j].lorentz[ix].type=="E" and iloc !=parsed[j].lorentz[ix].value)) : summable.append(True) else : summable.append(False) sc = summable.count(True) # less than 3 contractable indices, leave for later if(sc<3) : continue # can contract to a vector elif(sc==3) : offLoc = -1 for i in range(0,len(summable)): if(not summable[i]) : offLoc = i break else : offLoc = 0 indices=[] dTemp=0 for ix in range(0,len(parsed[j].lorentz)) : dTemp += parsed[j].lorentz[ix].dimension if(ix!=offLoc) : indices.append(parsed[j].lorentz[ix]) # contract all the indices if(sc==4) : dimension[2] += dTemp iTemp = (parsed[j].lorentz[0],parsed[j].lorentz[1], parsed[j].lorentz[2],parsed[j].lorentz[3]) if(iTemp in defns) : output += "*(%s)" % defns[iTemp][0] parsed[j]="" else : name = "dot%s" % (len(defns)+1) unit = computeUnit(dTemp) defns[iTemp] = [name,"complex<%s> %s =-%s*epsilon(%s,%s,%s);" % (unit,name,parsed[j].lorentz[0], indices[0],indices[1],indices[2]) ] output += "*(%s)" % name parsed[j]="" # contract 3 indices leaving a vector else : iTemp = (indices[0],indices[1],indices[2]) sign = "1" if(offLoc%2!=0) : sign="-1" if(iTemp in defns) : name = defns[iTemp][0] else : name = "V%s" % (len(defns)+1) unit = computeUnit(dTemp) defns[iTemp] = [name,"LorentzVector > %s =-epsilon(%s,%s,%s);" % (unit,name, indices[0],indices[1],indices[2]) ] newIndex = LorentzIndex(int(name[1:])) newIndex.type="V" newIndex.dimension=dTemp output += "*(%s)" % (sign) oi = parsed[j].lorentz[offLoc] if(oi.type!="D") : parsed[j].name="Metric" parsed[j].spins=[] parsed[j].value=0 parsed[j].lorentz=[newIndex,oi] else : found=False for k in range(0,len(parsed)): if(k==j or parsed[k]=="") : continue for ll in range(0,len(parsed[k].lorentz)) : if(parsed[k].lorentz[ll]==oi) : found=True parsed[k].lorentz[ll]=newIndex break if(found) : break if(not found) : - print "Problem contracting indices of Epsilon tensor" + print("Problem contracting indices of Epsilon tensor") raise SkipThisVertex() parsed[j]="" elif(parsed[j].name=="Tensor") : # not an external tensor if(parsed[j].value!=iloc) : # now check the lorentz indices con=[] uncon=[] dtemp=0 for li in parsed[j].lorentz : if(li.type=="P" or li.type=="V") : con.append(li) dtemp+=li.dimension elif(li.type=="E") : if(li.value!=iloc) : con.append(li) else : uncon.append(li) else : - print "Can't handle index ",li,"in tensor",parsed[j] + print("Can't handle index ",li,"in tensor",parsed[j]) raise SkipThisVertex() if(len(con)==2) : iTemp = ("T%s%s%s"% (parsed[j].value,con[0],con[1])) dimension[2]+=dtemp if(iTemp in defns) : output += "*(%s)" % defns[iTemp][0] else : unit=computeUnit(dtemp) name = "dot%s" % (len(defns)+1) defns[iTemp] = [name,"complex<%s> %s = T%s.preDot(%s)*%s;" % (unit,name,parsed[j].value,con[0],con[1])] output += "*(%s)" % name parsed[j]="" # handled in final stage else : continue elif(parsed[j].name.find("Proj")>=0 or parsed[j].name.find("Gamma")>=0 or parsed[j].name.find("Identity")>=0) : continue elif(parsed[j].name=="P" and parsed[j].lorentz[0].type=="R") : continue else : - print 'Lorentz structure',parsed[j],'not handled' + print('Lorentz structure',parsed[j],'not handled') raise SkipThisVertex() # remove leading * if(output!="" and output[0]=="*") : output = output[1:] # remove any (now) empty elements parsed = [x for x in parsed if x != ""] return (output,parsed,dimension,eps) def finalContractions(output,parsed,dimension,lorentztag,iloc,defns) : if(len(parsed)==0) : return (output,dimension) elif(len(parsed)!=1) : - print "Summation can't be handled",parsed + print("Summation can't be handled",parsed) raise SkipThisVertex() if(parsed[0].name=="Tensor") : # contracted with off-shell vector if(parsed[0].value!=iloc) : found = False for ll in parsed[0].lorentz : if(ll.type=="E" and ll.value==iloc) : found = True else : lo=ll if(found) : dimension[2]+= lo.dimension unit=computeUnit(lo.dimension) if(lo==parsed[0].lorentz[0]) : name="T%s%sF" % (parsed[0].value,lo) defns[name] = [name,"LorentzVector > %s = T%s.preDot(%s);" % (unit,name,parsed[0].value,lo)] else : name="T%s%sS" % (parsed[0].value,lo) defns[name] = [name,"LorentzVector > %s = T%s.postDot(%s);" % (unit,name,parsed[0].value,lo)] parsed[0]="" if(output=="") : output="1." output = "(%s)*(%s)" %(output,name) else : - print "Can\'t contract tensor",lo,iloc + print("Can\'t contract tensor",lo,iloc) raise SkipThisVertex() # off-shell tensor else : if(len(parsed[0].lorentz)!=0) : dimension[2]+=parsed[0].lorentz[0].dimension+parsed[0].lorentz[1].dimension tensor = tensorPropagator(parsed[0],defns) if(output=="") : output="1." output = [output,tensor,()] elif(parsed[0].name=="Metric") : found = False for ll in parsed[0].lorentz : if(ll.type=="E" and ll.value==iloc) : found = True else : lo=ll if(found) : parsed[0]="" dimension[2] += lo.dimension if(output=="") : output="1." output = "(%s)*(%s)" %(output,lo) else : - print "Structure can't be handled",parsed,iloc + print("Structure can't be handled",parsed,iloc) raise SkipThisVertex() return (output,dimension) def tensorPropagator(struct,defns) : # dummy index i0 = LorentzIndex(-1000) # index for momentum of propagator ip = LorentzIndex(struct.value) ip.type="P" ip.dimension=1 # the metric tensor terms=[] if(len(struct.lorentz)==0) : (dp,dTemp) = constructDotProduct(ip,ip,defns) pre = "-1./3.*(1.-%s*OM%s)" % (dp,struct.value) terms.append((pre,i0,i0)) pre = "-2./3.*(1.-%s*OM%s)" % (dp,struct.value) terms.append(("%s*OM%s" %(pre,struct.value),ip,ip)) else : # indices of the tensor ind1 = struct.lorentz[0] ind2 = struct.lorentz[1] # the dot products we need (d1,dtemp) = constructDotProduct(ind1,ip,defns) (d2,dtemp) = constructDotProduct(ind2,ip,defns) (d3,dtemp) = constructDotProduct(ind1,ind2,defns) # various terms in the propagator terms.append(("0.5",ind1,ind2)) terms.append(("-0.5*OM%s*%s"%(struct.value,d1),ip,ind2)) terms.append(("-0.5*OM%s*%s"%(struct.value,d2),ind1,ip)) terms.append(("0.5",ind2,ind1)) terms.append(("-0.5*OM%s*%s"%(struct.value,d2),ip,ind1)) terms.append(("-0.5*OM%s*%s"%(struct.value,d1),ind2,ip)) terms.append(("-1./3.*"+d3,i0,i0)) terms.append(("1./3.*OM%s*%s*%s"%(struct.value,d1,d2),i0,i0)) terms.append(("1./3.*OM%s*%s"%(struct.value,d3),ip,ip)) terms.append(("2./3.*OM%s*OM%s*%s*%s"%(struct.value,struct.value,d1,d2),ip,ip)) # compute the output as a dict output={} for i1 in imap: for i2 in imap : val="" for term in terms: if(term[0][0]!="-") : pre = "+"+term[0] else : pre = term[0] if(term[1]==i0) : if(i1==i2) : if(i1=="t") : val += pre else : if(pre[0]=="+") : val +="-"+pre[1:] else : val +="+"+pre[1:] else : val += "%s*%s%s*%s%s" % (pre,term[1],i1,term[2],i2) output["%s%s" % (i1,i2) ] = val.replace("+1*","+").replace("-1*","-") return output def generateVertex(iloc,L,parsed,lorentztag,vertex,defns) : # try to import sympy and exit if required try : import sympy from sympy import Matrix,Symbol except : - print 'ufo2herwig uses the python sympy module to translate general lorentz structures.' - print 'This must be installed if you wish to use this option.' - print 'EXITTING' + print('ufo2herwig uses the python sympy module to translate general lorentz structures.') + print('This must be installed if you wish to use this option.') + print('EXITTING') quit() eps=False # parse the lorentz structures output = [1.]*len(parsed) dimension=[] for i in range(0,len(parsed)) : dimension.append([0,0,0]) for i in range (0,len(parsed)) : (output[i],parsed[i],dimension[i],eps) = finishParsing(parsed[i],dimension[i],lorentztag,iloc,defns,eps) # still need to process gamma matrix strings for fermions if(lorentztag[0] in ["F","R"] ) : return convertDirac(output,dimension,eps,iloc,L,parsed,lorentztag,vertex,defns) # return the answer else : handled=True for i in range (0,len(parsed)) : if(len(parsed[i])!=0) : handled = False break if(not handled) : for i in range (0,len(parsed)) : (output[i],dimension[i]) = finalContractions(output[i],parsed[i],dimension[i],lorentztag,iloc,defns) return (output,dimension,eps) def convertDirac(output,dimension,eps,iloc,L,parsed,lorentztag,vertex,defns) : for i in range(0,len(parsed)): # skip empty elements if(len(parsed[i])==0 or (len(parsed[i])==1 and parsed[i][0]=="")) : continue # parse the string (output[i],dimension[i],defns) = convertDiracStructure(parsed[i],output[i],dimension[i], defns,iloc,L,lorentztag,vertex) return (output,dimension,eps) # parse the gamma matrices def convertMatrix(structure,spins,unContracted,Symbols,dtemp,defns,iloc) : i1 = structure.spin[0] i2 = structure.spin[1] if(structure.name=="Identity") : output = DiracMatrix() output.value = I4 output.index=0 output.name="M" structure="" elif(structure.name=="Gamma5") : output = DiracMatrix() output.value = G5 output.index=0 output.name="M" structure="" elif(structure.name=="ProjM") : output = DiracMatrix() output.value = PM output.index=0 output.name="M" structure="" elif(structure.name=="ProjP") : output = DiracMatrix() output.value = PP output.index=0 output.name="M" structure="" elif(structure.name=="Gamma") : # gamma(mu) lorentz matrix contracted with dummy index if(structure.lorentz[0].type=="D" or structure.lorentz[0].type=="R") : if(structure.lorentz[0] not in unContracted) : unContracted[structure.lorentz[0]] = 0 output = DiracMatrix() output.value=0 output.index=structure.lorentz[0] output.name="GMU" structure="" elif(structure.lorentz[0].type == "E" and structure.lorentz[0].value == iloc ) : if(structure.lorentz[0] not in unContracted) : unContracted[structure.lorentz[0]] = 0 output = DiracMatrix() output.value=0 output.index=structure.lorentz[0] output.name="GMU" structure="" elif(structure.lorentz[0].type == "T1" or structure.lorentz[0].type == "T2") : if(structure.lorentz[0] not in unContracted) : unContracted[structure.lorentz[0]] = 0 output = DiracMatrix() output.value=0 output.index=structure.lorentz[0] output.name="GMU" structure="" else : output=DiracMatrix() output.name="M" output.value = vslash.substitute({ "v" : structure.lorentz[0]}) Symbols += vslashS.substitute({ "v" : structure.lorentz[0]}) variable = computeUnit(structure.lorentz[0].dimension) #if(structure.lorentz[0].type!="V" or # structure.lorentz[0].type=="V") : dtemp[2] += structure.lorentz[0].dimension defns["vv%s" % structure.lorentz[0] ] = \ ["vv%s" % structure.lorentz[0], vslashD.substitute({ "var" : variable, "v" : structure.lorentz[0]})] structure="" else : - print 'Unknown Gamma matrix structure',structure + print('Unknown Gamma matrix structure',structure) raise SkipThisVertex() return (i1,i2,output,structure,Symbols) def checkRSContract(parsed,loc,dtemp) : rindex=LorentzIndex(loc) rindex.type="R" contract="" for i in range(0,len(parsed)) : if(parsed[i]=="") : continue found = False for ll in range(0,len(parsed[i].lorentz)) : if(parsed[i].lorentz[ll]==rindex) : found=True break if(not found) : continue if(parsed[i].name=="P") : dtemp[2]+=1 contract = LorentzIndex(parsed[i].value) contract.type="P" contract.dimension=1 parsed[i]="" break elif(parsed[i].name=="Metric") : for ll in parsed[i].lorentz : if(ll==rindex) : continue else : break if(ll.type=="P") : dtemp[2]+=1 contract=ll parsed[i]="" break elif(parsed[i].name=="Epsilon") : continue else : - print "Unkonwn type contracted with RS spinor",parsed[i] + print("Unkonwn type contracted with RS spinor",parsed[i]) raise SkipThisVertex() return contract def processChain(dtemp,parsed,spins,Symbols,unContracted,defns,iloc) : # piece of dimension which is common (0.5 for sbar and spinor) dtemp[0]+=1 # set up the spin indices sind = 0 lind = 0 expr=[] # now find the next thing in the matrix string ii = 0 index=0 while True : # already handled if(parsed[ii]=="") : ii+=1 continue # start of the chain elif(sind==0 and len(parsed[ii].spin)==2 and parsed[ii].spin[0]>0 ) : (sind,index,matrix,parsed[ii],Symbols) \ = convertMatrix(parsed[ii],spins,unContracted,Symbols,dtemp,defns,iloc) expr.append(matrix) # next element in the chain elif(index!=0 and len(parsed[ii].spin)==2 and parsed[ii].spin[0]==index) : (crap,index,matrix,parsed[ii],Symbols) \ = convertMatrix(parsed[ii],spins,unContracted,Symbols,dtemp,defns,iloc) expr.append(matrix) # check the index to see if we're at the end if(index>0) : lind=index break ii+=1 if(ii>=len(parsed)) : - print "Can't parsed the chain of dirac matrices" - print parsed + print("Can't parsed the chain of dirac matrices") + print(parsed) raise SkipThisVertex() # start and end of the spin chains # first particle spin 1/2 if(spins[sind-1]==2) : start = DiracMatrix() endT = DiracMatrix() start.index=0 endT .index=0 # start of chain and end of transpose # off shell if(sind==iloc) : start.name="M" endT .name="M" start.value = vslashM .substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind} ) Symbols += vslashMS.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind} ) endT.value = vslashM2.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind} ) defns["vvP%s" % sind ] = ["vvP%s" % sind , vslashD.substitute({ "var" : "Energy", "v" : "P%s" % sind })] dtemp[1]+=1 # onshell else : start.name="S" endT .name="S" subs = {'s' : ("sbar%s" % sind)} start.value = sbar .substitute(subs) Symbols += sline.substitute(subs) subs = {'s' : ("s%s" % sind)} endT.value = spinor.substitute(subs) Symbols += sline.substitute(subs) # spin 3/2 fermion elif spins[sind-1]==4 : # check if we can easily contract contract=checkRSContract(parsed,sind,dtemp) # off-shell if(sind==iloc) : oindex = LorentzIndex(sind) oindex.type="O" unContracted[oindex]=0 Symbols += vslashMS.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind} ) Symbols += momCom.substitute({"v" : "P%s" %sind }) defns["vvP%s" % sind ] = ["vvP%s" % sind , vslashD.substitute({ "var" : "Energy", "v" : "P%s" % sind })] dtemp[1] += 1 if(contract=="") : rindex = LorentzIndex(sind) rindex.type="R" start=DiracMatrix() start.value = Template(rslash.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind, "loc" : sind })) start.name = "RP" start.index = (oindex,rindex) endT=DiracMatrix() endT.value = Template(rslash2.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind, "loc" : sind })) endT.name = "RP" endT.index = (rindex,oindex) else : # construct dot product pi = LorentzIndex(sind) pi.type="P" pi.dimension=1 (name,dummy) = constructDotProduct(pi,contract,defns) Symbols += momCom.substitute({"v" : contract }) RB = vslash.substitute({ "v" : contract}) Symbols += vslashS.substitute({ "v" : contract }) Symbols += "%s = Symbol('%s')\n" % (name,name) defns["vv%s" % contract ] = ["vv%s" % contract, vslashD.substitute({ "var" : computeUnit(contract.dimension), "v" : "%s" % contract })] start=DiracMatrix() start.name="RQ" start.index = oindex start.value =Template( rslashB.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind, "loc" : sind, "DB" : RB, "dot" : name , "v2" : contract})) endT=DiracMatrix() endT.name = "RQ" endT.index = oindex endT.value = Template(rslash2B.substitute({ "v" : "P%s" % sind, "m" : "M%s" % sind, "loc" : sind , "DA" : RB, "dot" : name , "v2" : contract})) # on-shell else : start = DiracMatrix() endT = DiracMatrix() # no contraction if(contract=="" or (contract.type=="E" and contract.value==iloc) ) : if contract == "" : contract = LorentzIndex(sind) contract.type="R" # start of matrix string start.value = Template(sbar .substitute({'s' : ("Rsbar%s${L}" % sind)})) start.name="RS" start.index = contract # end of transpose string endT.value=Template(spinor.substitute({'s' : ("Rs%s${L}" % sind)})) endT.name="RS" endT.index = contract unContracted[contract]=0 # variables for sympy for LI in imap : Symbols += sline.substitute({'s' : ("Rs%s%s" % (sind,LI))}) Symbols += sline.substitute({'s' : ("Rsbar%s%s" % (sind,LI))}) else : # start of matrix string start.name="S" start.value = "Matrix([[%s,%s,%s,%s]])" % (RSDotProduct.substitute({'s' : ("Rsbar%s" % sind), 'v':contract, 'si' : 1}), RSDotProduct.substitute({'s' : ("Rsbar%s" % sind), 'v':contract, 'si' : 2}), RSDotProduct.substitute({'s' : ("Rsbar%s" % sind), 'v':contract, 'si' : 3}), RSDotProduct.substitute({'s' : ("Rsbar%s" % sind), 'v':contract, 'si' : 4})) endT.name="S" endT.value = "Matrix([[%s],[%s],[%s],[%s]])" % (RSDotProduct.substitute({'s' : ("Rs%s" % sind), 'v':contract, 'si' : 1}), RSDotProduct.substitute({'s' : ("Rs%s" % sind), 'v':contract, 'si' : 2}), RSDotProduct.substitute({'s' : ("Rs%s" % sind), 'v':contract, 'si' : 3}), RSDotProduct.substitute({'s' : ("Rs%s" % sind), 'v':contract, 'si' : 4})) Symbols += momCom.substitute({"v" : contract }) for LI in ["x","y","z","t"] : Symbols += sline.substitute({'s' : ("Rs%s%s" % (sind,LI))}) Symbols += sline.substitute({'s' : ("Rsbar%s%s" % (sind,LI))}) # last particle spin 1/2 if( spins[lind-1]==2 ) : end = DiracMatrix() startT = DiracMatrix() end .index=0 startT.index=0 # end of chain if(lind==iloc) : end.name ="M" startT.name="M" end.value = vslashM2.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind} ) startT.value = vslashM.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind} ) Symbols += vslashMS.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind} ) defns["vvP%s" % lind ] = ["vvP%s" % lind , vslashD.substitute({ "var" : "Energy", "v" : "P%s" % lind })] dtemp[1] += 1 else : startT.name="S" end .name="S" subs = {'s' : ("s%s" % lind)} end.value = spinor.substitute(subs) Symbols += sline.substitute(subs) subs = {'s' : ("sbar%s" % lind)} startT.value = sbar .substitute(subs) Symbols += sline.substitute(subs) # last particle spin 3/2 elif spins[lind-1]==4 : # check if we can easily contract contract=checkRSContract(parsed,lind,dtemp) # off-shell if(lind==iloc) : oindex = LorentzIndex(lind) oindex.type="O" unContracted[oindex]=0 Symbols += vslashMS.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind} ) Symbols += momCom.substitute({"v" : "P%s" %lind }) defns["vvP%s" % lind ] = ["vvP%s" % lind , vslashD.substitute({ "var" : "Energy", "v" : "P%s" % lind })] dtemp[1] += 1 if(contract=="") : rindex = LorentzIndex(lind) rindex.type="R" end=DiracMatrix() end.value = Template(rslash2.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind, "loc" : lind })) end.name = "RP" end.index = (rindex,oindex) startT=DiracMatrix() startT.value=Template(rslash.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind, "loc" : lind })) startT.name = "RP" startT.index = (oindex,rindex) else : # construct dot product pi = LorentzIndex(lind) pi.type="P" pi.dimension=1 (name,unit) = constructDotProduct(pi,contract,defns) Symbols += momCom.substitute({"v" : contract }) RB = vslash.substitute({ "v" : contract}) Symbols += vslashS.substitute({ "v" : contract }) Symbols += "%s = Symbol('%s')\n" % (name,name) defns["vv%s" % contract ] = ["vv%s" % contract, vslashD.substitute({ "var" : computeUnit(contract.dimension), "v" : "%s" % contract })] end=DiracMatrix() end.name="RQ" end.index = oindex end.value =Template( rslash2B.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind, "loc" : lind, "DA" : RB, "dot" : name , "v2" : contract})) startT=DiracMatrix() startT.name = "RQ" startT.index = oindex startT.value = Template(rslashB.substitute({ "v" : "P%s" % lind, "m" : "M%s" % lind, "loc" : lind , "DB" : RB, "dot" : name , "v2" : contract})) # on-shell else : end = DiracMatrix() startT = DiracMatrix() # no contraction if(contract=="" or (contract.type=="E" and contract.value==iloc) ) : if contract == "" : contract = LorentzIndex(lind) contract.type="R" # end of matrix string end.value = Template(spinor.substitute({'s' : ("Rs%s${L}" % lind)})) end.name = "RS" end.index = contract # start of matrix string startT.value = Template(sbar .substitute({'s' : ("Rsbar%s${L}" % lind)})) startT.name = "RS" startT.index = contract unContracted[contract]=0 # variables for sympy for LI in imap : Symbols += sline.substitute({'s' : ("Rs%s%s" % (lind,LI))}) Symbols += sline.substitute({'s' : ("Rsbar%s%s" % (lind,LI))}) # contraction else : # end of the matrix string end.name = "S" end.value = "Matrix([[%s],[%s],[%s],[%s]])" % (RSDotProduct.substitute({'s' : ("Rs%s" % lind), 'v':contract, 'si' : 1}), RSDotProduct.substitute({'s' : ("Rs%s" % lind), 'v':contract, 'si' : 2}), RSDotProduct.substitute({'s' : ("Rs%s" % lind), 'v':contract, 'si' : 3}), RSDotProduct.substitute({'s' : ("Rs%s" % lind), 'v':contract, 'si' : 4})) startT.name = "S" startT.value = "Matrix([[%s,%s,%s,%s]])" % (RSDotProduct.substitute({'s' : ("Rsbar%s" % lind), 'v':contract, 'si' : 1}), RSDotProduct.substitute({'s' : ("Rsbar%s" % lind), 'v':contract, 'si' : 2}), RSDotProduct.substitute({'s' : ("Rsbar%s" % lind), 'v':contract, 'si' : 3}), RSDotProduct.substitute({'s' : ("Rsbar%s" % lind), 'v':contract, 'si' : 4})) Symbols += momCom.substitute({"v" : contract }) for LI in ["x","y","z","t"] : Symbols += sline.substitute({'s' : ("Rs%s%s" % (lind,LI))}) Symbols += sline.substitute({'s' : ("Rsbar%s%s" % (lind,LI))}) return(sind,lind,expr,start,startT,end,endT,Symbols) def calculateDirac(expr,start,end,startT,endT,sind,lind,Symbols,iloc) : res=[] for ichain in range(0,len(start)) : # calculate the matrix string etemp="*".join(str(x) for x in expr[ichain]) temp={} exec("import sympy\nfrom sympy import Symbol,Matrix\n"+Symbols+"result="+ ( "%s*%s*%s" %(start[ichain],etemp,end[ichain]))) in temp res.append(temp["result"]) tempT={} exec("import sympy\nfrom sympy import Symbol,Matrix,Transpose\n"+Symbols+"result="+ ( "%s*%s*Transpose(%s)*%s*%s" %(startT[ichain],CC,etemp,CD,endT[ichain]))) in tempT res.append(tempT["result"]) if(len(start)==1) : if(iloc==0 or (iloc!=sind[ichain] and iloc!=lind[ichain])) : sVal = {'s' : temp ["result"][0,0],'sT' : tempT["result"][0,0]} else : sVal={} for jj in range(1,5) : sVal["s%s" % jj] = temp ["result"][jj-1] sVal["sT%s" % jj] = tempT["result"][jj-1] else : sVal={} sVal["s" ] = res[0][0,0]*res[2][0,0] sVal["sT2" ] = res[0][0,0]*res[3][0,0] sVal["sT1" ] = res[1][0,0]*res[2][0,0] sVal["sT12"] = res[1][0,0]*res[3][0,0] return sVal def addToOutput(res,nchain,sign,rTemp) : # 1 spin chain if(nchain==1) : for ii in range(0,2) : if(rTemp[ii][0].shape[0]==1) : # result is scalar if(rTemp[ii][0].shape[1]==1) : if(len(res[ii])==0) : res[ii].append(sign*rTemp [ii][0][0,0]) else : res[ii][0] += sign*rTemp [ii][0][0,0] # result is a spinor elif(rTemp[ii][0].shape[1]==4) : if(len(res[ii])==0) : for j in range(0,4) : res[ii].append(sign*rTemp[ii][0][0,j]) else : for j in range(0,4) : res[ii][j] += sign*rTemp[ii][0][0,j] else : - print "Size problem adding results A",sign,rTemp[ii].shape + print("Size problem adding results A",sign,rTemp[ii].shape) raise SkipThisVertex() # spinor elif(rTemp[ii][0].shape[0]==4 and rTemp[ii][0].shape[1]==1 ) : if(len(res[ii])==0) : for j in range(0,4) : res[ii].append(sign*rTemp[ii][0][j,0]) else : for j in range(0,4) : res[ii][j] += sign*rTemp[ii][0][j,0] else : - print "Size problem adding results A",sign,rTemp[ii][0].shape + print("Size problem adding results A",sign,rTemp[ii][0].shape) raise SkipThisVertex() # 2 spin chains, should only be for a vertex else : for j1 in range(0,2) : for j2 in range (0,2) : val = sign*rTemp[j1][0]*rTemp[j2][1] if(len(res[3])==0) : res[2*j1+j2].append(val[0,0]) else : res[2*j1+j2][0] += val[0,0] def calculateDirac2(expr,start,end,startT,endT,sind,lind,Symbols,defns, iloc,unContracted,spins,lorentz) : tDot="" # output sVal={} # no of chains nchain=len(expr) # now deal with the uncontracted cases contracted={} # sort out contracted and uncontracted indices keys = unContracted.keys() for key in keys: # summed dummy index if key.type=="D" : contracted[key]=0 del unContracted[key] # RS index elif key.type =="R" : contracted[key]=0 del unContracted[key] # tensor index elif key.type == "T1" or key.type=="T2" : contracted[key]=0 del unContracted[key] # external index elif key.type == "O" : continue # uncontracted vector index elif key.type=="E" or key.type=="Q": continue else : - print 'Unknown type of uncontracted index',key + print('Unknown type of uncontracted index',key) raise SkipThisVertex() # check the lorentz structures for lstruct in lorentz : if(lstruct.name=="Epsilon" or lstruct.name=="Vector") : for index in lstruct.lorentz : if(index.type=="E" and index.value==iloc) : unContracted[index]=0 elif(index.type=="P" or index.type=="E" or index.type=="R" or index.type=="D") : contracted[index]=0 else : - print 'Unknown index',index, 'in ',lstruct + print('Unknown index',index, 'in ',lstruct) raise SkipThisVertex() elif(lstruct.name=="Tensor") : if(iloc==lstruct.value) : Symbols += momCom.substitute({"v": "P%s"%lstruct.value}) Symbols += "OM%s = Symbol(\"OM%s\")\n" % (lstruct.value,lstruct.value) Symbols += "M%s2 = Symbol(\"M%s2\")\n" % (lstruct.value,lstruct.value) Symbols += "p2 = Symbol(\"p2\")\n" for ival in range(1,3) : newIndex=LorentzIndex(ival) newIndex.type="O" newIndex.dimension=0 lstruct.lorentz.append(newIndex) unContracted[newIndex]=0 # contracted with self if(len(lstruct.lorentz)==0) : pass # both indices uncontracted, deal with later elif lstruct.lorentz[0].type=="T1" and lstruct.lorentz[1].type=="T2": pass elif lstruct.lorentz[0].type=="T1": pIndex = LorentzIndex(lstruct.value) pIndex.dimension=1 pIndex.type="P" (tDot,dtemp) = constructDotProduct(pIndex,lstruct.lorentz[1],defns) Symbols+="%s = Symbol(\"%s\")\n" %(tDot,tDot) Symbols += momCom.substitute({"v": lstruct.lorentz[1]}) elif lstruct.lorentz[1].type=="T2" : pIndex = LorentzIndex(lstruct.value) pIndex.dimension=1 pIndex.type="P" (tDot,dtemp) = constructDotProduct(pIndex,lstruct.lorentz[0],defns) Symbols+="%s = Symbol(\"%s\")\n" %(tDot,tDot) Symbols += momCom.substitute({"v": lstruct.lorentz[0]}) # both indices still to be contracted else : contracted[lstruct.lorentz[0].type]=0 contracted[lstruct.lorentz[1].type]=0 else : - print 'Unknown lorentz object in calculateDirac2',lstruct,iloc + print('Unknown lorentz object in calculateDirac2',lstruct,iloc) raise SkipThisVertex() # iterate over the uncontracted indices while True : # loop over the unContracted indices res = [] for i in range(0,nchain) : res.append([]) res.append([]) # loop over the contracted indices while True : # sign from metric tensor in contractions sign = 1 - for key,val in contracted.iteritems() : + for key,val in contracted.items() : if(val>0) : sign *=-1 eTemp =[] sTemp =[] fTemp =[] sTTemp=[] fTTemp=[] # make the necessary replacements for remaining indices for ichain in range(0,nchain) : # compute the main expression eTemp.append([]) for val in expr[ichain] : # already a matrix if(val.name=="M") : eTemp[ichain].append(val) # gamma(mu), replace with correct dirac matrix elif(val.name=="GMU") : if(val.index in contracted) : eTemp[ichain].append(dirac[contracted[val.index]]) elif(val.index in unContracted) : eTemp[ichain].append(dirac[unContracted[val.index]]) else : - print 'Unknown index for gamma matrix',val + print('Unknown index for gamma matrix',val) raise SkipThisVertex() # unknown to be sorted out else : - print 'Unknown type in expr',val + print('Unknown type in expr',val) raise SkipThisVertex() # start and end # start if(start[ichain].name=="S" or start[ichain].name=="M" ) : sTemp.append(start[ichain].value) elif(start[ichain].name=="RS") : if(start[ichain].index in contracted) : sTemp.append(start[ichain].value.substitute({"L" : imap[ contracted[start[ichain].index]] })) else : sTemp.append(start[ichain].value.substitute({"L" : imap[unContracted[start[ichain].index]] })) elif(start[ichain].name=="RQ") : i1 = unContracted[start[ichain].index] sTemp.append(start[ichain].value.substitute({"A" : imap[i1], "DA" : dirac[i1] })) elif(start[ichain].name=="RP") : i1 = unContracted[start[ichain].index[0]] i2 = contracted[start[ichain].index[1]] eta=0 if(i1==i2) : if(i1==0) : eta = 1 else : eta = -1 sTemp.append(start[ichain].value.substitute({"eta" : eta, "A":imap[i1] , "B":imap[i2] , "DA": dirac[i1], "DB": dirac[i2]})) else : - print 'Barred spinor not a spinor',start[ichain] + print('Barred spinor not a spinor',start[ichain]) raise SkipThisVertex() if(startT[ichain].name=="S" or startT[ichain].name=="M" ) : sTTemp.append(startT[ichain].value) elif(startT[ichain].name=="RS") : if(startT[ichain].index in contracted) : sTTemp.append(startT[ichain].value.substitute({"L" : imap[ contracted[startT[ichain].index]] })) else : sTTemp.append(startT[ichain].value.substitute({"L" : imap[unContracted[startT[ichain].index]] })) elif(startT[ichain].name=="RQ") : i1 = unContracted[startT[ichain].index] sTTemp.append(startT[ichain].value.substitute({"A" : imap[i1], "DA" : dirac[i1] })) elif(startT[ichain].name=="RP") : i1 = unContracted[startT[ichain].index[0]] i2 = contracted[startT[ichain].index[1]] eta=0 if(i1==i2) : if(i1==0) : eta = 1 else : eta = -1 sTTemp.append(startT[ichain].value.substitute({"eta" : eta, "A":imap[i1] , "B":imap[i2] , "DA": dirac[i1], "DB": dirac[i2]})) else : - print 'barred spinorT not a spinor',startT[ichain] + print('barred spinorT not a spinor',startT[ichain]) raise SkipThisVertex() # end if(end[ichain].name=="S" or end[ichain].name=="M" ) : fTemp.append(end[ichain].value) elif(end[ichain].name=="RS") : if(end[ichain].index in contracted) : fTemp.append(end[ichain].value.substitute({"L" : imap[ contracted[end[ichain].index]] })) else : fTemp.append(end[ichain].value.substitute({"L" : imap[unContracted[end[ichain].index]] })) elif(end[ichain].name=="RQ") : i1 = unContracted[end[ichain].index] fTemp.append(end[ichain].value.substitute({"B" : imap[i1], "DB": dirac[i1] })) elif(end[ichain].name=="RP") : i1 = contracted[end[ichain].index[0]] i2 = unContracted[end[ichain].index[1]] eta=0 if(i1==i2) : if(i1==0) : eta = 1 else : eta = -1 fTemp.append(end[ichain].value.substitute({"eta" : eta, "A":imap[i1] , "B":imap[i2] , "DA": dirac[i1], "DB": dirac[i2]})) else : - print 'spinor not a spinor',end[ichain] + print('spinor not a spinor',end[ichain]) raise SkipThisVertex() if(endT[ichain].name=="S" or endT[ichain].name=="M" ) : fTTemp.append(endT[ichain].value) elif(endT[ichain].name=="RS") : if(endT[ichain].index in contracted) : fTTemp.append(endT[ichain].value.substitute({"L" : imap[ contracted[endT[ichain].index]] })) else : fTTemp.append(endT[ichain].value.substitute({"L" : imap[unContracted[endT[ichain].index]] })) elif(endT[ichain].name=="RQ") : i1 = unContracted[endT[ichain].index] fTTemp.append(endT[ichain].value.substitute({"B" : imap[i1], "DB": dirac[i1] })) elif(endT[ichain].name=="RP") : i1 = contracted[endT[ichain].index[0]] i2 = unContracted[endT[ichain].index[1]] eta=0 if(i1==i2) : if(i1==0) : eta = 1 else : eta = -1 fTTemp.append(endT[ichain].value.substitute({"eta" : eta, "A":imap[i1] , "B":imap[i2] , "DA": dirac[i1], "DB": dirac[i2]})) else : - print 'spinorT not a spinor',endT[ichain] + print('spinorT not a spinor',endT[ichain]) raise SkipThisVertex() # and none dirac lorentz structures isZero = False for li in lorentz: # uncontracted vector if(li.name=="Vector") : index = unContracted[li.lorentz[0]] Symbols += momCom.substitute({"v":li.value}) for ichain in range(0,nchain) : eTemp[ichain].append("%s%s"% (li.value,imap[index]) ) elif(li.name=="Epsilon") : value="" ival=[] for index in li.lorentz : if(index in contracted) : if(index.type=="P" or index.type=="E") : value += "*%s%s" % (index,imap[contracted[index]]) ival.append(contracted[index]) elif(index.type=="R" or index.type=="D") : ival.append(contracted[index]) else : - print 'Unknown index in Epsilon Tensor',index + print('Unknown index in Epsilon Tensor',index) raise SkipThisVertex() elif(index in unContracted) : ival.append(unContracted[index]) if(len(value)!=0 and value[0]=="*") : value = value[1:] eVal = epsValue[ival[0]][ival[1]][ival[2]][ival[3]] if(eVal==0) : isZero = True else : for ichain in range(0,nchain) : eTemp[ichain].append("(%s*%s)"% (eVal,value) ) elif(li.name=="Tensor") : if(li.lorentz[0] in unContracted and li.lorentz[1] in unContracted) : value='0.5*(%s)'%tPropA[unContracted[li.lorentz[0]]][unContracted[li.lorentz[0]]].substitute({"iloc" : li.value}) for ichain in range(0,nchain) : eTemp[ichain].append("(%s)"% (value) ) elif(len(li.lorentz)==4) : if li.lorentz[0].type=="T1" and li.lorentz[1].type=="T2" : value='0.5*(%s)'%tPropC[unContracted[li.lorentz[2]]][unContracted[li.lorentz[3]]][contracted[li.lorentz[0]]][contracted[li.lorentz[1]]].substitute({"iloc" : li.value}) elif li.lorentz[0].type=="T1": value='0.5*(%s)'%tPropB[unContracted[li.lorentz[2]]][unContracted[li.lorentz[3]]][contracted[li.lorentz[0]]].substitute({"iloc" : li.value, "V" : li.lorentz[1], "dot" : tDot}) elif li.lorentz[1].type=="T2" : value= '0.5*(%s)'%tPropB[unContracted[li.lorentz[2]]][unContracted[li.lorentz[3]]][contracted[li.lorentz[1]]].substitute({"iloc" : li.value, "V" : li.lorentz[0], "dot" : tDot}) else : - print 'Both contracted tensor indices contracted' + print('Both contracted tensor indices contracted') raise SkipThisVertex() for ichain in range(0,nchain) : eTemp[ichain].append("(%s)"% (value) ) else : - print 'Uncontracted on-shell tensor' + print('Uncontracted on-shell tensor') raise SkipThisVertex() # unknown else : - print 'Unknown expression in lorentz loop',li + print('Unknown expression in lorentz loop',li) raise SkipThisVertex() # now evaluate the result if(not isZero) : rTemp =[[],[]] for ichain in range(0,nchain) : core = "*".join(str(x) for x in eTemp[ichain]) temp={} exec("import sympy\nfrom sympy import Symbol,Matrix\n"+Symbols+"result="+ ( "(%s)*(%s)*(%s)" %(sTemp[ichain],core,fTemp[ichain]))) in temp rTemp[0].append(temp["result"]) temp={} exec("import sympy\nfrom sympy import Symbol,Matrix,Transpose\n"+Symbols+"result="+ ( "(%s)*(%s)*(Transpose(%s))*(%s)*(%s)" %(sTTemp[ichain],CC,core,CD,fTTemp[ichain]))) in temp rTemp[1].append(temp["result"]) # and add it to the output addToOutput(res,nchain,sign,rTemp) #### END OF THE CONTRACTED LOOP ##### # increment the indices being summed over keys=contracted.keys() ii = len(keys)-1 while ii >=0 : if(contracted[keys[ii]]<3) : contracted[keys[ii]]+=1 break else : contracted[keys[ii]]=0 ii-=1 nZero=0 - for (key,val) in contracted.iteritems() : + for (key,val) in contracted.items() : if(val==0) : nZero+=1 if(nZero==len(contracted)) : break ###### END OF THE UNCONTRACTED LOOP ###### # no uncontracted indices if(len(unContracted)==0) : if(len(res[0])==1) : # scalar if(len(res)==2) : sVal["s" ] = res[0] sVal["sT"] = res[1] # 4 fermion else : sVal["s" ] = res[0] sVal["sT2" ] = res[1] sVal["sT1" ] = res[2] sVal["sT12"] = res[3] # spinor elif(len(res[0])==4) : for k in range(0,4) : sVal[ "s%s" % (k+1) ] = res[0][k] sVal[ "sT%s" % (k+1) ] = res[1][k] else : - print 'Sum problem',len(res),len(res[0]) + print('Sum problem',len(res),len(res[0])) raise SkipThisVertex() break # uncontracted indices else : istring = "" - for (key,val) in unContracted.iteritems() : + for (key,val) in unContracted.items() : istring +=imap[val] if(len(istring)>2) : - print 'Index problem',istring + print('Index problem',istring) raise SkipThisVertex() sVal[istring] = res[0] sVal[istring+"T"] = res[1] # increment the unsummed indices keys=unContracted.keys() ii = len(keys)-1 while ii >=0 : if(unContracted[keys[ii]]<3) : unContracted[keys[ii]]+=1 break else : unContracted[keys[ii]]=0 ii-=1 nZero=0 - for (key,val) in unContracted.iteritems() : + for (key,val) in unContracted.items() : if(val==0) : nZero+=1 if(nZero==len(unContracted)) : break # handle the vector case if( "tt" in sVal) : if(len(sVal)==32 and "tt" in sVal and len(sVal["tt"])==1) : for key in sVal: sVal[key] = sVal[key][0] else : - print 'Tensor sum problem',len(sVal) + print('Tensor sum problem',len(sVal)) raise SkipThisVertex() elif( "t" in sVal ) : # deal with pure vectors if(len(sVal)==8 and "t" in sVal and len(sVal["t"])==1) : pass # RS spinors elif(len(sVal)==8 and "t" in sVal and len(sVal["t"])==4) : pass else : - print 'Value problem',len(sVal) + print('Value problem',len(sVal)) raise SkipThisVertex() else : if("s" in sVal) : for key in sVal: sVal[key] = sVal[key][0] return sVal def convertDiracStructure(parsed,output,dimension,defns,iloc,L,lorentztag,vertex) : # get the spins of the particles spins = vertex.lorentz[0].spins # check if we have one or two spin chains nchain = (lorentztag.count("F")+lorentztag.count("R"))/2 # storage of the intermediate results expr =[] start =[] end =[] startT=[] endT =[] sind=[0]*nchain lind=[0]*nchain unContracted={} Symbols="" dtemp=[0,0,0] # parse the dirac matrix strings for ichain in range(0,nchain) : expr .append([]) start .append("") startT.append("") end .append("") endT .append("") sind[ichain],lind[ichain],expr[ichain],start[ichain],startT[ichain],end[ichain],endT[ichain],Symbols =\ processChain(dtemp,parsed,spins,Symbols,unContracted,defns,iloc) # standard ordering of the chains for ichain in range(0,nchain) : for jchain in range(ichain+1,nchain) : if(sind[jchain] s%s = sW%s.wave();" % (fIndex[i-1],fIndex[i-1])) nf+=1 else : decls.append("SpinorBarWaveFunction & sbarW%s" % (fIndex[i-1])) momenta.append([False,"Lorentz5Momentum P%s =-sbarW%s.momentum();" % (fIndex[i-1],fIndex[i-1])]) waves.append("LorentzSpinorBar sbar%s = sbarW%s.wave();" % (fIndex[i-1],fIndex[i-1])) nf+=1 elif(spin==3) : decls.append("VectorWaveFunction & vW%s" % (i)) momenta.append([False,"Lorentz5Momentum P%s =-vW%s.momentum();" % (i,i)]) waves.append("LorentzPolarizationVector E%s = vW%s.wave();" % (i,i)) elif(spin==4) : if(i%2==1) : decls.append("RSSpinorWaveFunction & RsW%s" % (i)) momenta.append([False,"Lorentz5Momentum P%s =-RsW%s.momentum();" % (i,i)]) waves.append("LorentzRSSpinor Rs%s = RsW%s.wave();" % (i,i)) nf+=1 else : decls.append("RSSpinorBarWaveFunction & RsbarW%s" % (i)) momenta.append([False,"Lorentz5Momentum P%s =-RsbarW%s.momentum();" % (i,i)]) waves.append("LorentzRSSpinorBar Rsbar%s = RsbarW%s.wave();" % (i,i)) nf+=1 elif(spin==5) : decls.append("TensorWaveFunction & tW%s" % (i)) momenta.append([False,"Lorentz5Momentum P%s =-tW%s.momentum();" % (i,i)]) waves.append("LorentzTensor T%s = tW%s.wave();" % (i,i)) else : - print 'Unknown spin',spin + print('Unknown spin',spin) raise SkipThisVertex() poff += "-P%s" % (i) # ensure unbarred spinor first ibar=-1 isp=-1 for i in range(0,len(decls)) : if(decls[i].find("Bar")>0 and ibar==-1) : ibar=i elif(decls[i].find("Spinor")>=0 and isp==-1) : isp=i if(isp!=-1 and ibar!=-1 and isp>ibar) : decls[ibar],decls[isp] = decls[isp],decls[ibar] # constrct the signature poff = ("Lorentz5Momentum P%s = " % iloc ) + poff sig="" if(iloc==0) : sig="%s evaluate(Energy2, const %s)" % (offType,", const ".join(decls)) # special for VVT vertex if(len(vertex.lorentz[0].spins)==3 and vertex.lorentz[0].spins.count(3)==2 and vertex.lorentz[0].spins.count(5)==1) : sig = sig[0:-1] + ", Energy vmass=-GeV)" else : sig="%s evaluate(Energy2, int iopt, tcPDPtr out, const %s, complex mass=-GeV, complex width=-GeV)" % (offType,", const ".join(decls)) momenta.append([True,poff+";"]) # special for VVT vertex if(len(vertex.lorentz[0].spins)==3 and vertex.lorentz[0].spins.count(3)==2 and vertex.lorentz[0].spins.count(5)==1 and vertex.lorentz[0].spins[iloc-1]==5) : sig=sig.replace("complex mass=-GeV","Energy vmass=-GeV, complex mass=-GeV") for i in range(0,len(momenta)) : momenta[i][0]=True return offType,nf,poff,sig def combineResult(res,nf,ispin,vertex) : # extract the vals and dimensions (vals,dim) = res # construct the output structure # vertex and off-shell scalars if(ispin<=1) : otype={'res':""} # spins elif(ispin==2) : otype={'s1':"",'s2':"",'s3':"",'s4':""} # vectors elif(ispin==3) : if( "t" in vals[0][1] ) : otype={'t':"",'x':"",'y':"",'z':""} else : otype={"res":""} # RS spinors elif(ispin==4) : otype={} for i1 in imap : for i in range(1,5) : otype["%ss%s"% (i1,i)]="" # off-shell tensors elif(ispin==5) : otype={} for i1 in imap : for i2 in imap : otype["%s%s"%(i1,i2)]="" else : - print 'Unknown spin',ispin + print('Unknown spin',ispin) raise SkipThisVertex() expr=[otype] for i in range(0,nf-1) : expr.append(copy.copy(otype)) # dimension for checking dimCheck=dim[0] for i in range(0,len(vals)) : # simple signs if(vals[i]=="+" or vals[i]=="-") : for ii in range(0,len(expr)) : - for(key,val) in expr[ii].iteritems() : + for(key,val) in expr[ii].items() : expr[ii][key] = expr[ii][key]+vals[i] continue # check the dimensions if(dimCheck[0]!=dim[i][0] or dimCheck[1]!=dim[i][1] or dimCheck[2]!=dim[i][2]) : - print "Dimension problem in result",i,dimCheck,dim,vertex - print vertex.lorentz + print("Dimension problem in result",i,dimCheck,dim,vertex) + print(vertex.lorentz) for j in range(0,len(vals)) : - print j,dim[j],vals[j] + print(j,dim[j],vals[j]) raise SkipThisVertex() # simplest case if(isinstance(vals[i], basestring)) : for ii in range(0,len(expr)) : - for(key,val) in expr[ii].iteritems() : + for(key,val) in expr[ii].items() : expr[ii][key] = expr[ii][key]+vals[i] continue # more complex structures pre = vals[i][0] if(pre=="(1.0)") : pre="" if(not isinstance(vals[i][1],dict)) : - print 'must be a dict here' + print('must be a dict here') raise SkipThisVertex() # tensors if("tt" in vals[i][1]) : for i1 in imap : for i2 in imap : key="%s%s"%(i1,i2) if(pre=="") : expr[0][key] += "(%s)" % vals[i][1][key] else : expr[0][key] += "%s*(%s)" % (pre,vals[i][1][key]) if(len(expr)==2) : if(pre=="") : expr[1][key] +="(%s)" % vals[i][1][key+"T"] else : expr[1][key] +="%s*(%s)" % (pre,vals[i][1][key+"T"]) # standard fermion vertex case elif(len(vals[i][1])==2 and "s" in vals[i][1] and "sT" in vals[i][1]) : if(pre=="") : expr[0]["res"] += "(%s)" % vals[i][1]["s"] expr[1]["res"] += "(%s)" % vals[i][1]["sT"] else : expr[0]["res"] += "%s*(%s)" % (pre,vals[i][1]["s"]) expr[1]["res"] += "%s*(%s)" % (pre,vals[i][1]["sT"]) # spinor case elif(len(vals[i][1])==8 and "s1" in vals[i][1]) : for jj in range(1,5) : if(pre=="") : expr[0]["s%s" % jj] += "(%s)" % vals[i][1]["s%s" % jj] expr[1]["s%s" % jj] += "(%s)" % vals[i][1]["sT%s" % jj] else : expr[0]["s%s" % jj] += "%s*(%s)" % (pre,vals[i][1]["s%s" % jj]) expr[1]["s%s" % jj] += "%s*(%s)" % (pre,vals[i][1]["sT%s" % jj]) # vector elif(len(vals[i][1])%4==0 and "t" in vals[i][1] and len(vals[i][1]["t"])==1 ) : for i1 in imap : if(pre=="") : expr[0][i1] += "(%s)" % vals[i][1][i1][0] else : expr[0][i1] += "%s*(%s)" % (pre,vals[i][1][i1][0]) if(len(expr)==2) : if(pre=="") : expr[1][i1] +="(%s)" % vals[i][1][i1+"T"][0] else : expr[1][i1] +="%s*(%s)" % (pre,vals[i][1][i1+"T"][0]) # 4 fermion vertex case elif(len(vals[i][1])==4 and "sT12" in vals[i][1]) : if(pre=="") : expr[0]["res"] += "(%s)" % vals[i][1]["s"] expr[1]["res"] += "(%s)" % vals[i][1]["sT2"] expr[2]["res"] += "(%s)" % vals[i][1]["sT1"] expr[3]["res"] += "(%s)" % vals[i][1]["sT12"] else : expr[0]["res"] += "%s*(%s)" % (pre,vals[i][1]["s"]) expr[1]["res"] += "%s*(%s)" % (pre,vals[i][1]["sT2"]) expr[2]["res"] += "(%s)" % vals[i][1]["sT1"] expr[3]["res"] += "(%s)" % vals[i][1]["sT12"] # RS spinor elif(len(vals[i][1])%4==0 and "t" in vals[i][1] and len(vals[i][1]["t"])==4 ) : for i1 in imap : for k in range(1,5) : key = "%ss%s" % (i1,k) if(pre=="") : expr[0][key] += "(%s)" % vals[i][1][i1][k-1] expr[1][key] += "(%s)" % vals[i][1][i1+"T"][k-1] else : expr[0][key] += "%s*(%s)" % (pre,vals[i][1][i1][k-1]) expr[1][key] += "%s*(%s)" % (pre,vals[i][1][i1+"T"][k-1]) else : - print 'problem with type',vals[i] + print('problem with type',vals[i]) raise SkipThisVertex() # no of particles in the vertex vDim = len(vertex.lorentz[0].spins) # compute the unit and apply it unit = computeUnit2(dimCheck,vDim) if(unit!="") : for ii in range(0,len(expr)) : - for (key,val) in expr[ii].iteritems() : + for (key,val) in expr[ii].items() : expr[ii][key] = "(%s)*(%s)" % (val,unit) return expr def headerReplace(inval) : return inval.replace("virtual","").replace("ScalarWaveFunction","").replace("SpinorWaveFunction","") \ .replace("SpinorBarWaveFunction","").replace("VectorWaveFunction","").replace("TensorWaveFunction","") \ .replace("Energy2","q2").replace("int","").replace("complex","").replace("Energy","").replace("=-GeV","") \ .replace("const &","").replace("tcPDPtr","").replace(" "," ").replace("Complex","") def combineComponents(result,offType,RS) : for i in range(0,len(result)) : - for (key,value) in result[i].iteritems() : + for (key,value) in result[i].items() : output=py2cpp(value.strip(),True) result[i][key]=output[0] # simplest case, just a value if(len(result[0])==1 and "res" in result[0]) : for i in range(0,len(result)) : result[i] = result[i]["res"] result[i]=result[i].replace("1j","ii") return # calculate the substitutions if(not isinstance(result[0],basestring)) : subs=[] for ii in range(0,len(result)) : subs.append({}) - for (key,val) in result[ii].iteritems() : + for (key,val) in result[ii].items() : subs[ii]["out%s" % key]= val # spinors if("s1" in result[0]) : stype = "LorentzSpinor" sbtype = "LorentzSpinorBar" if(offType.find("Bar")>0) : (stype,sbtype)=(sbtype,stype) subs[0]["type"] = stype result[0] = SpinorTemplate.substitute(subs[0]) subs[1]["type"] = sbtype result[1] = SpinorTemplate.substitute(subs[1]) # tensors elif("tt" in result[0]) : for ii in range(0,len(result)) : result[ii] = Template("LorentzTensor(${outxx},\n${outxy},\n${outxz},\n${outxt},\n${outyx},\n${outyy},\n${outyz},\n${outyt},\n${outzx},\n${outzy},\n${outzz},\n${outzt},\n${outtx},\n${outty},\n${outtz},\n${outtt})").substitute(subs[ii]) result[ii]=result[ii].replace("(+","(") # vectors elif("t" in result[0]) : for ii in range(0,len(result)) : result[ii] = Template("LorentzVector(${outx},\n${outy},\n${outz},\n${outt})").substitute(subs[ii]) result[ii]=result[ii].replace("(+","(") # RS spinors elif("ts1" in result[0]) : stype = "LorentzRSSpinor" sbtype = "LorentzRSSpinorBar" if(offType.find("Bar")>0) : (stype,sbtype)=(sbtype,stype) subs[0]["type"] = stype result[0] = RSSpinorTemplate.substitute(subs[0]) subs[1]["type"] = sbtype result[1] = RSSpinorTemplate.substitute(subs[1]) else : - print 'Type not implemented',result + print('Type not implemented',result) raise SkipThisVertex() for i in range(0,len(result)) : result[i]=result[i].replace("1j","ii") def generateEvaluateFunction(model,vertex,iloc,values,defns,vertexEval,cf,order) : RS = "R" in vertex.lorentz[0].name FM = "F" in vertex.lorentz[0].name # extract the start and end of the spin chains if( RS or FM ) : fIndex = vertexEval[0][0][0][2] else : fIndex=0 # first construct the signature of the function decls=[] momenta=[] waves=[] offType,nf,poff,sig = constructSignature(vertex,order,iloc,decls,momenta,waves,fIndex) # combine the different terms in the result symbols=set() localCouplings=[] result=[] ispin = 0 if(iloc!=0) : ispin = vertex.lorentz[0].spins[iloc-1] # put the lorentz structures and couplings together for j in range(0,len(vertexEval)) : # get the lorentz structure piece expr = combineResult(vertexEval[j],nf,ispin,vertex) # get the coupling for this bit val, sym = py2cpp(values[j]) localCouplings.append("Complex local_C%s = %s;\n" % (j,val)) symbols |=sym # put them together vtype="Complex" if("res" in expr[0] and offType=="VectorWaveFunction") : vtype="LorentzPolarizationVector" if(len(result)==0) : for ii in range(0,len(expr)) : result.append({}) - for (key,val) in expr[ii].iteritems() : + for (key,val) in expr[ii].items() : result[ii][key] = " %s((local_C%s)*(%s)) " % (vtype,j,val) else : for ii in range(0,len(expr)) : - for (key,val) in expr[ii].iteritems(): + for (key,val) in expr[ii].items(): result[ii][key] += " + %s((local_C%s)*(%s)) " % (vtype,j,val) # for more complex types merge the spin/lorentz components combineComponents(result,offType,RS) # multiple by scalar wavefunctions scalars="" for i in range (0,len(vertex.lorentz[0].spins)) : if(vertex.lorentz[0].spins[i]==1 and i+1!=iloc) : scalars += "sca%s*" % (i+1) if(scalars!="") : for ii in range(0,len(result)) : result[ii] = "(%s)*(%s)" % (result[ii],scalars[0:-1]) # vertex, just return the answer if(iloc==0) : result[0] = "return (%s)*(%s);\n" % (result[0],py2cpp(cf)[0]) if(FM or RS) : for i in range(1,len(result)) : result[i] = "return (%s)*(%s);\n" % (result[i],py2cpp(cf)[0]) # off-shell particle else : # off-shell scalar if(vertex.lorentz[0].spins[iloc-1] == 1 ) : result[0] = scaTemplate.format(iloc=iloc,cf=py2cpp(cf[0])[0],res=result[0]) if(FM or RS) : result[1] = scaTemplate.format(iloc=iloc,cf=py2cpp(cf[0])[0],res=result[1]) # off-shell fermion elif(vertex.lorentz[0].spins[iloc-1] == 2 ) : result[0] = sTemplate.format(iloc=iloc,cf=py2cpp(cf[0])[0],offTypeA=offType.replace("WaveFunction",""), res=result[0].replace( "M%s" % iloc, "mass" ),offTypeB=offType) if(FM or RS) : if(offType.find("Bar")>0) : offTypeT=offType.replace("Bar","") else : offTypeT=offType.replace("Spinor","SpinorBar") result[1] = sTemplate.format(iloc=iloc,cf=py2cpp(cf[0])[0],offTypeA=offTypeT.replace("WaveFunction",""), res=result[1].replace( "M%s" % iloc, "mass" ),offTypeB=offTypeT) # off-shell vector elif(vertex.lorentz[0].spins[iloc-1] == 3 ) : result[0] = vecTemplate.format(iloc=iloc,res=result[0],cf=py2cpp(cf)[0]) if(FM or RS) : result[1] = vecTemplate.format(iloc=iloc,cf=py2cpp(cf[0])[0],res=result[1]) elif(vertex.lorentz[0].spins[iloc-1] == 4 ) : if(offType.find("Bar")>0) : offTypeT=offType.replace("Bar","") else : offTypeT=offType.replace("Spinor","SpinorBar") result[1] = RSTemplate.format(iloc=iloc,cf=py2cpp(cf[0])[0],offTypeA=offTypeT.replace("WaveFunction",""), res=result[1].replace( "M%s" % iloc, "mass" ),offTypeB=offTypeT) result[0] = RSTemplate.format(iloc=iloc,cf=py2cpp(cf[0])[0],offTypeA=offType.replace("WaveFunction",""), res=result[0].replace( "M%s" % iloc, "mass" ),offTypeB=offType) # tensors elif(vertex.lorentz[0].spins[iloc-1]) : if(RS) : - print "RS spinors and tensors not handled" + print("RS spinors and tensors not handled") raise SkipThisVertex() result[0] = tenTemplate.format(iloc=iloc,cf=py2cpp(cf)[0],res=result[0]) if(FM or RS) : result[1] = tenTemplate.format(iloc=iloc,cf=py2cpp(cf)[0],res=result[1]) else : - print 'Unknown spin for off-shell particle',vertex.lorentz[0].spins[iloc-1] + print('Unknown spin for off-shell particle',vertex.lorentz[0].spins[iloc-1]) raise SkipThisVertex() # check if momenta defns needed to clean up compile of code - for (key,val) in defns.iteritems() : + for (key,val) in defns.items() : if( isinstance(key, basestring)) : if(key.find("vvP")==0) : momenta[int(key[3])-1][0] = True else : for vals in key : if(vals.type=="P") : momenta[vals.value-1][0] = True # cat the definitions defString="" - for (key,value) in defns.iteritems() : + for (key,value) in defns.items() : if(value[0]=="") : continue if(value[0][0]=="V" or value[0][0]=="T") : defString+=" %s\n" %value[1] - for (key,value) in defns.iteritems() : + for (key,value) in defns.items() : if(value[0]=="") : continue if(value[0][0]!="V" and value[0][0]!="T") : defString+=" %s\n" %value[1] if(len(result)<=2) : sorder=swapOrder(vertex,iloc,momenta,fIndex) else : sorder="" momentastring="" for i in range(0,len(momenta)) : if(momenta[i][0] and momenta[i][1]!="") : momentastring+=momenta[i][1]+"\n " # special for 4-point VVVV if(vertex.lorentz[0].spins.count(3)==4 and iloc==0) : sig=sig.replace("Energy2","Energy2,int") header="virtual %s" % sig sig=sig.replace("=-GeV","") symboldefs = [ def_from_model(model,s) for s in symbols ] function = evaluateTemplate.format(decl=sig,momenta=momentastring,defns=defString, waves="\n ".join(waves),symbols='\n '.join(symboldefs), couplings="\n ".join(localCouplings), result=result[0],swap=sorder) # special for transpose if(FM or RS) : h2=header if(not RS) : h2=header.replace("evaluate","evaluateN") function=function.replace("evaluate(","evaluateN(") headers=[] newHeader="" for ifunction in range(1,len(result)) : waveNew=[] momentastring="" htemp = header.split(",") irs=-1 isp=-1 # RS case if(RS) : # sort out the wavefunctions for i in range(0,len(waves)) : if(waves[i].find("Spinor")<0) : waveNew.append(waves[i]) continue if(waves[i].find("Bar")>0) : waveNew.append(waves[i].replace("Bar","").replace("bar","")) else : waveNew.append(waves[i].replace("Spinor","SpinorBar").replace(" s"," sbar").replace("Rs","Rsbar")) # sort out the momenta definitions for i in range(0,len(momenta)) : if(momenta[i][0] and momenta[i][1]!="") : if(momenta[i][1].find("barW")>0) : momentastring+=momenta[i][1].replace("barW","W")+"\n " elif(momenta[i][1].find("sW")>0) : momentastring+=momenta[i][1].replace("sW","sbarW")+"\n " else : momentastring+=momenta[i][1]+"\n " # header string for i in range(0,len(htemp)) : if(htemp[i].find("RS")>0) : if(htemp[i].find("Bar")>0) : htemp[i]=htemp[i].replace("Bar","").replace("RsbarW","RsW") else : htemp[i]=htemp[i].replace("Spinor","SpinorBar").replace("RsW","RsbarW") if(i>0) : irs=i elif(htemp[i].find("Spinor")>0) : if(htemp[i].find("Bar")>0) : htemp[i]=htemp[i].replace("Bar","").replace("sbarW","sW") else : htemp[i]=htemp[i].replace("Spinor","SpinorBar").replace("sW","sbarW") if(i>0) : isp=i if(irs>0 and isp >0) : htemp[irs],htemp[isp] = htemp[isp],htemp[irs] # fermion case else : htemp2 = header.split(",") # which fermions to exchange if(len(fIndex)==2) : isp = (fIndex[0],) ibar = (fIndex[1],) else : if(ifunction==1) : isp = (fIndex[2],) ibar = (fIndex[3],) elif(ifunction==2) : isp = (fIndex[0],) ibar = (fIndex[1],) elif(ifunction==3) : isp = (fIndex[0],fIndex[2]) ibar = (fIndex[1],fIndex[3]) # wavefunctions for i in range(0,len(waves)) : if(waves[i].find("Spinor")<0) : waveNew.append(waves[i]) continue if(waves[i].find("Bar")>0) : found=False for itest in range(0,len(ibar)) : if(waves[i].find("sbarW%s"%ibar[itest])>=0) : waveNew.append(waves[i].replace("Bar","").replace("bar","")) found=True break if(not found) : waveNew.append(waves[i]) else : found=False for itest in range(0,len(isp)) : if(waves[i].find("sW%s"%isp[itest])>=0) : waveNew.append(waves[i].replace("Spinor","SpinorBar").replace(" s"," sbar")) found=True break if(not found) : waveNew.append(waves[i]) # momenta definitions for i in range(0,len(momenta)) : if(momenta[i][0] and momenta[i][1]!="") : if(momenta[i][1].find("barW")>0) : found = False for itest in range(0,len(ibar)) : if(momenta[i][1].find("barW%s"%ibar[itest])>=0) : momentastring+=momenta[i][1].replace("barW","W")+"\n " found=True break if(not found) : momentastring+=momenta[i][1]+"\n " elif(momenta[i][1].find("sW")>0) : found=False for itest in range(0,len(isp)) : if(momenta[i][1].find("sW%s"%isp[itest])>=0) : momentastring+=momenta[i][1].replace("sW","sbarW")+"\n " found=True break if(not found) : momentastring+=momenta[i][1]+"\n " else : momentastring+=momenta[i][1]+"\n " # header for i in range(0,len(htemp)) : if(htemp[i].find("Spinor")<0) : continue if(htemp[i].find("Bar")>0) : if(i==0) : htemp[i] =htemp [i].replace("Bar","") htemp2[i]=htemp2[i].replace("Bar","") continue for itest in range(0,len(ibar)) : if(htemp[i].find("sbarW%s"%ibar[itest])>=0) : htemp[i] =htemp [i].replace("Bar","").replace("sbarW","sW") htemp2[i]=htemp2[i].replace("Bar","").replace("sbarW%s"%ibar[itest],"sW%s"%isp[itest]) break else : if(i==0) : htemp [i]=htemp [i].replace("Spinor","SpinorBar") htemp2[i]=htemp2[i].replace("Spinor","SpinorBar") continue for itest in range(0,len(isp)) : if(htemp[i].find("sW%s"%isp[itest])>=0) : htemp [i]=htemp [i].replace("Spinor","SpinorBar").replace("sW","sbarW") htemp2[i]=htemp2[i].replace("Spinor","SpinorBar").replace("sW%s"%isp[itest],"sbarW%s"%ibar[itest]) break # header for transposed function hnew = ','.join(htemp) hnew = hnew.replace("virtual ","").replace("=-GeV","") if(not RS) : theader = ','.join(htemp2) theader = theader.replace("virtual ","").replace("=-GeV","") if(len(result)==2) : hnew =hnew .replace("evaluate","evaluateT") theader=theader.replace("evaluate","evaluateT") else : hnew =hnew .replace("evaluate","evaluateT%s" % ifunction) theader=theader.replace("evaluate","evaluateT%s" % ifunction) if(iloc not in fIndex) : theader = headerReplace(theader) else : theader = headerReplace(h2).replace("evaluateN","evaluateT") headers.append(theader) newHeader += hnew +";\n" else : newHeader += hnew fnew = evaluateTemplate.format(decl=hnew,momenta=momentastring,defns=defString, waves="\n ".join(waveNew),symbols='\n '.join(symboldefs), couplings="\n ".join(localCouplings), result=result[ifunction],swap=sorder) function +="\n" + fnew if(FM and not RS) : if(len(result)==2) : if(iloc!=fIndex[1]) : fi=1 stype="sbar" else : fi=0 stype="s" header = vTemplateT.format(header=header.replace("Energy2,","Energy2 q2,"), normal=headerReplace(h2), transpose=theader,type=stype, iloc=fIndex[fi],id=vertex.particles[fIndex[fi]-1].pdg_code) \ +newHeader+h2.replace("virtual","") else : sorder = swapOrderFFFF(vertex,iloc,fIndex) header = vTemplate4.format(header=header.replace("Energy2,","Energy2 q2,"), iloc1=fIndex[1],iloc2=fIndex[3],swap=sorder, id1=vertex.particles[fIndex[1]-1].pdg_code, id2=vertex.particles[fIndex[3]-1].pdg_code, cf=py2cpp(cf)[0], res1=headerReplace(h2).replace("W",""), res2=headers[0].replace("W",""), res3=headers[1].replace("W",""), res4=headers[2].replace("W",""))\ +newHeader+h2.replace("virtual","") else : header=header + ";\n" + newHeader return (header,function) evaluateMultiple = """\ {decl} {{ {code} }} """ def multipleEvaluate(vertex,spin,defns) : if(spin==1) : name="scaW" elif(spin==3) : name="vW" else : - print 'Evaluate multiple problem',spin + print('Evaluate multiple problem',spin) raise SkipThisVertex() if(len(defns)==0) : return ("","") header = defns[0] ccdefn = header.replace("=-GeV","").replace("virtual ","").replace("Energy2","Energy2 q2") code="" spins=vertex.lorentz[0].spins iloc=1 waves=[] for i in range(0,len(spins)) : if(spins[i]==spin) : waves.append("%s%s" %(name,i+1)) for i in range(0,len(spins)) : if(spins[i]==spin) : if(iloc==1) : el="" else : el="else " call = headerReplace(defns[iloc]) if(iloc!=1) : call = call.replace(waves[0],waves[iloc-1]) pdgid = vertex.particles[i].pdg_code code += " %sif(out->id()==%s) return %s;\n" % (el,pdgid,call) iloc+=1 code+=" else assert(false);\n" return (header,evaluateMultiple.format(decl=ccdefn,code=code)) 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,266 +1,312 @@ +from __future__ import print_function from string import Template from os import path -import sys,cmath +import sys,cmath,glob 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 = '{name}.template'.format(name=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 coupling_orders(vertex, coupling, defns): # if more than one take QCD over QED and then lowest order in QED if type(coupling) is list: - print 'not sure this happens' + print('not sure this happens') quit() 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: output={} for ctype in defns : output[ctype]=coupling.order.get(ctype,0) return output def def_from_model(FR,s): """Return a C++ line that defines parameter s as coming from the model file.""" if("hw_kine" in s) :return "" 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'({symb})(\W|$)'.format(symb=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] class SkipThisVertex(Exception): pass def extractAntiSymmetricIndices(instring,funct) : terms = instring.strip(funct).strip(")").split(",") sign=1. for iy in range(0,len(terms)) : for ix in range(-1,-len(terms)+iy,-1) : swap = False if(len(terms[ix])==1 and len(terms[ix-1])==1) : swap = int(terms[ix]) items + line=line.replace("iteritems","items") + # fix imports + if("import" in line) : + for val in names : + if("import %s" %val in line) : + line=line.replace("import %s" %val, + "from . import %s" %val) + if(val in tNames) : tNames.remove(val) + if("from %s" %val in line) : + line=line.replace("from %s" %val, + "from .%s" %val) + # add brackets to print statements + if("print" in line) : + if line.strip()[0:5] == "print" : + line=line.strip("\n").replace("print","print(")+")\n" + output += line + line=inFile.readline() + inFile.close() + if(isInit) : + if "__init__" in tNames : tNames.remove("__init__") + temp="" + for val in tNames : + temp+= "from . import %s\n" % val + output = temp+output + with open(fName,'w') as dest: + dest.write(output) diff --git a/Models/Feynrules/python/ufo2peg/particles.py b/Models/Feynrules/python/ufo2peg/particles.py --- a/Models/Feynrules/python/ufo2peg/particles.py +++ b/Models/Feynrules/python/ufo2peg/particles.py @@ -1,243 +1,243 @@ +from __future__ import print_function from string import Template - # ignore these, they're in Hw++ already # TODO reset Hw++ settings instead SMPARTICLES = { 1:'d', 2:'u', 3:'s', 4:'c', 5:'b', 6:'t', # think about this one later 11:'e-', 12:'nu_e', 13:'mu-', 14:'nu_mu', 15:'tau-', 16:'nu_tau', 21:'g', 22:'gamma', 23:'Z0', 24:'W+', -1:'dbar', -2:'ubar', -3:'sbar', -4:'cbar', -5:'bbar', -6:'tbar', -11:'e+', -12:'nu_ebar', -13:'mu+', -14:'nu_mubar', -15:'tau+', -16:'nu_taubar', -24:'W-', } particleT = Template( """ create ThePEG::ParticleData $name # values set to 999999 are recalculated later from other model parameters setup $name $pdg_code $name $mass $width $wcut $ctau $charge $color $spin 0 """ ) class ParticleConverter: 'Convert a FR particle to extract the information ThePEG needs.' def __init__(self,p,parmsubs,modelparameters): self.name = p.name self.pdg_code = p.pdg_code self.spin = p.spin self.color = p.color if self.color == 1: self.color = 0 self.selfconjugate = 0 self.mass = parmsubs[str(p.mass)] if type(self.mass) == str: value = modelparameters[self.mass] try: value = value.real except: pass newname = '%s_ABS' % self.mass self.mass = '${%s}' % newname modelparameters[newname] = abs(value) else: try: self.mass = self.mass.real except: pass self.mass = 999999. # abs(self.mass) hbarc = 197.3269631e-15 # GeV mm (I hope ;-) ) self.width = parmsubs[str(p.width)] if type(self.width) == str: width = modelparameters[self.width] ctau = (hbarc / width) if width != 0 else 0 newname = '%s_CTAU' % self.width self.ctau = '${%s}' % newname modelparameters[newname] = ctau wcut = 10 * width newname = '%s_WCUT' % self.width self.wcut = '${%s}' % newname modelparameters[newname] = wcut self.width = '${%s}' % self.width else: self.ctau = 999999. # (hbarc / self.width) if self.width != 0 else 0 self.wcut = 999999. #10.0 * self.width self.width = 999999. # was blank line before self.charge = int(3 * p.charge) def subs(self): return self.__dict__ def check_effective_vertex(FR,p,ig) : for vertex in FR.all_vertices: if(len(vertex.particles) != 3) : continue if(p not in vertex.particles ) : continue ng=0 for part in vertex.particles : if(part.pdg_code==ig) : ng+=1 if(ng==2) : return False return True def thepeg_particles(FR,parameters,modelname,modelparameters,forbidden_names,hw_higgs): plist = [] antis = {} names = [] splittings = [] done_splitting_QCD=[] done_splitting_QED=[] for p in FR.all_particles: if p.spin == -1: continue gsnames = ['goldstone', 'goldstoneboson', 'GoldstoneBoson'] def gstest(name): try: return getattr(p,name) except AttributeError: return False if any(map(gstest, gsnames)): continue if p.pdg_code in SMPARTICLES: continue if p.pdg_code == 25 and not hw_higgs: plist.append( """ set /Herwig/Particles/h0:Mass_generator NULL set /Herwig/Particles/h0:Width_generator NULL rm /Herwig/Masses/HiggsMass rm /Herwig/Widths/hWidth """ ) if p.name in forbidden_names: - print 'RENAMING PARTICLE',p.name,'as ',p.name+'_UFO' + print('RENAMING PARTICLE',p.name,'as ',p.name+'_UFO') p.name +="_UFO" subs = ParticleConverter(p,parameters,modelparameters).subs() if not (p.pdg_code == 25 and hw_higgs) : plist.append( particleT.substitute(subs) ) pdg, name = subs['pdg_code'], subs['name'] names.append(name) if -pdg in antis: plist.append( 'makeanti %s %s\n' % (antis[-pdg], name) ) elif not (p.pdg_code == 25 and hw_higgs) : plist.append( 'insert /Herwig/NewPhysics/NewModel:DecayParticles 0 %s\n' % name ) plist.append( 'insert /Herwig/Shower/ShowerHandler:DecayInShower 0 %s # %s' % (abs(pdg), name) ) antis[pdg] = name selfconjugate = 1 class SkipMe(Exception): pass def spin_name(s): spins = { 1 : 'Zero', 2 : 'Half', 3 : 'One' } if s not in spins: raise SkipMe() else: return spins[s] def col_name(c): cols = { 3 : 'Triplet', 6 : 'Sextet', 8 : 'Octet' } return cols[c] try: # QCD splitting functions if p.color in [3,6,8] and abs(pdg) not in done_splitting_QCD: # which colors? done_splitting_QCD.append(abs(pdg)) splitname = '{name}SplitFnQCD'.format(name=p.name) sudname = '{name}SudakovQCD'.format(name=p.name) splittings.append( """ create Herwig::{s}{s}OneSplitFn {name} set {name}:InteractionType QCD set {name}:ColourStructure {c}{c}Octet cp /Herwig/Shower/SudakovCommon {sudname} set {sudname}:SplittingFunction {name} do /Herwig/Shower/SplittingGenerator:AddFinalSplitting {pname}->{pname},g; {sudname} """.format(s=spin_name(p.spin), name=splitname, c=col_name(p.color), pname=p.name, sudname=sudname) ) except SkipMe: pass # QED splitting functions try: if p.charge != 0 and abs(pdg) not in done_splitting_QED: done_splitting_QED.append(abs(pdg)) splitname = '{name}SplitFnQED'.format(name=p.name) sudname = '{name}SudakovQED'.format(name=p.name) splittings.append( """ create Herwig::{s}{s}OneSplitFn {name} set {name}:InteractionType QED set {name}:ColourStructure ChargedChargedNeutral cp /Herwig/Shower/SudakovCommon {sudname} set {sudname}:SplittingFunction {name} set {sudname}:Alpha /Herwig/Shower/AlphaQED do /Herwig/Shower/SplittingGenerator:AddFinalSplitting {pname}->{pname},gamma; {sudname} """.format(s=spin_name(p.spin), name=splitname, pname=p.name, sudname=sudname) ) except SkipMe: pass if p.charge == 0 and p.color == 1 and p.spin == 1 and not (p.pdg_code == 25 and hw_higgs) : if(check_effective_vertex(FR,p,21)) : plist.append( """ insert /Herwig/{ModelName}/V_GenericHGG:Bosons 0 {pname} """.format(pname=p.name, ModelName=modelname) ) if(check_effective_vertex(FR,p,22)) : plist.append( """ insert /Herwig/{ModelName}/V_GenericHPP:Bosons 0 {pname} """.format(pname=p.name, ModelName=modelname) ) return ''.join(plist)+''.join(splittings), names diff --git a/Models/Feynrules/python/ufo2peg/vertices.py b/Models/Feynrules/python/ufo2peg/vertices.py --- a/Models/Feynrules/python/ufo2peg/vertices.py +++ b/Models/Feynrules/python/ufo2peg/vertices.py @@ -1,894 +1,889 @@ +from __future__ import print_function import sys,pprint from .helpers import CheckUnique,getTemplate,writeFile,coupling_orders,def_from_model from .converter import py2cpp from .collapse_vertices import collapse_vertices from .check_lorentz import tensorCouplings,VVVordering,lorentzScalar,\ processTensorCouplings,scalarCouplings,processScalarCouplings,scalarVectorCouplings,\ processScalarVectorCouplings,vectorCouplings,processVectorCouplings,fermionCouplings,processFermionCouplings,\ RSCouplings from .general_lorentz import convertLorentz,generateEvaluateFunction,multipleEvaluate from .helpers import SkipThisVertex,extractAntiSymmetricIndices,isGoldstone # prefactors for vertices lfactors = { 'FFV' : '-complex(0,1)', # ok 'VVV' : 'complex(0,1)', # changed to fix ttbar 'VVVS' : 'complex(0,1)', # should be as VVV '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)', 'RFS' : 'complex(0,1)', 'RFV' : 'complex(0,1)', } genericVertices=['FFFF','FFVV','FFSS','FFVS','VVVV','VVVT','FFVT', 'RFVV','RFVS','RFSS','SSST','VVST','FFST'] skipped5Point=False # template for the header for a vertex VERTEXHEADER = """\ #include "ThePEG/Helicity/Vertex/{spindirectory}/{lorentztag}Vertex.h" """ GENERALVERTEXHEADER = """\ #include "ThePEG/Helicity/Vertex/Abstract{lorentztag}Vertex.h" """ # template for the implmentation for a vertex VERTEXCLASS = getTemplate('Vertex_class') GENERALVERTEXCLASS = getTemplate('GeneralVertex_class') # template for the .cc file for vertices VERTEX = getTemplate('Vertex.cc') vertexline = """\ create Herwig::FRModel{classname} /Herwig/{modelname}/{classname} insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/{classname} """ def get_lorentztag(spin): """Produce a ThePEG spin tag for the given numeric FR spins.""" spins = { 1 : 'S', 2 : 'F', 3 : 'V', 4 : 'R', 5 : 'T', -1 : 'U' } result=[] for i in range(0,len(spin)) : result.append((spins[spin[i]],i+1)) - def spinsort(a,b): - """Helper function for ThePEG's FVST spin tag ordering.""" - (a1,a2) = a - (b1,b2) = b - if a1 == b1: return 0 - for letter in 'URFVST': - if a1 == letter: return -1 - if b1 == letter: return 1 + def spinsort(inVal): + output=0 + vals=['U','R','F','V','S','T'] + for val in vals: + if inVal==val : break + output+=1 + return output - result = sorted(result, cmp=spinsort) + result = sorted(result, key=spinsort) order=[] output="" for i in range(0,len(result)) : (a,b) = result[i] order.append(b) output+=a return (output,order) def unique_lorentztag(vertex): """Check and return the Lorentz tag of the vertex.""" unique = CheckUnique() for l in vertex.lorentz: (lorentztag,order) = get_lorentztag(l.spins) unique( lorentztag ) lname = l.name[:len(lorentztag)] if sorted(lorentztag) != sorted(lname): raise Exception("Lorentztags: %s is not %s in %s" % (lorentztag,lname,vertex)) return (lorentztag,order) 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 coloursort(a,b) : - if a == b: return 0 - i1=int(a[4]) - i2=int(b[4]) - if(i1==i2) : return 0 - elif(i1id()=={id1} && p2->id()=={id3} && p3->id()=={id2}) || (p1->id()=={id2} && p2->id()=={id1} && p3->id()=={id3}) || (p1->id()=={id3} && p2->id()=={id2} && p3->id()=={id1})) {{ sign *= -1.; }} norm(norm()*sign); """ if(not "p1" in couplingptrs[0]) : couplingptrs[0] += ' p1' if(not "p2" in couplingptrs[1]) : couplingptrs[1] += ' p2' if(not "p3" in couplingptrs[2]) : couplingptrs[2] += ' p3' if("Bar" not in vertex.color[0]) : order,sign = extractAntiSymmetricIndices(vertex.color[0],"Epsilon(") else : order,sign = extractAntiSymmetricIndices(vertex.color[0],"EpsilonBar(") subs = {"id1" : vertex.particles[int(order[0])-1].pdg_code, "id2" : vertex.particles[int(order[1])-1].pdg_code, "id3" : vertex.particles[int(order[2])-1].pdg_code, "epssign" : sign } append+=EPSSIGN.format(**subs) return couplingptrs,append class VertexConverter: 'Convert the vertices in a FR model to extract the information ThePEG needs.' def __init__(self,model,parmsubs,defns) : 'Initialize the parameters' self.verbose=False self.vertex_skipped=False self.ignore_skipped=False self.model=model self.all_vertices= [] self.vertex_names = {} self.modelname="" self.no_generic_loop_vertices = False self.parmsubs = parmsubs self.couplingDefns = defns self.genericTensors = False self.hw_higgs = False def readArgs(self,args) : 'Extract the relevant command line arguments' self.ignore_skipped = args.ignore_skipped self.verbose = args.verbose self.modelname = args.name self.no_generic_loop_vertices = args.no_generic_loop_vertices self.include_generic = args.include_generic self.genericTensors = args.use_generic_for_tensors self.hw_higgs = args.use_Herwig_Higgs def should_print(self) : 'Check if we should output the results' return not self.vertex_skipped or self.ignore_skipped def convert(self) : 'Convert the vertices' if(self.verbose) : - print 'verbose mode on: printing all vertices' - print '-'*60 + print('verbose mode on: printing all vertices') + print('-'*60) labels = ('vertex', 'particles', 'Lorentz', 'C_L', 'C_R', 'norm') pprint.pprint(labels) # extract the vertices self.all_vertices = self.model.all_vertices # find the SM Higgs boson higgs=None if(self.hw_higgs) : for part in self.model.all_particles: if(part.pdg_code==25) : higgs=part break # convert the vertices vertexclasses, vertexheaders = [], set() ifile=1 icount=0 for vertexnumber,vertex in enumerate(self.all_vertices,1) : # process the vertex (skip,vertexClass,vertexHeader) = \ self.processVertex(vertexnumber,vertex) # check it can be handled if(skip) : continue # check if Higgs and skip if using Hw higgs sector if higgs in vertex.particles : nH = vertex.particles.count(higgs) # skip trilinear and quartic higgs vertices if ( nH == len(vertex.particles) ) : vertex.herwig_skip_vertex = True continue elif (len(vertex.particles)-nH==2) : p=[] for part in vertex.particles : if(part!=higgs) : p.append(part) if(nH==1 and p[0].pdg_code==-p[1].pdg_code and abs(p[0].pdg_code) in [1,2,3,4,5,6,11,13,15,24]) : vertex.herwig_skip_vertex = True continue elif((nH==1 or nH==2) and p[0].pdg_code==p[1].pdg_code and p[0].pdg_code ==23) : vertex.herwig_skip_vertex = True continue elif(nH==2 and p[0].pdg_code==-p[1].pdg_code and abs(p[0].pdg_code)==24) : vertex.herwig_skip_vertex = True continue # add to the list icount +=1 vertexclasses.append(vertexClass) vertexheaders.add(vertexHeader) WRAP = 25 if icount % WRAP == 0 or vertexHeader.find("Abstract")>=0: write_vertex_file({'vertexnumber' : ifile, 'vertexclasses' : '\n'.join(vertexclasses), 'vertexheaders' : ''.join(vertexheaders), 'ModelName' : self.modelname}) vertexclasses = [] vertexheaders = set() icount=0 ifile+=1 # exit if there's vertices we can't handle if not self.should_print(): sys.stderr.write( """ Error: The conversion was unsuccessful, some vertices could not be generated. The new --include-generic option should be able to generate these. Otherwise, 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 still stuff to output it do it if vertexclasses: write_vertex_file({'vertexnumber' : ifile + 1, 'vertexclasses' : '\n'.join(vertexclasses), 'vertexheaders' : ''.join(vertexheaders), 'ModelName' : self.modelname}) - print '='*60 + print('='*60) def setCouplingPtrs(self,lorentztag,qcd,append,prepend) : couplingptrs = [',tcPDPtr']*len(lorentztag) if lorentztag == 'VSS': couplingptrs[1] += ' p2' elif lorentztag == 'FFV': couplingptrs[0] += ' p1' elif (lorentztag == 'VVV' or lorentztag == 'VVVS' or lorentztag == "SSS" or lorentztag == "VVVT" ) \ and (append or prepend ) : couplingptrs[0] += ' p1' couplingptrs[1] += ' p2' couplingptrs[2] += ' p3' elif (lorentztag == 'VVVV' and (qcd < 2 or append)) or\ (lorentztag == "SSSS" and prepend ): couplingptrs[0] += ' p1' couplingptrs[1] += ' p2' couplingptrs[2] += ' p3' couplingptrs[3] += ' p4' return couplingptrs def processVertex(self,vertexnumber,vertex) : global skipped5Point # get the Lorentz tag for the vertex lorentztag,order = unique_lorentztag(vertex) # check if we should skip the vertex vertex.herwig_skip_vertex = checkGhostGoldstoneVertex(lorentztag,vertex) # check the order of the vertex and skip 5 points if(len(lorentztag)>=5) : vertex.herwig_skip_vertex = True if(not skipped5Point) : skipped5Point = True - print "Skipping 5 point vertices which aren\'t used in Herwig7" + print("Skipping 5 point vertices which aren\'t used in Herwig7") if(vertex.herwig_skip_vertex) : return (True,"","") # check if we support this at all if( lorentztag not in lfactors and lorentztag not in genericVertices) : 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 self.vertex_skipped=True return (True,"","") # get the factor for the vertex generic = False try: lf = lfactors[lorentztag] if( self.genericTensors and "T" in lorentztag) : raise KeyError() except KeyError: if(not self.include_generic) : 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 self.vertex_skipped=True return (True,"","") else : lf=1. generic=True # get the ids of the particles at the vertex plistarray = [ ','.join([ str(vertex.particles[o-1].pdg_code) for o in order ]) ] # parse the colour structure for the vertex try: L,pos = colors(vertex) cs,cf = colorfactor(vertex,L,pos,lorentztag) except SkipThisVertex: msg = 'Warning: Color structure for vertex ( {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 self.vertex_skipped=True return (True,"","") ### classname classname = 'V_%s' % vertex.name if(not generic) : try : return self.extractGeneric(vertex,order,lorentztag,classname,plistarray,pos,lf,cf,cs) except SkipThisVertex: if(not self.include_generic) : msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \ 'is not supported, may have a non-perturbative form.\n'.format(tag=lorentztag, name=vertex.name, ps=' '.join(map(str,vertex.particles))) sys.stderr.write(msg) vertex.herwig_skip_vertex = True self.vertex_skipped=True return (True,"","") else : try : return self.extractGeneral(vertex,order,lorentztag,classname,plistarray,pos,cf,cs) except SkipThisVertex: msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \ 'is not supported, may have a non-perturbative form.\n'.format(tag=lorentztag, name=vertex.name, ps=' '.join(map(str,vertex.particles))) sys.stderr.write(msg) vertex.herwig_skip_vertex = True self.vertex_skipped=True return (True,"","") else : try : return self.extractGeneral(vertex,order,lorentztag,classname,plistarray,pos,cf,cs) except SkipThisVertex: msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \ 'is not supported, may have a non-perturbative form.\n'.format(tag=lorentztag, name=vertex.name, ps=' '.join(map(str,vertex.particles))) sys.stderr.write(msg) vertex.herwig_skip_vertex = True self.vertex_skipped=True return (True,"","") def extractGeneric(self,vertex,order,lorentztag,classname,plistarray,pos,lf,cf,cs) : classes="" headers="" # identify the maximum colour flow and orders of the couplings maxColour=0 couplingOrders=[] self.vertex_names[vertex.name] = [classname] - for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems(): + for (color_idx,lorentz_idx),coupling in vertex.couplings.items(): maxColour=max(maxColour,color_idx) orders = coupling_orders(vertex, coupling, self.couplingDefns) if(orders not in couplingOrders) : couplingOrders.append(orders) # loop the order of the couplings iorder = 0 for corder in couplingOrders : iorder +=1 cname=classname if(iorder!=1) : cname= "%s_%s" % (classname,iorder) self.vertex_names[vertex.name].append(cname) header = "" prepend="" append="" all_couplings=[] for ix in range(0,maxColour+1) : all_couplings.append([]) # loop over the colour structures for colour in range(0,maxColour+1) : - for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems() : + for (color_idx,lorentz_idx),coupling in vertex.couplings.items() : # check colour structure and coupling order if(color_idx!=colour) : continue if(coupling_orders(vertex, coupling, self.couplingDefns)!=corder) : continue # get the prefactor for the lorentz structure L = vertex.lorentz[lorentz_idx] prefactors = calculatePrefactor(lf,cf[color_idx]) # calculate the value of the coupling value = couplingValue(coupling) # handling of the different types of couplings if lorentztag in ['FFS','FFV']: all_couplings[color_idx] = fermionCouplings(value,prefactors,L,all_couplings[color_idx],order) elif 'T' in lorentztag : append, all_couplings[color_idx] = tensorCouplings(vertex,value,prefactors,L,lorentztag,pos, all_couplings[color_idx],order) elif 'R' in lorentztag : all_couplings[color_idx] = RSCouplings(value,prefactors,L,all_couplings[color_idx],order) elif lorentztag == 'VVS' or lorentztag == "VVSS" or lorentztag == "VSS" : all_couplings[color_idx] = scalarVectorCouplings(value,prefactors,L,lorentztag, all_couplings[color_idx],order) elif lorentztag == "SSS" or lorentztag == "SSSS" : prepend, header, all_couplings[color_idx] = scalarCouplings(vertex,value,prefactors,L,lorentztag, all_couplings[color_idx],prepend,header) elif "VVV" in lorentztag : all_couplings[color_idx],append = vectorCouplings(vertex,value,prefactors,L,lorentztag,pos, all_couplings[color_idx],append,corder["QCD"],order) else: raise SkipThisVertex() # final processing of the couplings symbols = set() if(lorentztag in ['FFS','FFV']) : (normcontent,leftcontent,rightcontent,append) = processFermionCouplings(lorentztag,vertex, self.model,self.parmsubs, all_couplings,order) elif('T' in lorentztag) : (leftcontent,rightcontent,normcontent) = processTensorCouplings(lorentztag,vertex,self.model, self.parmsubs,all_couplings,order) elif(lorentztag=="SSS" or lorentztag=="SSSS") : normcontent = processScalarCouplings(self.model,self.parmsubs,all_couplings) elif(lorentztag=="VVS" or lorentztag =="VVSS" or lorentztag=="VSS") : normcontent,append,lorentztag,header,sym = processScalarVectorCouplings(lorentztag,vertex, self.model,self.parmsubs, all_couplings,header,order) symbols |=sym elif("VVV" in lorentztag) : normcontent,append,header =\ processVectorCouplings(lorentztag,vertex,self.model,self.parmsubs,all_couplings,append,header) else : SkipThisVertex() # set the coupling ptrs in the setCoupling call couplingptrs = self.setCouplingPtrs(lorentztag.replace("General",""), corder["QCD"],append != '',prepend != '') ### do we need left/right? 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 = 'Complex((%s) * GeV / UnitRemoval::E)' % normcalc elif lorentztag in ['GeneralVVS']: normcalc = 'Complex(-(%s) * UnitRemoval::E / GeV )' % normcalc elif lorentztag in ['FFT','VVT', 'SST', 'FFVT', 'VVVT' , 'VVVS' ]: normcalc = 'Complex((%s) / GeV * UnitRemoval::E)' % normcalc norm = 'norm(' + normcalc + ');' # finally special handling for eps tensors if(len(vertex.color)==1 and vertex.color[0].find("Epsilon")>=0) : couplingptrs, append = epsilonSign(vertex,couplingptrs,append) # define unkown symbols from the model symboldefs = [ def_from_model(self.model,s) for s in symbols ] couplingOrder="" - for coupName,coupVal in corder.iteritems() : + for coupName,coupVal in corder.items() : couplingOrder+=" orderInCoupling(CouplingType::%s,%s);\n" %(coupName,coupVal) ### assemble dictionary and fill template subs = { 'lorentztag' : lorentztag, # ok 'classname' : cname, # 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 'addToPlist' : '\n'.join([ 'addToList(%s);'%s for s in plistarray]), 'parameters' : '', 'couplingOrders' : couplingOrder, 'colourStructure' : cs, 'couplingptrs' : ''.join(couplingptrs), 'spindirectory' : spindirectory(lorentztag), 'ModelName' : self.modelname, 'prepend' : prepend, 'append' : append, 'header' : header } # ok # print info if required if self.verbose: - print '-'*60 + print('-'*60) pprint.pprint(( classname, plistarray, leftcalc, rightcalc, normcalc )) headers+=VERTEXHEADER.format(**subs) classes+=VERTEXCLASS.substitute(subs) return (False,classes,headers) def extractGeneral(self,vertex,order,lorentztag,classname,plistarray,pos,cf,cs) : eps=False classes="" headers="" # check the colour flows, three cases supported either 1 flow or 3 in gggg # or multiple wierd ones in FFFF cidx=-1 gluon4point = (len(pos[8])==4 and vertex.lorentz[0].spins.count(3)==4) FFFF = (len(pos[3])==2 and len(pos[-3])==2 and vertex.lorentz[0].spins.count(2)==4) couplingOrders=[] colours={} - for (color_idx,lorentz_idx),coupling in vertex.couplings.iteritems() : + for (color_idx,lorentz_idx),coupling in vertex.couplings.items() : orders = coupling_orders(vertex, coupling, self.couplingDefns) if(orders not in couplingOrders) : couplingOrders.append(orders) if(gluon4point) : color = vertex.color[color_idx] f = color.split("*") (o1,s1) = extractAntiSymmetricIndices(f[0],"f(") (o2,s2) = extractAntiSymmetricIndices(f[1],"f(") if(o2[0]1 and i!=2) : mult[i] = [] for i in range(0,imax) : (evalHeader,evalCC) = generateEvaluateFunction(self.model,vertex,i,values,defns[i],vertexEval[i],cfactor,order) if(i!=0 and spins[i-1] in mult) : if(len(mult[spins[i-1]])==0) : mult[spins[i-1]].append(evalHeader) evalHeader=evalHeader.replace("evaluate(","evaluate%s(" % i) evalCC =evalCC .replace("evaluate(","evaluate%s(" % i) mult[spins[i-1]].append(evalHeader) header+=" "+evalHeader+";\n" impls+=evalCC # combine the multiple defn if needed - for (key,val) in mult.iteritems() : + for (key,val) in mult.items() : (evalHeader,evalCC) = multipleEvaluate(vertex,key,val) if(evalHeader!="") : header += " "+evalHeader+";\n" if(evalCC!="") : impls += evalCC impls=impls.replace("evaluate", "FRModel%s::evaluate" % cname) couplingOrder="" - for coupName,coupVal in corder.iteritems() : + for coupName,coupVal in corder.items() : couplingOrder+=" orderInCoupling(CouplingType::%s,%s);\n" %(coupName,coupVal) ### assemble dictionary and fill template subs = { 'lorentztag' : lorentztag, 'classname' : cname, 'addToPlist' : '\n'.join([ 'addToList(%s);'%s for s in plistarray]), 'ModelName' : self.modelname, 'couplingOrders' : couplingOrder, 'colourStructure' : cstruct, 'evaldefs' : header, 'evalimpls' : impls} newHeader = GENERALVERTEXHEADER.format(**subs) if(eps) : newHeader +="#include \"ThePEG/Helicity/epsilon.h\"\n" headers+=newHeader classes+=GENERALVERTEXCLASS.substitute(subs) return (False,classes,headers) def get_vertices(self,libname): vlist = ['library %s\n' % libname] for v in self.all_vertices: if v.herwig_skip_vertex: continue for name in self.vertex_names[v.name] : vlist.append( vertexline.format(modelname=self.modelname, classname=name) ) if( not self.no_generic_loop_vertices and not self.hw_higgs ) : vlist.append('insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_GenericHPP\n'.format(modelname=self.modelname) ) vlist.append('insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_GenericHGG\n'.format(modelname=self.modelname) ) # add Hw higgs vertices if required if(self.hw_higgs) : for vertex in ["FFH","HGG","HHH","HPP","HZP","WWHH","WWH"]: vlist.append('insert %s:ExtraVertices 0 /Herwig/Vertices/%sVertex\n' % (self.modelname,vertex) ) return ''.join(vlist)