(* Program to Calculate the Helicity Amplitudes in Bhabha scattering,
   Eq. (6.17) of QFT 

 Compatible with FC 9.3 and FC 10.1.0

 Last Version: 07/12/2040
 Author: Jorge C. Romao
 email: jorge.romao@tecnico.ulisboa.pt
*)

Remove["Global`*"]

(* Definitions *)
SetOptions[DiracSimplify,DiracSubstitute67->False]


dp=Function[s,If[s==1,GA[6],If[s==-1,GA[7]]]]
	    
U[p_,s_]:= dp[s] . SpinorU[p,0]

UBar[p_,s_]:=  SpinorUBar[p,0] . dp[-s]

gvga= gm dp[-1] + gp dp[1]

delta[s1_,s2_]:=If[s1==s2,1,0]

(* 02.05.2020
 
To be compatible with FeynCalc 9.3 we have to separate the spinor and 
spinorBar in the definitions of GS[p] and Polarization sums.

The code below works for FeynCalc 9.3 and 10.1.0
*)


PolSConjU=Function[{k,p,l,i}, If[i==1, 2 U[k,l],If[i==2, 2 U[p,-l]]]]

PolSConjUBar=Function[{k,p,l,i}, If[i==1,UBar[p,l],If[i==2, UBar[k,-l]]]]

MydsU=Function[{p,i},If[i==1, U[p,1],If[i==2, U[p,-1]]]] 

MydsUBar=Function[{p,i},If[i==1, UBar[p,1],If[i==2, UBar[p,-1]]]] 

(* Line Ubar(p1,s1).GA[mu].U(p2,s2) Ubar[p3,s3].GA[mu].U[p4,s4] *)  
  
UBarGammaU = Function[{p1, s1, p2, s2, p3, s3, p4, s4},
X1 = delta[s1,s2] DiracSimplify[2 UBar[p3, s3].U[p1, -s1] UBar[p2, -s1].U[p4, s4]];
X2 = delta[s1,s2] DiracSimplify[2 UBar[p3, s3].U[p2, s2] UBar[p1, s1].U[p4, s4]];
X1 + X2]
  
(* Amplitudes *)

  
M1[s1_,s2_,s3_,s4_]:=UBarGammaU[p2,s2,p1,s1,p3,s3,p4,s4]

res1[s1_,s2_,s3_,s4_]:=DiracSimplify[M1[s1,s2,s3,s4]]/s  

M2[s1_,s2_,s3_,s4_]:=UBarGammaU[p3,s3,p1,s1,p2,s2,p4,s4] 

res2[s1_,s2_,s3_,s4_]:=-DiracSimplify[M2[s1,s2,s3,s4]]/t 


(* Simplify *)
vlist={p1,p2,p3,p4}

simp1=Table[SpinorU[vlist[[i]],0] . GA[7]. SpinorU[vlist[[j]],0] -> sp[vlist[[i]],vlist[[j]]],{i,1,4},{j,1,4}] /. {sp[p_, p_] -> 0}

simp2=Table[SpinorU[vlist[[i]],0] . DiracMatrix[6] . SpinorU[vlist[[j]],0] -> spc[vlist[[j]],vlist[[i]]],{i,1,4},{j,1,4}] /. {spc[q_, q_] -> 0}

simp=Flatten[{simp1,simp2}]//FCI;


M[s1_,s2_,s3_,s4_]:=Expand[res1[s1,s2,s3,s4]+res2[s1,s2,s3,s4] /. simp]


