-- 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 sin(theta) cos(phi),
      rho sin(theta) sin(phi),
      rho cos(theta))

p4 = (E,
      -rho sin(theta) cos(phi),
      -rho sin(theta) sin(phi),
      -rho cos(theta))

u11 = (E + m, 0, p1[4], p1[2] + i p1[3]) / sqrt(E + m)
u12 = (0, E + m, p1[2] - i p1[3], -p1[4]) / sqrt(E + m)

v21 = (p2[4], p2[2] + i p2[3], E + m, 0) / sqrt(E + m)
v22 = (p2[2] - i p2[3], -p2[4], 0, E + m) / sqrt(E + m)

u31 = (E + M, 0, p3[4], p3[2] + i p3[3]) / sqrt(E + M)
u32 = (0, E + M, p3[2] - i p3[3], -p3[4]) / sqrt(E + M)

v41 = (p4[4], p4[2] + i p4[3], E + M, 0) / sqrt(E + M)
v42 = (p4[2] - i p4[3], -p4[4], 0, E + M) / sqrt(E + M)

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))

"Verify Casimir trick"

u1 = (u11,u12)
v2 = (v21,v22)
u3 = (u31,u32)
v4 = (v41,v42)

v2bar = dot(conj(v2),gamma0) -- adjoint of v2
u3bar = dot(conj(u3),gamma0) -- adjoint of u3

M1(a,b,c,d) = dot(dot(u3bar[c],gammaL,v4[d]),dot(v2bar[b],gammaT,u1[a]))

MM = sum(a,1,2,sum(b,1,2,sum(c,1,2,sum(d,1,2, M1(a,b,c,d) conj(M1(a,b,c,d)) ))))

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)
f = contract(dot(T1,transpose(T2)))

check(f == MM)

"ok"

"Verify probability density"

s = dot(p1 + p2, gmunu, p1 + p2)

check(4 s^2 == 64 E^4)

check(f / (64 E^4) == 1 + cos(theta)^2 + (m^2 + M^2) / E^2 sin(theta)^2 + m^2 M^2 cos(theta)^2 / E^4)

-- verify integral

f = 1 + cos(theta)^2
I = -1/3 cos(theta)^3 - cos(theta)
check(f sin(theta) == d(I,theta))

-- verify cdf

F = (I - eval(I,theta,0)) / (eval(I,theta,pi) - eval(I,theta,0))
check(F == -1/8 cos(theta)^3 - 3/8 cos(theta) + 1/2)

"ok"