Define spin up particle spinor:

(uup[px_, py_, pz_, e_, 
            m_] = (Sqrt[e + m]) {1, 0, 
              (pz/(e + m)), ((px + I py)/(e + m))}) // 
      MatrixForm

Spin down particle spinor:

(udown[px_, py_, pz_, e_, 
            m_] = (Sqrt[e + m]) {0, 
              1, ((px - I py)/(e + m)), -pz/(e + m)}) // 
      MatrixForm

Spin down antiparticle spinor:

(vdown[px_, py_, pz_, e_, 
            m_] = (Sqrt[e + m]) {((px - I py)/(e + m)), -pz/(e + 
                    m), 0, 1}) // MatrixForm

Spin up antiparticle spinor:

(vup[px_, py_, pz_, e_, 
            m_] = -(Sqrt[e + m]) {pz/(e + 
                    m), ((px + I py)/(e + m)), 1, 0}) // 
      MatrixForm

Define gamma matrices:

MatrixForm[
      gamma0 = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 
            0, -1}}]

MatrixForm[
      gamma1 = {{0, 0, 0, 1}, {0, 0, 1, 0}, {0, -1, 0, 0}, {-1, 0, 0, 
            0}}]

MatrixForm[
      gamma2 = {{0, 0, 0, -I}, {0, 0, I, 0}, {0, I, 0, 0}, {-I, 0, 0, 
            0}}]

MatrixForm[
      gamma3 = {{0, 0, 1, 0}, {0, 0, 0, -1}, {-1, 0, 0, 0}, {0, 1, 0, 
            0}}]

Define Current:

current[a_, 
        b_] := {(a /. Complex[x_, y_] -> Complex[x, -y]) . 
          gamma0 . gamma0 . 
          b, (a /. Complex[x_, y_] -> Complex[x, -y]) . gamma0 . 
          gamma1 . b, (a /. Complex[x_, y_] -> Complex[x, -y]) . 
          gamma0 . gamma2 . 
          b, (a /. Complex[x_, y_] -> Complex[x, -y]) . gamma0 . 
          gamma3 . b}

Vector of the two electron spin states:

ee  =  {uup[0, 0, p, e, me], udown[0, 0, p, e, me]}

Vector of the two positron spin states:

pp = {vup[0, 0, -p, e, me], vdown[0, 0, -p, e, me]}

Vector of the four possible spin combinations:

matrix = 
      Simplify[Flatten[
            Table[current[pp[[i]], ee[[j]]], {i, 1, 2}, {j, 1, 2}], 
            1] /. p -> Sqrt[e^2 - me^2]]

Vector of the muon minus two spin states:

mum  =  {uup[pm Sin[theta], 0, pm Cos[theta], e, mmu], 
        udown[pm Sin[theta], 0, pm Cos[theta], e, mmu]}

Vector of the muon plus two spin states:

mup  =  {vup[-pm Sin[theta], 0, -pm Cos[theta], e, mmu], 
        vdown[-pm Sin[theta], 0, -pm Cos[theta], e, 
          mmu]}

Vector of the four possible muon spin combinations:

mmatrix = 
      Simplify[Flatten[
            Table[current[mum[[i]], mup[[j]]], {i, 1, 2}, {j, 1, 2}], 
            1] /. pm -> Sqrt[e^2 - mmu^2]]

currentcurrent = 
      Flatten[Simplify[
          Table[matrix[[i]] . mmatrix[[k]], {i, 1, 4}, {k, 1, 
              4}]]]

msqp = 
      Table[Simplify[
          Simplify[
              currentcurrent[[i]]*currentcurrent[[i]]/(2 e)^4] /. 
            Cos[2 theta] -> 2 Cos[theta]^2 - 1, 
          Trig -> False], {i, 1, 16}]

conditions = {e ->  .25, me ->  .0005, 
        mmu ->  .105}

Table[
      Plot[msqp[[i]] /. conditions, {theta, 0, Pi}, 
        PlotRange -> {0, 4}], {i, 1, 16}]

msq = 
      Simplify[Simplify[Apply[Plus, msqp]/4] /. 
          Cos[2 theta] -> 2 Cos[theta]^2 - 1, 
        Trig -> False]

aa = 
      Plot[msq /. conditions, {theta, 0, Pi}, 
        PlotRange -> {0, 4}]

bb = 
      Plot[1 + Cos[theta]^2, {theta, 0, Pi}, 
        PlotRange -> {0, 2}]

Show[aa, bb]