(* Program to Calculate the Traces of e E -> G G without neglecting
   the electron mass, as in Problem 5.13 of QFT

 Last Version: 06/12/2024
 Compatible with FeynCalc 10.1.0
 Author: Jorge C. Romao
 email: jorge.romao@tecnico.ulisboa.pt
*)


(* Definitions *)

prop[p_,m_]:=GS[p]+m


(* Traces *)
Line1:= prop[p2,-m] . GA[mu] .  prop[p1-k1,m] . GA[nu] . prop[p1,m] . GA[nu] . 
prop[p1-k1,m] . GA[mu]

Line2:= prop[p2,-m] . GA[mu] .  prop[p1-k2,m] . GA[nu] . prop[p1,m] . GA[nu] . 
prop[p1-k2,m] . GA[mu]

Line12:= prop[p2,-m] . GA[nu] .  prop[p1-k1,m] . GA[mu] . prop[p1,m] . GA[nu] .  prop[p1-k2,m] . GA[mu]

Line21:= prop[p2,-m] . GA[mu] .  prop[p1-k2,m] . GA[nu] . prop[p1,m] . GA[mu] .  prop[p1-k1,m] . GA[nu]


ans1= Simplify[Contract[TR[Line1]]/(t-m^2)^2]
ans2= Simplify[Contract[TR[Line2]]/(u-m^2)^2]
ans12= Simplify[Contract[TR[Line12]]/((t-m^2) (u-m^2))]
ans21= Simplify[Contract[TR[Line21]]/((t-m^2) (u-m^2))]


(* Results for the |Amplitude|^2  *)
  ans=  (ans1 + ans2 + ans12 + ans21) //FCE

onshell={SP[p1,p1]->m^2,SP[p2,p2]->m^2,SP[k1,k1]->0,SP[k2,k2]->0}

res = ans/4 /. onshell  

kin={SP[p1,p2]->(s-2 m^2)/2 ,SP[p1,k1]->(m^2-t)/2 ,SP[p1,k2]->(m^2-u)/2,SP[p2,k1]->(m^2-u)/2,SP[p2,k2]->(m^2-t)/2,SP[k1,k2]->s/2}

aux=Expand[res /. kin]
kin2={s->2 m^2 -u -t} 
aux= Expand[aux /. kin2]

(* Comparing with Eq. (4.50) *)
Eq5271=aux /. {u -> m^2 - U, t -> m^2 - T} // Simplify

Print["\!\(\*UnderscriptBox[\(\[Sum]\), \(spins\)]\)|M\!\(\*SuperscriptBox[\(|\), \(2\)]\) = ",Eq5271]

Print["where:    U->\*SuperscriptBox[m, 2] - u, T-> \*SuperscriptBox[m, 2] - t \n"] 

(* Total Cross Section *)
(* Change to the variable theta and do the integration *)

final = Expand[aux /. {t-> Dt+m^2,u->Du+m^2}]
final2 = final /. {Du -> -s/2 (1 + beta x), Dt -> -s/2 (1 - beta x)}

xs = Integrate[final2, {x, -1, 1}, Assumptions -> {m > 0, s > 0, 
               beta > 0,beta<=1}];

(* 1/2 for identical particles *)
xs = xs/32/Pi/s/(s-4 m^2)^(1/2)*(s)^(1/2) *(4 Pi \[Alpha])^2 *(1/2) 
xs = xs /. beta-> (s-4 m^2)^(1/2)/s^(1/2) // Expand // Simplify

Print["Cross Section = ",xs]

(********************* Fim do Programa **********************)
