diff --git a/Models/Feynrules/python/ufo2herwig b/Models/Feynrules/python/ufo2herwig --- a/Models/Feynrules/python/ufo2herwig +++ b/Models/Feynrules/python/ufo2herwig @@ -1,475 +1,478 @@ #! /usr/bin/env python from __future__ import division from __future__ import print_function import os, sys, pprint, argparse, re,copy from string import Template # add path to the ufo conversion modules modulepath = os.path.join("@PKGLIBDIR@",'python') sys.path.append(modulepath) from ufo2peg import * -# set up the option parser for command line input +# 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', + 'ufodir', + metavar='UFO_directory', help='the UFO model directory' ) parser.add_argument( - '-v', '--verbose', - action="store_true", + '-v', '--verbose', + action="store_true", help="print verbose output" ) parser.add_argument( - '-n','--name', - default="FRModel", + '-n','--name', + default="FRModel", help="set custom nametag for the model" ) parser.add_argument( - '--ignore-skipped', - action="store_true", + '--ignore-skipped', + action="store_true", help="ignore skipped vertices and produce output anyway" ) parser.add_argument( - '--split-model', - action="store_true", + '--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() # convert model to python 3 if needed if(args.convert) : prepForConversion(args.ufodir) convertToPython3(args.ufodir) # import the model try : import importlib,importlib.util try : path,mod = os.path.split(os.path.abspath(args.ufodir)) spec = importlib.util.spec_from_file_location(mod,os.path.join(os.path.abspath(args.ufodir),"__init__.py")) FR = spec.loader.load_module() except : print ("Could not load the UFO python module.\n", "This is usually because you are using python3 and the UFO model is in python2.\n", "The --convert option can be used to try and convert it, but there is no guarantee.\n", "As this modifies the model you should make sure you have a backup copy.\n") quit() except : print ("Newer python import did not work, reverting to python 2 imp module") import imp path,mod = os.path.split(os.path.abspath(args.ufodir)) FR = imp.load_module(mod,*imp.find_module(mod,[path])) ################################################## ################################################## # get the Model name from the arguments modelname = args.name libname = modelname + '.so' -# define arrays and variables +# 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' ] ) +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) + import cmath + return eval(x, + {'sec':FR.function_library.sec, + 'csc':FR.function_library.csc, + '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 +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 + + +# more arrays used for substitution in templates paramsforstream = [] parmmodelconstr = [] # loop over parameters and fill in template stuff according to internal/external and complex/real # WARNING: Complex external parameter input not tested! if args.verbose: print ('verbose mode on: printing all parameters') print ('-'*60) paramsstuff = ('name', 'expression', 'default value', 'nature') pprint.pprint(paramsstuff) interfacedecl_T = """\ static Parameter<{modelname}, {type}> interface{pname} ("{pname}", "The interface for parameter {pname}", &{modelname}::{pname}_, {value}, 0, 0, false, false, Interface::nolimits); """ # sort out the couplings couplingDefns = { "QED" : 99, "QCD" : 99 } try : for coupling in FR.all_orders: name = coupling.name.upper() couplingDefns[name]= coupling.expansion_order except: for coupling in FR.all_couplings: for name,value in coupling.order.items(): if(name not in couplingDefns) : couplingDefns[name]=99 # sort out the particles massnames = {} widthnames = {} for particle in FR.all_particles: # skip ghosts and goldstones if(isGhost(particle) or isGoldstone(particle)) : continue if particle.mass != 'ZERO' and particle.mass.name != 'ZERO': if(particle.mass in massnames) : if(abs(particle.pdg_code) not in massnames[particle.mass]) : massnames[particle.mass].append(abs(particle.pdg_code)) else : massnames[particle.mass] = [abs(particle.pdg_code)] if particle.width != 'ZERO' and particle.width.name != 'ZERO': if(particle.width in widthnames) : if(abs(particle.pdg_code) not in widthnames[particle.width]) : widthnames[particle.width].append(abs(particle.pdg_code)) else : widthnames[particle.width] = [abs(particle.pdg_code)] interfaceDecls = [] modelparameters = {} for p in FR.all_parameters: value = parmsubs[p.name] if p.type == 'real': assert( value.imag < 1.0e-16 ) value = value.real if p.nature == 'external': if p not in massnames and p not in widthnames: - interfaceDecls.append( + interfaceDecls.append( interfacedecl_T.format(modelname=modelname, pname=p.name, value=value, - type=typemap(p.type)) + 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('()()','()') + text=text.replace('sec()(','sec(') + text=text.replace('csc()(','csc(') doinit.append(' %s_ = %s;' % (p.name, text) ) if(p.type == 'complex') : doinit.append(' if(std::isnan(%s_.real()) || std::isnan(%s_.imag()) || std::isinf(%s_.real()) || std::isinf(%s_.imag())) {throw InitException() << "Calculated parameter %s is nan or inf in Feynrules model. Check your input parameters.";}' % (p.name,p.name,p.name,p.name,p.name) ) else : doinit.append(' if(std::isnan(%s_) || std::isinf(%s_)) {throw InitException() << "Calculated parameter %s is nan or inf in Feynrules model. Check your input parameters.";}' % (p.name,p.name,p.name) ) if p in massnames: for idCode in massnames[p] : doinit.append(' resetMass(%s,%s_ * GeV);' % (idCode, p.name) ) if p in widthnames: for idCode in widthnames[p] : doinit.append(' getParticleData(%s)->width(%s_ * GeV);' % (idCode, p.name) ) doinit.append(' getParticleData(%s)->cTau (%s_ == 0.0 ? Length() : hbarc/(%s_*GeV));' % (idCode, p.name, p.name) ) doinit.append(' getParticleData(%s)->widthCut(10. * %s_ * GeV);' % (idCode, p.name) ) elif p.nature == 'external': if p in massnames: for idCode in massnames[p] : doinit.append(' %s_ = getParticleData(%s)->mass() / GeV;' % (p.name, idCode) ) if p in widthnames: for idCode in widthnames[p] : doinit.append(' %s_ = getParticleData(%s)->width() / GeV;' % (p.name, idCode) ) if args.verbose: pprint.pprint((p.name,p.value, value, p.nature)) pcwriter = ParamCardWriter(FR.all_parameters) paramcard_output = '\n'.join(pcwriter.output) ### special treatment # if p.name == 'aS': # expression = '0.25 * sqr(strongCoupling(q2)) / Constants::pi' # elif p.name == 'aEWM1': # expression = '4.0 * Constants::pi / sqr(electroMagneticCoupling(q2))' # elif p.name == 'Gf': # expression = 'generator()->standardModel()->fermiConstant() * GeV2' paramconstructor=': ' for ncount in range(0,len(parmconstr)) : paramconstructor += parmconstr[ncount] if(ncount != len(parmconstr) -1) : paramconstructor += ',' if(ncount%5 == 0 ) : paramconstructor += "\n" paramout="" paramin ="" for ncount in range(0,len(paramsforstream)) : if(ncount !=0 ) : paramout += "<< " + paramsforstream[ncount] paramin += ">> " + paramsforstream[ncount] else : paramout += paramsforstream[ncount] paramin += paramsforstream[ncount] if(ncount%5 == 0 ) : paramout += "\n" paramin += "\n" parmtextsubs = { 'parmgetters' : '\n'.join(parmgetters), 'parmdecls' : '\n'.join(parmdecls), 'parmconstr' : paramconstructor, 'getters' : '', 'decls' : '', 'addVertex' : '', 'doinit' : '\n'.join(doinit), 'ostream' : paramout, 'istream' : paramin , 'refs' : '', 'parmextinter': ''.join(interfaceDecls), 'ModelName': modelname, 'calcfunctions': '', 'param_card_data': paramcard_output } ################################################## ################################################## ################################################## # set up the conversion of the vertices vertexConverter = VertexConverter(FR,parmvalues,couplingDefns) vertexConverter.readArgs(args) # convert the vertices vertexConverter.convert() cdefs="" couplingOrders="" ncount=2 for name,value in couplingDefns.items() : if(name=="QED") : couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" %(name,1,value) elif (name=="QCD") : couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" %(name,2,value) else : ncount+=1 cdefs +=" const T %s = %s;\n" % (name,ncount) couplingOrders+=" setCouplings(\"%s\",make_pair(%s,%s));\n" % (name,ncount,value) # coupling definitions couplingTemplate= """\ namespace ThePEG {{ namespace Helicity {{ namespace CouplingType {{ typedef unsigned T; /** * Enums for couplings */ {coup} }} }} }} """ if(cdefs!="") : cdefs = couplingTemplate.format(coup=cdefs) parmtextsubs['couplings'] = cdefs parmtextsubs['couplingOrders'] = couplingOrders # particles plist, names = thepeg_particles(FR,parmsubs,modelname,modelparameters,args.forbidden_particle_name,args.use_Herwig_Higgs) particlelist = [ "# insert HPConstructor:Outgoing 0 /Herwig/{n}/Particles/{p}".format(n=modelname,p=p) for p in names ] # make the first one active to have something runnable in the example .in file particlelist[0] = particlelist[0][2:] particlelist = '\n'.join(particlelist) modelfilesubs = { 'plist' : plist, 'vlist' : vertexConverter.get_vertices(libname), 'setcouplings': '\n'.join(parmmodelconstr), 'ModelName': modelname } # write the files from templates according to the above subs if vertexConverter.should_print(): MODEL_HWIN = getTemplate('LHC-FR.in') if(not args.split_model) : MODEL_CC = [getTemplate('Model.cc')] else : MODEL_EXTRA_CC=getTemplate('Model6.cc') extra_names=[] extra_calcs=[] parmtextsubs['doinit']="" for i in range(0,len(doinit)) : if( i%20 == 0 ) : function_name = "initCalc" +str(int(i/20)) parmtextsubs['doinit'] += function_name +"();\n" parmtextsubs['calcfunctions'] += "void " + function_name + "();\n" extra_names.append(function_name) extra_calcs.append("") extra_calcs[-1] += doinit[i] + "\n" for i in range(0,len(extra_names)) : ccname = '%s_extra_%s.cc' % (modelname,i) writeFile( ccname, MODEL_EXTRA_CC.substitute({'ModelName' : modelname, 'functionName' : extra_names[i], 'functionCalcs' : extra_calcs[i] }) ) MODEL_CC = [getTemplate('Model1.cc'),getTemplate('Model2.cc'),getTemplate('Model3.cc'), getTemplate('Model4.cc'),getTemplate('Model5.cc')] MODEL_H = getTemplate('Model.h') print ('LENGTH',len(MODEL_CC)) MODELINFILE = getTemplate('FR.model') - writeFile( 'LHC-%s.in' % modelname, + writeFile( 'LHC-%s.in' % modelname, MODEL_HWIN.substitute({ 'ModelName' : modelname, 'Particles' : particlelist }) ) modeltemplate = Template( MODELINFILE.substitute(modelfilesubs) ) writeFile( '%s.h' % modelname, MODEL_H.substitute(parmtextsubs) ) for i in range(0,len(MODEL_CC)) : if(len(MODEL_CC)==1) : ccname = '%s.cc' % modelname else : ccname = '%s.cc' % (modelname + str(i)) writeFile( ccname, MODEL_CC[i].substitute(parmtextsubs) ) writeFile( modelname +'.template', modeltemplate.template ) writeFile( modelname +'.model', modeltemplate.substitute( modelparameters ) ) # copy the Makefile-FR to current directory, # replace with the modelname for compilation with open(os.path.join(modulepath,'Makefile-FR'),'r') as orig: with open('Makefile','w') as dest: dest.write(orig.read().replace("FeynrulesModel.so", libname)) print('finished generating model:\t', modelname) print('model directory:\t\t', args.ufodir) print('generated:\t\t\t', len(FR.all_vertices), 'vertices') print('='*60) print('library:\t\t\t', libname) print('input file:\t\t\t', 'LHC-' + modelname +'.in') print('model file:\t\t\t', modelname +'.model') print('='*60) print("""\ To complete the installation, compile by typing "make". -An example input file is provided as LHC-FRModel.in, +An example input file is provided as LHC-FRModel.in, you'll need to change the required particles in there. """) print('DONE!') print('='*60) - diff --git a/Models/Feynrules/python/ufo2peg/Model.cc.template b/Models/Feynrules/python/ufo2peg/Model.cc.template --- a/Models/Feynrules/python/ufo2peg/Model.cc.template +++ b/Models/Feynrules/python/ufo2peg/Model.cc.template @@ -1,94 +1,102 @@ // -*- C++ -*- // // ${ModelName}.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2007 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ${ModelName} class. // #include "${ModelName}.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/Vertex/VertexBase.h" #include "ThePEG/Utilities/DescribeClass.h" #include //#include "Herwig/Models/General/ModelGenerator.h" using namespace ThePEG; using namespace Herwig; ${ModelName}::${ModelName}() ${parmconstr} {} IBPtr ${ModelName}::clone() const { return new_ptr(*this); } IBPtr ${ModelName}::fullclone() const { return new_ptr(*this); } +double sec(double x) { + return 1./cos(x); +} + +double csc(double x) { + return 1./sin(x); +} + void ${ModelName}::doinit() { ${doinit} BSMModel::doinit(); ${addVertex} writeParamCard(); ${couplingOrders} } void ${ModelName}::doinitrun() { BSMModel::doinitrun(); writeParamCard(); } void ${ModelName}::persistentOutput(PersistentOStream & os) const { os << ${ostream} ; } void ${ModelName}::persistentInput(PersistentIStream & is, int) { is >> ${istream} ; } void ${ModelName}::writeParamCard() const { ofstream card("param_card.dat"); card << "#####################################################\n" << "## DO NOT EDIT - GENERATED BY HERWIG UFO CONVERTER ##\n" << "#####################################################\n\n"; card ${param_card_data} << '\n'; card.close(); } // Static variable needed for the type description system in ThePEG. DescribeClass<${ModelName},BSMModel> describeThePEG${ModelName}("Herwig::${ModelName}", "${ModelName}.so"); void ${ModelName}::Init() { ${refs} static ClassDocumentation<${ModelName}> documentation ("The ${ModelName} class inherits from BSMModel" "and supplies additional couplings and access to the ${ModelName}" "vertices for helicity amplitude calculations" ); ${parmextinter} } diff --git a/Models/Feynrules/python/ufo2peg/Model1.cc.template b/Models/Feynrules/python/ufo2peg/Model1.cc.template --- a/Models/Feynrules/python/ufo2peg/Model1.cc.template +++ b/Models/Feynrules/python/ufo2peg/Model1.cc.template @@ -1,38 +1,47 @@ // -*- C++ -*- // // ${ModelName}.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2007 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ${ModelName} class. // #include "${ModelName}.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/Vertex/VertexBase.h" #include "ThePEG/Utilities/DescribeClass.h" //#include "Herwig/Models/General/ModelGenerator.h" using namespace ThePEG; using namespace Herwig; ${ModelName}::${ModelName}() ${parmconstr} {} IBPtr ${ModelName}::clone() const { return new_ptr(*this); } IBPtr ${ModelName}::fullclone() const { return new_ptr(*this); } + + +double sec(double x) { + return 1./cos(x); +} + +double csc(double x) { + return 1./sin(x); +} 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,3030 +1,3028 @@ 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) raise SkipThisVertex() else : 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) + 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) 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) 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) 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) 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) 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) 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) 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) 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") : + 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) : + 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) : + 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]) 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 : + 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]) 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") 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]) 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') 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) 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) 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) 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') 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) 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]) 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) 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 + ( "%s*%s*%s" %(start[ichain],etemp,end[ichain])), 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 + ( "%s*%s*Transpose(%s)*%s*%s" %(startT[ichain],CC,etemp,CD,endT[ichain])), 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) 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) 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) 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) 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 += "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) 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.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) raise SkipThisVertex() # unknown to be sorted out else : 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]) 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]) 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]) 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]) 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) 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') raise SkipThisVertex() for ichain in range(0,nchain) : eTemp[ichain].append("(%s)"% (value) ) else : print('Uncontracted on-shell tensor') raise SkipThisVertex() # unknown else : 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.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])) raise SkipThisVertex() break # uncontracted indices else : istring = "" for (key,val) in unContracted.items() : istring +=imap[val] if(len(istring)>2) : 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.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)) 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)) 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 + nchain = int((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) 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 : + else : 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].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) for j in range(0,len(vals)) : print(j,dim[j],vals[j]) raise SkipThisVertex() - # simplest case - if(isinstance(vals[i], basestring)) : + # simplest case + if(isinstance(vals[i], str)) : for ii in range(0,len(expr)) : 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') 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]) 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].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].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)) : + if(not isinstance(result[0],str)) : subs=[] for ii in range(0,len(result)) : subs.append({}) 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 + 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) 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].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].items(): + 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") 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]) raise SkipThisVertex() # check if momenta defns needed to clean up compile of code for (key,val) in defns.items() : - if( isinstance(key, basestring)) : + if( isinstance(key, str)) : 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.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.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) 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)) - -