………
defineMatrices 'G_a', 'G5', 'D[m, p_a]', Matrix1.matrix,
'b[p_a]', 'u[p_a]', 'cd[p_a]', Matrix1.vector,
'c[p_a]', 'cu[p_a]', 'd[p_a]', Matrix1.covector
//trace of Dirac matrices
def dTrace = DiracTrace['G_m', 'G5', 'e_abcd']
//simplification of Levi-Civita
def eSimplify = LeviCivitaSimplify.minkowski['e_abcd'.t]
//spinor propagator
def propagator = 'D[m,q_a] = (m + q_a*G^a)/(m**2 - q_a*q^a)'.t
//spin projection for Bc
def spinProjection = 'b[p1_m]*c[p2_m] = F*(p_m*G^m + M)*G5'.t
//masses and other scala combinations
def masses = 'p_m*p^m = M**2'.t & 'p1_m*p1^m = mb**2'.t &
'p2_m*p2^m = mc**2'.t & 'k1_m*k1^m = 0'.t & 'k2_m*k2^m = mu**2'.t & 'k3_m*k3^m = md**2'.t
def scalars = 'k1_a*k2^a = k12'.t & 'k1_a*k3^a = k13'.t & 'k2_a*k3^a = k23'.t
scalars = masses & scalars
//matrix element (ud - denotes ud leg)
//a) - photon emission from b-quark
def Ma = '(-1/3)*c[p2_a]*G_m*(1-G5)*D[mb, p1_a - k1_a]*G_n*b[p1_a] * 1/mW**2 * d[k3_a]*G^m*(1-G5)*u[k2_a]'.t
//b) - photon emission from c-antiquark
def Mb = '(2/3)*c[p2_a]*G_n*D[mc, -p2_a + k1_a]*G_m*(1-G5)*b[p1_a] * 1/mW**2 * d[k3_a]*G^m*(1-G5)*u[k2_a]'.t
//c) - photon emission from d-quark
def Mf = '(-1/3)*c[p2_a]*G_m*(1-G5)*b[p1_a] * 1/mW**2 * d[k3_a]*G_n*D[md, k1_a + k3_a]*G^m*(1-G5)*u[k2_a]'.t
//d) - photon emission from u-antiquark
def Me = '(2/3)*c[p2_a]*G_m*(1-G5)*b[p1_a] * 1/mW**2 * d[k3_a]*G^m*(1-G5)*D[mu, -k1_a - k2_a]*G_n*u[k2_a]'.t
//total matrix element
def M = Ma + Mb + Mf + Me
//substitutions
M = (spinProjection & propagator) >> M
M = 'p1_a = mb/M*p_a' >> M
M = 'p2_a = mc/M*p_a' >> M
//simplification
M = (ExpandAll & EliminateMetrics & scalars) >> M
M = dTrace >> M
M = (ExpandAll & EliminateMetrics & eSimplify & scalars) >> M
//complex conjugation of matrix element
def Mc = (Conjugate & Reverse[Matrix1] & 'u[k2_a]*d[k3_a] = cd[k3_a]*cu[k2_a]'.t & 'G5 = -G5'.t) >> M
//squared matrix element summed over photon polarizations
def M2 = M * ('{_n -> _m}'.mapping >> Mc) * 'g^mn'.t
M2 = ExpandAndEliminate >> M2
//u and d quarks spin sums
def uq = 'u[k2_a]*cu[k2_a] = -mu + k2_j*G^j'.t
def dq = 'cd[k3_a]*d[k3_a] = md + k3_j*G^j'.t
M2 = (uq & dq) >> M2
M2 = dTrace >> M2
//final simplifications
M2 = (ExpandAndEliminate & scalars & 'd_i^i = 4'.t & eSimplify & scalars) >> M2
//momentum conservation
M2 = 'p_a = k1_a + k2_a + k3_a'.t >> M2
M2 = (ExpandAndEliminate & scalars & 'd_i^i = 4'.t & eSimplify & scalars) >> M2
println 'done'
M2 = (ExpandAll & scalars) >> M2
println M2.toString(OutputFormat.WolframMathematica)