Muon decay matrix element

Needs["Tracer`", 
      ToString[ToFileName[{$TopDirectory, "Addons", "Applications", "\
dirac"}, "Tracer.m"]]]

VectorDimension[4]

Spur[1]

muontensor = 
      G[le, pmu + mmu U] ** G[le, {alpha}] ** (G[le, U] - G[le, G5]) ** 
        G[le, pnumu] ** G[le, {mu}] ** (G[le, U] - G[le, G5])

mutensor = muontensor /. le -> 1

electrontensor = 
      G[le, pe + me U] ** G[le, {mu}] ** (G[le, U] - G[le, G5]) ** 
        G[le, pnue] ** 
        G[le, {alpha}] ** (G[le, U] - G[le, G5])

etensor = Simplify[electrontensor /. le -> 1]

msq = 
      Simplify[(1/2) (Gf^2) ContractEpsGamma[etensor mutensor]/2]

Assume electron is massless:

metric = {{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 
          0, -1}}
ppmu = {mmu, 0, 0, 0}
ppe  =  {p, sint p ,  0,  cost p}
ppnue = {k, 0, 0, k}
ppnumu  =  ppmu-ppe-ppnue

msq1 = 
      Simplify[msq /. {pe . pnumu -> ppe . metric . ppnumu, 
              pmu . pnue -> ppmu . metric . ppnue}] /. 
        sint^2 -> 1 - cost^2

Cost is determined by mass shell constraint of neutrino:

eq1 = 
      Simplify[ 
          ppnumu . metric . ppnumu /. sint^2 -> 1 - cost^2]  ==  
        0

sol1  =  Solve[eq1,  cost][[1]]

msq2 = Simplify[msq1 /. sol1]

phasespace  =  1/(8 mmu) 1/(2 Pi)^3

dlambda  =  msq2 * phasespace

With the assumption that the electron is massless, the available area of the Daltiz plot is just the triangle (bordered by k=mmu/2, p=mmu/2, and p=mmu/2-k) within the energy conservation triangle (p=0,k=0, p=mmu-k).  

Integrate with both dp and dk.

dl1  =  Integrate[dlambda, {p, mmu/2 - k, mmu/2}]

lambda = Integrate[dl1, {k, 0, mmu/2}]