(* Program to Calculate the Helicity Amplitudes in Bhabha scattering,
   Eq. (5.99) of ITC 2020 

Compatible with FC 9.2 and 9.3

 Last Version: 01/05/2020
 Author: Jorge C. Romao
 email: jorge.romao@tecnico.ulisboa.pt
*)

Remove["Global`*"]

(* Definitions *)

dp[s_]:= (1 + s DiracMatrix[5])/2
U[p_,s_]:= dp[s] . Spinor[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.2 and 9.3
*)


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]]]] 
  
UBarGammaU = Function[{p3, s3, p2, p1, s1, p4, s4},
    X1 = DiracSimplify[2 UBar[p3, s3].U[p1, s1] UBar[p2, s1].U[p4, s4]];
    X2 = DiracSimplify[2 UBar[p3, s3].U[p2, -s1] UBar[p1, -s1].U[p4, s4]];
    X1 + X2]
  
(* Amplitudes *)

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

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

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

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


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

simp1=Table[Spinor[vlist[[i]],0] . Spinor[vlist[[j]],0] -> sp[vlist[[i]],vlist[[j]]] + spc[vlist[[j]],vlist[[i]]],{i,1,4},{j,1,4}] /. {sp[p_, p_] -> 0, spc[q_, q_] -> 0}

simp2=Table[Spinor[vlist[[i]],0] . DiracMatrix[5] . Spinor[vlist[[j]],0] -> -sp[vlist[[i]],vlist[[j]]] + spc[vlist[[j]],vlist[[i]]],{i,1,4},{j,1,4}] /. {sp[p_, p_] -> 0, spc[q_, q_] -> 0}


simp=Flatten[{simp1,simp2}];


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


