SWIG/6.py at main · HepLib/SWIG · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#!/usr/bin/env python3

# python version for 6.cpp

from HepLib import *

A = Symbol("A")
e = Symbol("e")
ebar = Symbol("ebar")
mu = Symbol("mu")
mubar = Symbol("mubar")

p = Vector("p")
P = Vector("P")
k = Vector("k")
K = Vector("K")

me = Symbol("me")
mm = Symbol("mm")

proc = Process()
proc.Model = """
    [e, ebar, -]
    [mu, mubar, -]
    [A, A, +]
    [ebar, e, A]
    [mubar, mu, A]
"""
proc.In = "e[p],ebar[P]"
proc.Out = "mubar[k],mu[K]"
proc.Options = "onshell"
proc.Loops = 0

st = {}
st["p"] = p
st["P"] = P
st["k"] = k
st["K"] = K

amps = proc.Amplitudes(st)

set_InOutTeX(-1,"$e^-(p)$")
set_InOutTeX(-3,"$e^+(P)$")
set_InOutTeX(-2,"$\\mu^+(k)$")
set_InOutTeX(-4,"$\\mu^-(K)$")

set_LineTeX(Symbol("e"),"fermion, edge label=$e$")
set_LineTeX(Symbol("ebar"),"anti fermion, edge label=$e$")
set_LineTeX(Symbol("mu"),"fermion, edge label=$\\mu$")
set_LineTeX(Symbol("mbar"),"anti fermion, edge label=$\\mu$")
set_LineTeX(Symbol("A"),"photon, edge label=$\\gamma$")

Process.DrawPDF(amps, "amps.pdf")

class ClassFR(MapFunction):
    def __init__(self):
        MapFunction.__init__(self)

    def map(self, e):
        if(isFunction(e,"OutField") or isFunction(e,"InField")):
            return expr(1)
        elif(isFunction(e, "Propagator")):
            fi1 = e.op(0).op(1)
            fi2 = e.op(1).op(1)
            mom = e.op(2)
            if(e.op(0).op(0)==A):
                return (-I) * SP(LI(fi1),LI(fi2)) / SP(mom); # Feynman Gauge
            elif(e.op(0).op(0)==ebar):
                return I * Matrix(GAS(mom)+GAS(1)*me, DI(fi1),DI(fi2)) / (SP(mom)-me*me)
            elif(e.op(0).op(0)==mubar):
                return I * Matrix(GAS(mom)+GAS(1)*mm, DI(fi1),DI(fi2)) / (SP(mom)-mm*mm)
        elif(isFunction(e, "Vertex")):
            fi1 = e.op(0).op(1)
            fi2 = e.op(1).op(1)
            fi3 = e.op(2).op(1)
            if(e.op(0).op(0)==ebar):
                return I*Symbol("e")*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2))
            elif(e.op(0).op(0)==mubar):
                return I*Symbol("e")*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2))
        else:
            return e.map(self)

amps_FR = ClassFR()(amps[0])

co << "amps_FR: " << endl
co << amps_FR << endl << endl

ampL = amps_FR
ampR = IndexL2R(conjugate(ampL));
def SS1(p,m,i):
    return Matrix(GAS(p)+m*GAS(1),DI(i),RDI(i))
def SS2(p,m,i):
    return Matrix(GAS(p)+m*GAS(1),RDI(i),DI(i))

M2 = ampL * ampR * SS1(p,me,-1) * SS2(P,-me,-3) * SS1(k,-mm,-2) * SS2(K,mm,-4);
M2 = MatrixContract(M2);

co << "M2: " << endl
co << M2 << endl << endl

set_form_using_dim4(True)
letSP(p,me*me)
letSP(P,me*me)
letSP(k,mm*mm)
letSP(K,mm*mm)

res = form(M2);
co << "Final M2:" << endl
co << factor(res.subs(me>>0)) << endl