(* Program to Calculate the Helicity Amplitudes Z-> f fbar decay 
   with massive fermions, Eq. (6.45) of ITC 

 Compatible with 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]]]];


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

(* Fermions with masss *)

Um[p1_,p2_,m_,s_]:=(If[s==1,sp[p1,p2],spc[p2,p1]]/m U[p1,s] + U[p2,-s])
Vm[p1_,p2_,m_,s_]:=(If[s==1,sp[p1,p2],spc[p2,p1]]/m U[p1,s] - U[p2,-s])

UmBar[p1_,p2_,m_,s_]:=(If[s==1,spc[p1,p2],sp[p2,p1]]/m UBar[p1,s] 
		       + UBar[p2,-s])
VmBar[p1_,p2_,m_,s_]:=(If[s==1,spc[p1,p2],sp[p2,p1]]/m UBar[p1,s] 
		       - UBar[p2,-s])

(* Amplitudes *)

M1[s1_,s2_]:=-gZ Sum[(UmBar[w1,w2,m,s1] . UBarGammaUGammaU[r1,r2,-1,i1]  UBarGammaUGammaUBar[r1,r2,-1,i1] . gvga . Vm[t1,t2,m,s2]),{i1,1,2}];
  
M1c[s1_,s2_]:=-gZ Sum[(VmBar[t1,t2,m,s2] . UBarGammaUGammaU[r1,r2,-1,i1]  UBarGammaUGammaUBar[r1,r2,-1,i1] . gvga . Um[w1,w2,m,s1]),{i1,1,2}];

res[s1_,s2_]:=DiracSimplify[M1[s1,s2]];
resc[s1_,s2_]:=DiracSimplify[M1c[s1,s2]];

vlist={w1,w2,t1,t2,r1,r2}

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

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



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



Maux[s1_,s2_]:=Expand[res[s1,s2] /. simp]
Mauxc[s1_,s2_]:=Expand[resc[s1,s2] /. simp]

M[s1_,s2_]:= Expand[Maux[s1,s2]]
Mc[s1_,s2_]:=Expand[Mauxc[s1,s2]]

(* Change to True to write OutputFortran *)

If[False,
ANS=Expand[M[1,1]*Mc[1,1]+M[-1,1]*Mc[-1,1]+M[1,-1]*Mc[1,-1]+M[-1,-1]*Mc[-1,-1]];


simp4={sp[p_,q_] spc[p_,q_]->2 dot[p,q]};


ANS=ANS /. simp4;


simp5={dot[t1,t2]->m^2/2,dot[t2,t1]->m^2/2,dot[w1,w2]->m^2/2,dot[w2,w1]->m^2/2};


ANS=ANS /. simp5 ;
ANS=ANS /. simp4 ;
ANS=ANS /. simp4 ;
ANS=ANS /. simp4 ;
ANS=ANS /. simp5 ;


SetOptions[$Output,PageWidth->65];
(* << Format.m *)


mpp=M[1,1];
mpm=M[1,-1];
mmp=M[-1,1];
mmm=M[-1,-1];


stmp=OpenWrite["ZfF.f",FormatType -> FortranForm];
Write[stmp,"MPP=",mpp];
Write[stmp,"MPM=",mpm];
Write[stmp,"MMP=",mmp];
Write[stmp,"MMM=",mmm];
Close[stmp];
]
