-- Verify formulas for muon pair production p = sqrt(E^2 - m^2) rho = sqrt(E^2 - M^2) p1 = (E, 0, 0, p) p2 = (E, 0, 0, -p) p3 = (E, rho expsin(theta) expcos(phi), rho expsin(theta) expsin(phi), rho expcos(theta)) p4 = (E, -rho expsin(theta) expcos(phi), -rho expsin(theta) expsin(phi), -rho expcos(theta)) 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)) 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) X1 = pslash1 + m I X2 = pslash2 - m I X3 = pslash3 + M I X4 = pslash4 - M I T1 = contract(dot(X3,gammaT,X4,gammaT),1,4) T2 = contract(dot(X2,gammaL,X1,gammaL),1,4) T = contract(dot(T1,transpose(T2))) "Verify probability density (1=ok)" f = 1 + expcos(theta)^2 + (m^2 + M^2) E^(-2) expsin(theta)^2 + m^2 M^2 E^(-4) expcos(theta)^2 T == 64 E^4 f "Verify momentum formula (1=ok)" f = 32 dot(p1,gmunu,p3) dot(p2,gmunu,p4) + 32 dot(p1,gmunu,p4) dot(p2,gmunu,p3) + 32 m^2 dot(p3,gmunu,p4) + 32 M^2 dot(p1,gmunu,p2) + 64 M^2 m^2 T == f "Verify Mandelstam formula (1=ok)" s = dot(p1 + p2,gmunu,p1 + p2) t = dot(p1 - p3,gmunu,p1 - p3) u = dot(p1 - p4,gmunu,p1 - p4) f = 16 t^2 - 16 s u + 32 u (m^2 + M^2) - 8 s^2 + 64 M^2 s + 64 m^2 s - 48 (m^2 + M^2)^2 T == f