====== B_c \to u \bar d \gamma ====== ---- ====Theory==== ......... ====Code==== 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)