-- www.eigenmath.org/muon-decay-2.txt -- Verify dot product formula for probability density. p1 = (E1, p1x, p1y, p1z) -- muon p2 = (E2, p2x, p2y, p2z) -- muon neutrino p3 = (E3, p3x, p3y, p3z) -- electron antineutrino p4 = (E4, p4x, p4y, p4z) -- electron I = ((1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)) gmunu = ((1,0,0,0),(0,-1,0,0),(0,0,-1,0),(0,0,0,-1)) gamma0 = ((1,0,0,0),(0,1,0,0),(0,0,-1,0),(0,0,0,-1)) gamma1 = ((0,0,0,1),(0,0,1,0),(0,-1,0,0),(-1,0,0,0)) gamma2 = ((0,0,0,-i),(0,0,i,0),(0,i,0,0),(-i,0,0,0)) gamma3 = ((0,0,1,0),(0,0,0,-1),(-1,0,0,0),(0,1,0,0)) gamma5 = i dot(gamma0,gamma1,gamma2,gamma3) gamma = (gamma0,gamma1,gamma2,gamma3) gammaT = transpose(gamma) gammaL = transpose(dot(gmunu,gamma)) pslash1 = dot(p1,gmunu,gamma) pslash2 = dot(p2,gmunu,gamma) pslash3 = dot(p3,gmunu,gamma) pslash4 = dot(p4,gmunu,gamma) -- T1 is the first trace matrix T1 = contract(dot(pslash4,gammaT,I - gamma5,pslash3,gammaT,I - gamma5),1,4) -- T2 is the second trace matrix T2 = contract(dot(pslash2,gammaL,I - gamma5,pslash1,gammaL,I - gamma5),1,4) -- T is the product of T1 and T2 T = contract(dot(T1,transpose(T2))) -- print 1 if equal T == 256 dot(p1,gmunu,p3) dot(p2,gmunu,p4)