(* Program to Calculate the Traces of the bremstrahlung and compare
   with the Bethe-Heitler formula. See ITC Eqs 5.97, and 5.99, 
   pag 186.

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

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

  
GammaM[mu_]:=-1/2/w/Di  GS[z] . (GS[pi]-GS[k] +m). GA[mu] + \
  1/2/w/Df GA[mu] . ( GS[pf] + GS[k] +m) . GS[z]
GammaMBar[mu_]:= -1/2/w/Di  GA[mu] . (GS[pi]-GS[k] +m). GS[z] + \
  1/2/w/Df GS[z] . ( GS[pf] + GS[k] +m) . GA[mu]


Line1:= - prop[pf,m] . GammaM[mu] . prop[pi,m] . GammaMBar[mu]
  ans1= Simplify[Contract[TR[Line1]]]//FCE


kin={SP[z,z]->1,SP[z,pi]->Ei,SP[z,pf]->Ef,SP[z,k]->w,SP[pi,pi]->m^2,SP[pf,pf]->m^2,SP[k,pi]->w Di ,SP[k,pf]->w Df,SP[k,k]->0,SP[pi,pf]-> Ei Ef + 1/2 (q2 -Ei^2 -Ef^2 +2 m^2 -w^2 +2 w (Ei-Di) -2 w (Ef-Df))}

kin2={SP[z,z]->1,SP[z,pi]->Ei,SP[z,pf]->Ef,SP[z,k]->w,SP[pi,pi]->m^2,SP[pf,pf]->m^2,SP[k,pi]->w Di ,SP[k,pf]->w Df,SP[k,k]->0}

Eq597=ans1 /. kin2//Expand

Print["Eq 5.97 =", Eq597]

ans2= ans1 /. kin//Simplify
MyBH=Expand[ans2]


(* Bethe-Heitler *)
BH1=pf^2 sf^2 (4 Ei^2 -q2)/ Df^2 
BH2=pi^2 si^2 (4 Ef^2 -q2)/ Di^2 
BH3=2 w^2/Di/Df (pi^2 si^2 + pf^2 sf^2) 
BH4=-2 pi pf /Di/Df si sf cp (2 Ei^2 + 2 Ef^2 -q2) 

BH=2/w^2*(BH1+BH2+BH3+BH4) //Simplify


(* Compare with Bethe-Heitler *)
sub1={pi^2 si^2 -> Ei^2 - m^2 - (Ei - Di)^2}
sub2={pf^2 sf^2 -> Ef^2 - m^2 - (Ef - Df)^2}
sub3 ={pi pf si sf cp-> -1/2 (q2 - Ei^2 - Ef^2 + 2 m^2 - w^2 + 2 w (Ei - Di) -  2 w (Ef - Df))-(Ei-Di) (Ef-Df)};
sub4={Ef->Ei-w}

TestBH:= Simplify[BH-MyBH /.sub1 /. sub2 /. sub3 /. sub4]

Print["\n Eq. (5.99) - Eq. (5.97) =", TestBH]
