(* Program to Calculate the Helicity Amplitudes Z-> f fbar decay,
   Eq. (6.39) of ITC 2020 

Works for Feyncalc 9.3 and 10.1.0
 
 Last Version: 08/12/2024
 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]


UBarGammaUGammaU=Function[{p1,p2,si,i},If[i==1, 2 U[p2,si],If[i==2, 2 U[p1,-si]]]];

UBarGammaUGammaUBar=Function[{p1,p2,si,i},If[i==1, UBar[p1,si],If[i==2, UBar[p2,-si]]]];

  
(* Amplitudes *)

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

M1[s1_,s2_]:= -gZ Sum[(UBar[q1,s1] . UBarGammaUGammaU[r1,r2,-1,i1] UBarGammaUGammaUBar[r1,r2,-1,i1] . gvga . U[q2,s2]),{i1,1,2}]

res1[s1_,s2_]:=DiracSimplify[M1[s1,s2]]

vlist={q1,q2,r1,r2}

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] . GA[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_]:=Expand[(res1[s1,s2] /. simp)]


