(* ------------------------------------------------------------
Input program for FeynCalc for the e^- e^+ -> Higgs gamma. 
The notation and other explanations can be  found at the page:

http://porthos.ist.utl.pt/CTQFT

See also Nuclear Physics B267(1986)509

Author: Jorge C. Romão
email: jorge.romao@ist.utl.pt
---------------------------------------------------------------*)
(* Definitons *)
Remove["Global`*"] 
dm[mu_]:=DiracMatrix[mu]
dm[5]:=DiracMatrix[5]
ds[p_]:=DiracSlash[p]
mt[mu_,nu_]:=MetricTensor[mu,nu]
fv[p_,mu_]:=FourVector[p,mu]
epsilon[a_,b_,c_,d_]:=LeviCivita[a,b,c,d]
id[n_]:=IdentityMatrix[n]
sp[p_,q_]:=ScalarProduct[p,q]
li[mu_]:=LorentzIndex[mu]
prop[p_,m_]:=m + ds[p]

PV[p_,mu_]:=PolarizationVector[p,mu]
WMWPA[p_,k_,q_,a_,b_,m_]:=mt[a,b] (fv[p,m]-fv[k,m]) +\
 mt[b,m] (fv[k,a]-fv[q,a]) + mt[m,a] (fv[q,b]-fv[p,b])
VScalar[p_,k_,m_]:=fv[p-k,m]


(* Kinematics *)
kin={ScalarProduct[k,q]->(q2-mH^2)/2}
kinxs={sp[pm,pp]->s/2,sp[k,pp]->(s-2 kpm -mH^2)/2,sp[k,pm]->kpm,q2->s}

(* --------------------------------------

       The ZHG 3-point Function 

----------------------------------------*)

(* Diagrams F1 and F1a *)

numF1ZHG:=Tr[prop[p,mf] . dm[nu] . (gv -ga dm[5]) . prop[p-q,mf] . prop[p-k,mf] . dm[mu]]
numF1aZHG:=Tr[prop[-p,mf] . dm[mu] . prop[-p+k,mf] .  prop[-p+q,mf] . dm[nu] . (gv-ga dm[5])]


ampFZHG:= Qf mf/cw2 (1/2 ) Contract[(numF1ZHG+numF1aZHG) PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mf], \
PropagatorDenominator[p-q,mf],PropagatorDenominator[p-k,mf]]


onshell={sp[k,k]->0, sp[q,q]->q2}

resFZHG:= (-I/Pi^2) OneLoop[p,ampFZHG] /. onshell


ansFZHG=PaVeReduce[resFZHG]


aux =(Pair[LorentzIndex[nu], Momentum[
       Polarization[k, I]]]*
     Pair[Momentum[k], Momentum[q]] - 
    Pair[LorentzIndex[nu], Momentum[k]]*
     Pair[Momentum[q], Momentum[
       Polarization[k, I]]])

aux5=PolarizationVector[k,nu]
aux2=FourVector[k,nu]  Pair[Momentum[q], Momentum[
       Polarization[k, I]]]
aux3=FourVector[q,nu]  Pair[Momentum[q], Momentum[
       Polarization[k, I]]]


CFZHG=Coefficient[ansFZHG,aux] 

TestCFZHG:= Simplify[ansFZHG - (aux) CFZHG]

(* Diagrams G8 e G8a *)

numG8ZHG:= (1/2 mH^2) (-1/2) (1-tw2) VScalar[p,k-p,mu] VScalar[p-q,-p,nu]
numG8aZHG:=(1/2 mH^2) (-1/2) (1-tw2) VScalar[-p+k,p,mu] VScalar[-p,p-q,nu]


ampG8ZHG:=Contract[(numG8ZHG+numG8aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG8ZHG:= (-I/Pi^2) OneLoop[p,ampG8ZHG] /. onshell


(* Diagrams G13 *)

numG13ZHG:= - mH^2 (-1/2) (1-tw2) mt[mu,nu]


ampG13ZHG:=Contract[numG13ZHG PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-q+k,mW]]


resG13ZHG:= (-I/Pi^2) OneLoop[p,ampG13ZHG] /. onshell



(* Diagrams G8+G8a+G13 *)

ansG8ZHG=PaVeReduce[resG8ZHG] 

ansG13ZHG=PaVeReduce[resG13ZHG] 

ansG8G13ZHG=Simplify[ansG8ZHG+ansG13ZHG  /. q2-> mH^2 + 2 ScalarProduct[q,k]]

auxG8G13ZHG = Coefficient[ansG8G13ZHG,aux] 

TestG8G13ZHG:= Simplify[ansG8G13ZHG - (aux) auxG8G13ZHG]

XG8G13ZHG:= Simplify[auxG8G13ZHG /. ScalarProduct[q,k]-> (q2-mH^2)/2]


(* Diagram G1+G1a *)

numG1ZHG:= (mW^2) (-1) WMWPA[-p,p-q,q,a,b,nu] WMWPA[-p+k,p,-k,b,a,mu]
numG1aZHG:= (mW^2) (-1) WMWPA[p-q,-p,q,b,a,nu] WMWPA[p,-p+k,-k,a,b,mu]


ampG1ZHG:=Contract[(numG1ZHG+numG1aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG1ZHG:= (-I/Pi^2) OneLoop[p,ampG1ZHG] /. onshell

ansG1ZHG=PaVeReduce[resG1ZHG] 



(* Diagram G2+G2a *)

numG2ZHG:= (-mW^4) tw2 mt[mu,nu]
numG2aZHG:= (-mW^4) tw2 mt[mu,nu]


ampG2ZHG:=Contract[(numG2ZHG+numG2aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG2ZHG:= (-I/Pi^2) OneLoop[p,ampG2ZHG] /. onshell

ansG2ZHG=PaVeReduce[resG2ZHG] 



(* Diagram G3+G3a  *)

numG3ZHG:= (-1/2 mW^2) (-1) VScalar[p-k,k-q,b] \
WMWPA[-p,p-q,q,mu,b,nu] 

numG3aZHG:= (1/2 mW^2) (-1) VScalar[p-k,k-q,b]\
WMWPA[p-q,-p,q,b,mu,nu] 


ampG3ZHG:=Contract[(numG3ZHG+numG3aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG3ZHG:= (-I/Pi^2) OneLoop[p,ampG3ZHG] /. onshell

ansG3ZHG=PaVeReduce[resG3ZHG] 


(* Diagram G4+G4a  *)

numG4ZHG:= (1/2 mW^2) tw2 VScalar[-p+q,k-q,b] \
WMWPA[-p+k,p,-k,b,nu,mu] 

numG4aZHG:= (-1/2 mW^2) tw2 VScalar[-p+q,k-q,b] \
WMWPA[p,-p+k,-k,nu,b,mu] 


ampG4ZHG:=Contract[(numG4ZHG+numG4aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG4ZHG:= (-I/Pi^2) OneLoop[p,ampG4ZHG] /. onshell

ansG4ZHG=PaVeReduce[resG4ZHG] 


(* Diagram G5+G5a  *)

numG5ZHG:= (-1/2 mW^2 ) (-1/2) (1-tw2) VScalar[p-q,-p,nu] VScalar[-p+q,k-q,mu] 

numG5aZHG:= (1/2 mW^2 ) (-1/2) (1-tw2) VScalar[-p,p-q,nu] VScalar[-p+q,k-q,mu] 


ampG5ZHG:=Contract[(numG5ZHG+numG5aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG5ZHG:= (-I/Pi^2) OneLoop[p,ampG5ZHG] /. onshell

ansG5ZHG=PaVeReduce[resG5ZHG] 


(* Diagram G6+G6a  *)

numG6ZHG:= (1/2 mW^2) tw2 VScalar[p,-p+k,mu] VScalar[p-k,k-q,nu]

numG6aZHG:= (-1/2 mW^2 ) tw2  VScalar[-p+k,p,mu] VScalar[p-k,k-q,nu] 



ampG6ZHG:=Contract[(numG6ZHG+numG6aZHG) PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG6ZHG:= (-I/Pi^2) OneLoop[p,ampG6ZHG] /. onshell

ansG6ZHG=PaVeReduce[resG6ZHG] 


(* Diagram G7+G7a  OK *)

numG7ZHG:= (-1/2 mW^2 mH^2) tw2 mt[mu,nu]

numG7aZHG:= numG7ZHG 



ampG7ZHG:=Contract[(numG7ZHG+numG7aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG7ZHG:= (-I/Pi^2) OneLoop[p,ampG7ZHG] /. onshell

ansG7ZHG=PaVeReduce[resG7ZHG] 


(* Diagram G9+G9a  *)

numG9ZHG:= 2 ( - 1/2 mW^2) (-1) fv[p,nu] fv[p-k,mu]

numG9aZHG:= 2 ( - 1/2 mW^2) (-1) fv[-p,mu] fv[-p+q,nu]


ampG9ZHG:=Contract[(numG9ZHG+numG9aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG9ZHG:= (-I/Pi^2) OneLoop[p,ampG9ZHG] /. onshell

ansG9ZHG=PaVeReduce[resG9ZHG] 


(* Diagrams G10 and G10a *)

numG10ZHG:= (-1/2 mW^2) tw2  mt[mu,nu]
numG10aZHG:= (-1/2 mW^2) tw2 mt[mu,nu]


ampG10ZHG:=Contract[(numG10ZHG +numG10aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW]]


resG10ZHG:= (-I/Pi^2) OneLoop[p,ampG10ZHG] /. onshell


ansG10ZHG=PaVeReduce[resG10ZHG] 


(* Diagrams G11 and G11a *)

numG11ZHG:= (-1/2 mW^2) tw2 mt[mu,nu]
numG11aZHG:= (-1/2 mW^2) tw2 mt[mu,nu]


ampG11ZHG:=Contract[(numG11ZHG +numG11aZHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-q,mW]]


resG11ZHG:= (-I/Pi^2) OneLoop[p,ampG11ZHG] /. onshell


ansG11ZHG=PaVeReduce[resG11ZHG] 


(* Diagrams G12 and G12a *)

numG12ZHG:=(-1) (-mW^2) mt[a,b] (2 mt[a,b] mt[mu,nu]- mt[a,mu] mt[b,nu]\
 -mt[a,nu] mt[b,mu])

ampG12ZHG:=Contract[numG12ZHG PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-q+k,mW]]


resG12ZHG:= (-I/Pi^2) OneLoop[p,ampG12ZHG] /. onshell


ansG12ZHG=PaVeReduce[resG12ZHG] 


(* Test that the divergences cancell out *)

div={B0[m1_,m2_,m3_]->Div}

ansCounterTermZHG= mW^2 (1+ tw2) aux5 2 B0[0,mW^2,mW^2]

ansZHG=ansG1ZHG+ansG2ZHG+ansG3ZHG+ansG4ZHG+ansG5ZHG+ansG6ZHG+ansG7ZHG + \
       ansG9ZHG+ansG10ZHG+ansG11ZHG+ansG12ZHG+ansCounterTermZHG

ansGZHG=ansZHG+ansG8G13ZHG 

ansGZHGkin= ansGZHG /. kin

TestDivZHG:= Simplify[Coefficient[ansGZHGkin /.div ,Div]]

TestDivDiagZHG=Function[exp,test=exp /. kin; test=Simplify[Coefficient[test /.div ,Div]];test=Coefficient[test,aux5];test=Coefficient[test,mW^2]];


(* Separation of Gauge Invariant Part *)

(*------------------------------------------------------------------------

ansGZHG = d5 PV[k,nu] + d2 sp[PV(k),q] fv[k,nu] + d3 sp[PV(k),q] fv[q,nu] 

We want to show that 

 ansGZHG =  d1 [ sp[k,q] PV[k,nu] - sp[PV(k),q] fv[k,nu] ] 

         + d3  sp[PV(k),q] fv[q,nu] - d4 (q^2-mZ^2)  PV[k,nu] ]
 Solution:

   d1 = -d2 

   -d4 (q^2 -mZ^2) = d5 - d1 sp[k,q]
   

   Test :  d1 sp[k,q] = d5 + d4 (q^2-mZ^2)

--------------------------------------------------------------------------*)

d5ZHG = Simplify[Coefficient[ansGZHG, aux5] ];
d2ZHG = Simplify[Coefficient[ansGZHG, aux2] ];
d3ZHG = Simplify[Coefficient[ansGZHG, aux3] ];


d1ZHG = -d2ZHG 
d4auxZHG = Simplify[d5ZHG - d1ZHG sp[k,q] /. kin ]
d4ZHG= -Coefficient[d4auxZHG,q2]

Testd4:= Simplify[d4auxZHG + d4ZHG*(q2 -mW^2 (1+tw2))]



(* --------------------------------------

       The GHG 3-point Function 

----------------------------------------*)

(* Diagrams F1 and F1a *)

numF1GHG:=Tr[prop[p,mf] . dm[nu] . prop[p-q,mf] . prop[p-k,mf] .dm[mu]]
numF1aGHG:=Tr[prop[-p,mf] . dm[mu] . prop[-p+k,mf] .  prop[-p+q,mf] . dm[nu]]


ampFGHG:= Qf^2 mf (-1/2 ) Contract[(numF1GHG+numF1aGHG) PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mf], \
PropagatorDenominator[p-q,mf],PropagatorDenominator[p-k,mf]]


onshell={sp[k,k]->0, sp[q,q]->q2}

resFGHG:= (-I/Pi^2) OneLoop[p,ampFGHG] /. onshell


ansFGHG=PaVeReduce[resFGHG]


aux =(Pair[LorentzIndex[nu], Momentum[
       Polarization[k, I]]]*
     Pair[Momentum[k], Momentum[q]] - 
    Pair[LorentzIndex[nu], Momentum[k]]*
     Pair[Momentum[q], Momentum[
       Polarization[k, I]]])

aux5=PolarizationVector[k,nu]
aux2=FourVector[k,nu]  Pair[Momentum[q], Momentum[
       Polarization[k, I]]]
aux3=FourVector[q,nu]  Pair[Momentum[q], Momentum[
       Polarization[k, I]]]


CFGHG=Coefficient[ansFGHG,aux] 


(* Diagrams G8 e G8a *)

numG8GHG:= (1/2 mH^2) VScalar[p,k-p,mu] VScalar[p-q,-p,nu]
numG8aGHG:=(1/2 mH^2) VScalar[-p+k,p,mu] VScalar[-p,p-q,nu]


ampG8GHG:=Contract[(numG8GHG+numG8aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG8GHG:= (-I/Pi^2) OneLoop[p,ampG8GHG] /. onshell


(* Diagrams G13 *)

numG13GHG:= - mH^2 mt[mu,nu]


ampG13GHG:=Contract[numG13GHG PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-q+k,mW]]


resG13GHG:= (-I/Pi^2) OneLoop[p,ampG13GHG] /. onshell



(* Diagrams G8+G8a+G13 *)

ansG8GHG=PaVeReduce[resG8GHG] 

ansG13GHG=PaVeReduce[resG13GHG] 

ansG8G13GHG=Simplify[ansG8GHG+ansG13GHG  /. q2-> mH^2 + 2 ScalarProduct[q,k]]

auxG8G13GHG = Coefficient[ansG8G13GHG,aux] 

TestG8G13GHG:= Simplify[ansG8G13GHG - (aux) auxG8G13GHG]

XG8G13GHG:= Simplify[auxG8G13GHG /. ScalarProduct[q,k]-> (q2-mH^2)/2]


(* Diagram G1+G1a *)

numG1GHG:= (mW^2) WMWPA[-p,p-q,q,a,b,nu] WMWPA[-p+k,p,-k,b,a,mu]
numG1aGHG:= (mW^2) WMWPA[p-q,-p,q,b,a,nu] WMWPA[p,-p+k,-k,a,b,mu]


ampG1GHG:=Contract[(numG1GHG+numG1aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG1GHG:= (-I/Pi^2) OneLoop[p,ampG1GHG] /. onshell

ansG1GHG=PaVeReduce[resG1GHG] 



(* Diagram G2+G2a *)

numG2GHG:= (-mW^4) mt[mu,nu]
numG2aGHG:= (-mW^4) mt[mu,nu]


ampG2GHG:=Contract[(numG2GHG+numG2aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG2GHG:= (-I/Pi^2) OneLoop[p,ampG2GHG] /. onshell

ansG2GHG=PaVeReduce[resG2GHG] 


(* Diagram G3+G3a  *)

numG3GHG:= (-1/2 mW^2) VScalar[p-k,k-q,b] \
WMWPA[-p,p-q,q,mu,b,nu] 

numG3aGHG:= (1/2 mW^2) VScalar[p-k,k-q,b]\
WMWPA[p-q,-p,q,b,mu,nu] 


ampG3GHG:=Contract[(numG3GHG+numG3aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG3GHG:= (-I/Pi^2) OneLoop[p,ampG3GHG] /. onshell

ansG3GHG=PaVeReduce[resG3GHG] 


(* Diagram G4+G4a  *)

numG4GHG:= (1/2 mW^2) VScalar[-p+q,k-q,b] \
WMWPA[-p+k,p,-k,b,nu,mu] 

numG4aGHG:= (-1/2 mW^2) VScalar[-p+q,k-q,b] \
WMWPA[p,-p+k,-k,nu,b,mu] 


ampG4GHG:=Contract[(numG4GHG+numG4aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG4GHG:= (-I/Pi^2) OneLoop[p,ampG4GHG] /. onshell

ansG4GHG=PaVeReduce[resG4GHG] 


(* Diagram G5+G5a  *)

numG5GHG:= (-1/2 mW^2 ) VScalar[p-q,-p,nu] VScalar[-p+q,k-q,mu] 

numG5aGHG:= (1/2 mW^2 ) VScalar[-p,p-q,nu] VScalar[-p+q,k-q,mu] 


ampG5GHG:=Contract[(numG5GHG+numG5aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG5GHG:= (-I/Pi^2) OneLoop[p,ampG5GHG] /. onshell

ansG5GHG=PaVeReduce[resG5GHG] 


(* Diagram G6+G6a  *)

numG6GHG:= (1/2 mW^2) VScalar[p,-p+k,mu] VScalar[p-k,k-q,nu]

numG6aGHG:= (-1/2 mW^2 )  VScalar[-p+k,p,mu] VScalar[p-k,k-q,nu] 



ampG6GHG:=Contract[(numG6GHG+numG6aGHG) PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG6GHG:= (-I/Pi^2) OneLoop[p,ampG6GHG] /. onshell

ansG6GHG=PaVeReduce[resG6GHG] 


(* Diagram G7+G7a  OK *)

numG7GHG:= (-1/2 mW^2 mH^2) mt[mu,nu]

numG7aGHG:= numG7GHG 



ampG7GHG:=Contract[(numG7GHG+numG7aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG7GHG:= (-I/Pi^2) OneLoop[p,ampG7GHG] /. onshell

ansG7GHG=PaVeReduce[resG7GHG] 


(* Diagram G9+G9a  *)

numG9GHG:= 2 ( - 1/2 mW^2) fv[p,nu] fv[p-k,mu]

numG9aGHG:= 2 ( - 1/2 mW^2) fv[-p,mu] fv[-p+q,nu]


ampG9GHG:=Contract[(numG9GHG+numG9aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW],PropagatorDenominator[p-q,mW]]


resG9GHG:= (-I/Pi^2) OneLoop[p,ampG9GHG] /. onshell

ansG9GHG=PaVeReduce[resG9GHG] 


(* Diagrams G10 and G10a *)

numG10GHG:= (-1/2 mW^2) mt[mu,nu]
numG10aGHG:= (-1/2 mW^2) mt[mu,nu]


ampG10GHG:=Contract[(numG10GHG +numG10aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW]]


resG10GHG:= (-I/Pi^2) OneLoop[p,ampG10GHG] /. onshell


ansG10GHG=PaVeReduce[resG10GHG] 


(* Diagrams G11 and G11a *)

numG11GHG:= (-1/2 mW^2) mt[mu,nu]
numG11aGHG:= (-1/2 mW^2) mt[mu,nu]


ampG11GHG:=Contract[(numG11GHG +numG11aGHG) PV[k,mu]  ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-q,mW]]


resG11GHG:= (-I/Pi^2) OneLoop[p,ampG11GHG] /. onshell


ansG11GHG=PaVeReduce[resG11GHG] 


(* Diagrams G12 and G12a *)

numG12GHG:= (-mW^2) mt[a,b] (2 mt[a,b] mt[mu,nu]- mt[a,mu] mt[b,nu]\
 -mt[a,nu] mt[b,mu])

ampG12GHG:=Contract[numG12GHG PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-q+k,mW]]


resG12GHG:= (-I/Pi^2) OneLoop[p,ampG12GHG] /. onshell


ansG12GHG=PaVeReduce[resG12GHG] 


(* Test that the divergences cancell out *)

div={B0[m1_,m2_,m3_]->Div}

ansGHG=ansG1GHG+ansG2GHG+ansG3GHG+ansG4GHG+ansG5GHG+ansG6GHG+ansG7GHG + \
       ansG9GHG+ansG10GHG+ansG11GHG+ansG12GHG 

ansGGHG=ansGHG+ansG8G13GHG 

ansGkinGHG= ansGGHG /. kin

TestDivGHG:= Simplify[Coefficient[ansGkinGHG /.div ,Div]]

TestDivDiagGHG=Function[exp,test=exp /. kin; test=Simplify[Coefficient[test /.div ,Div]]; test=Coefficient[test,aux5 mW^2]];


(* Separation of Gauge Invariant Part *)
(* 

ansG = c5 PV[k,nu] + c2 sp[PV(k),q] fv[k,nu] + c3 sp[PV(k),q] fv[q,nu] 

We want to show that 

 ansG =  c3 [ sp[PV(k),q] fv[q,nu] - q^2 PV[k,nu] ]

       + c1 [ sp[k,q] PV[k,nu] - sp[PV(k),q] fv[q,nu] ]

 Solution:

   c1 = -c2 

   Test :  c1 sp[k,q] = c5 + c3 q^2


*)


c5GHG = Simplify[Coefficient[ansGGHG, aux5] ];
c2GHG = Simplify[Coefficient[ansGGHG, aux2] ];
c3GHG = Simplify[Coefficient[ansGGHG, aux3] ];


c1GHG = -c2GHG 

TestGS:= Simplify[c1GHG sp[k,q] - c5GHG - c3GHG q2 /. kin]



(* --------------------------------------

       The Boxes

----------------------------------------*)

gpgm:=(g6 dm[6] + g7 dm[7])


(* Diagram B4 *)

numB4:= Spinor[pp,0] . dm[a] . gpgm . prop[p-k,0] . dm[mu] . ds[p] . dm[a] . gpgm . Spinor[pm,0] 


ampB4:=  Contract[ numB4 PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,me], \
PropagatorDenominator[p-k,me], \
PropagatorDenominator[p-k+pp,mZ], \
PropagatorDenominator[p-pm,mZ]]


onshell={sp[k,k]->0, sp[pp,pp]->0, sp[pm,pm]->0}

resB4:= (-I/Pi^2) OneLoop[p,ampB4] /. onshell

ansB4=PaVeReduce[resB4]

var = Select[Variables[ansB4], (Head[#] === StandardMatrixElement) &]


(* Definition of the Matrix Elements

 ME[1] = u(pp) . dm[mu] . dm[6] . u(pm)  PV[k,mu]
 ME[2] = u(pp) . dm[mu] . dm[7] . u(pm)  PV[k,mu]
 ME[3] = u(pp) . dm[mu] . dm[6] . u(pm)  fv[k,mu] sp[pm,PV[k]]
 ME[4] = u(pp) . dm[mu] . dm[7] . u(pm)  fv[k,mu] sp[pm,PV[k]]
 ME[5] = u(pp) . dm[mu] . dm[6] . u(pm)  fv[k,mu] sp[pp,PV[k]]
 ME[6] = u(pp) . dm[mu] . dm[7] . u(pm)  fv[k,mu] sp[pp,PV[k]]

 *)



  Set @@ {var, {ME[1],ME[2],ME[3],ME[4],ME[5],ME[6]}} 

B4c1=Coefficient[ansB4,ME[1]]
B4c2=Coefficient[ansB4,ME[2]]
B4c3=Coefficient[ansB4,ME[3]]
B4c4=Coefficient[ansB4,ME[4]]
B4c5=Coefficient[ansB4,ME[5]]
B4c6=Coefficient[ansB4,ME[6]]

TestB4:=Simplify[ansB4 - B4c1* ME[1]- B4c2* ME[2]- B4c3* ME[3]- B4c4* ME[4]- B4c5* ME[5]- B4c6* ME[6]]


(* Diagram B5 *)

numB5:= Spinor[pp,0] . dm[a] . gpgm . prop[p-pp,0] . dm[a] . gpgm . prop[pm-k,0] . dm[mu] . Spinor[pm,0] 


ampB5:=  Contract[ numB5 PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mZ], \
PropagatorDenominator[p+k-pp-pm,mZ], \
PropagatorDenominator[p-pp,me], \
PropagatorDenominator[pm-k,0]]


resB5:= (-I/Pi^2) OneLoop[p,ampB5] /. onshell

ansB5=PaVeReduce[resB5]


B5c1=Coefficient[ansB5,ME[1]]
B5c2=Coefficient[ansB5,ME[2]]
B5c3=Coefficient[ansB5,ME[3]]
B5c4=Coefficient[ansB5,ME[4]]
B5c5=Coefficient[ansB5,ME[5]]
B5c6=Coefficient[ansB5,ME[6]]

TestB5:=Simplify[ansB5 - B5c1* ME[1]- B5c2* ME[2]- B5c3* ME[3]- B5c4* ME[4]- B5c5* ME[5]- B5c6* ME[6]]




(* Diagram B6 *)

numB6:= Spinor[pp,0] . dm[mu] . prop[-pp+k,0] . dm[a] . gpgm . prop[p+pm,0] . dm[a] . gpgm . Spinor[pm,0] 


ampB6:=  Contract[ numB6 PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mZ], \
PropagatorDenominator[p-k+pp+pm,mZ], \
PropagatorDenominator[p+pm,me], \
PropagatorDenominator[pp-k,0]]


resB6:= (-I/Pi^2) OneLoop[p,ampB6] /. onshell

ansB6=PaVeReduce[resB6]


B6c1=Coefficient[ansB6,ME[1]]
B6c2=Coefficient[ansB6,ME[2]]
B6c3=Coefficient[ansB6,ME[3]]
B6c4=Coefficient[ansB6,ME[4]]
B6c5=Coefficient[ansB6,ME[5]]
B6c6=Coefficient[ansB6,ME[6]]

TestB6:=Simplify[ansB6 - B6c1* ME[1]- B6c2* ME[2]- B6c3* ME[3]- B6c4* ME[4]- B6c5* ME[5]- B6c6* ME[6]]


c1Z=B4c1+B5c1+B6c1
c2Z=B4c2+B5c2+B6c2
c3Z=-(B4c3+B5c3+B6c3)
c4Z=-(B4c4+B5c4+B6c4)
c5Z=-(B4c5+B5c5+B6c5)
c6Z=-(B4c6+B5c6+B6c6);


(* Test Gauge Invariance of Z boxes*)

TestGIBoxZ:=Simplify[c1Z - c3Z sp[k,pm] - c5Z sp[k,pp] /. me->0]

(* Diagram B2 *)


numB2a:= Spinor[pp,0] . dm[a] . dm[7] . prop[p+pm,0] . dm[b] . dm[7] . Spinor[pm,0] WMWPA[-p,p+k,-k,b,a,mu]

numB2b:= Spinor[pp,0] . dm[a] . dm[7] . prop[p-pp,0] . dm[b] . dm[7] . Spinor[pm,0] WMWPA[-p+k,p,-k,b,a,mu]


ampB2a:=  Contract[ (1/2) numB2a  PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p+pm,me], \
PropagatorDenominator[p,mW], \
PropagatorDenominator[p+k,mW], \
PropagatorDenominator[p+pm+pp,mW]]

ampB2b:=  Contract[ (1/2) numB2b  PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p-pp,me], \
PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW], \
PropagatorDenominator[p-pm-pp,mW]]

resB2a:= (-I/Pi^2) OneLoop[p,ampB2a] /. onshell
resB2b:= (-I/Pi^2) OneLoop[p,ampB2b] /. onshell


ansB2=PaVeReduce[resB2a+resB2b]

B2c1=Coefficient[ansB2,ME[1]]
B2c2=Coefficient[ansB2,ME[2]]
B2c3=Coefficient[ansB2,ME[3]]
B2c4=Coefficient[ansB2,ME[4]]
B2c5=Coefficient[ansB2,ME[5]]
B2c6=Coefficient[ansB2,ME[6]]

TestB2:=Simplify[ansB2 - B2c1* ME[1]- B2c2* ME[2]- B2c3* ME[3]- B2c4* ME[4]- B2c5* ME[5]- B2c6* ME[6]]


(* Diagram B3 *)

numB3a:= Spinor[pp,0] . dm[a] . dm[7] . prop[p+pm,0] . dm[mu] . dm[7] . Spinor[pm,0] VScalar[-p-k,k-pp-pm,a]

numB3b:= Spinor[pp,0] . dm[mu] . dm[7] . prop[p-pp,0] . dm[a] . dm[7] . Spinor[pm,0] VScalar[p-k,k-pp-pm,a]


ampB3a:=  Contract[ (1/4) numB3a PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p+pm,me], \
PropagatorDenominator[p,mW], \
PropagatorDenominator[p+k,mW], \
PropagatorDenominator[p+pm+pp,mW]]
resB3a:= (-I/Pi^2) OneLoop[p,ampB3a] /. onshell

ampB3b:=  Contract[ (-1/4) numB3b PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p-pp,me], \
PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k,mW], \
PropagatorDenominator[p-pm-pp,mW]]
resB3b:= (-I/Pi^2) OneLoop[p,ampB3b] /. onshell


ansB3=PaVeReduce[resB3a+resB3b]


B3c1=Coefficient[ansB3,ME[1]]
B3c2=Coefficient[ansB3,ME[2]]
B3c3=Coefficient[ansB3,ME[3]]
B3c4=Coefficient[ansB3,ME[4]]
B3c5=Coefficient[ansB3,ME[5]]
B3c6=Coefficient[ansB3,ME[6]]

TestB3:=Simplify[ansB3 - B3c1* ME[1]- B3c2* ME[2]- B3c3* ME[3]- B3c4* ME[4]- B3c5* ME[5]- B3c6* ME[6]]


(* Diagram B7 *)

numB7:= Spinor[pp,0] . dm[a] . dm[7] . prop[p-pp,0] . dm[a] . dm[7] . prop[pm-k,0] . dm[mu] . Spinor[pm,0] 


ampB7:=  Contract[ (1/2) numB7 PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p+k-pp-pm,mW], \
PropagatorDenominator[p-pp,me], \
PropagatorDenominator[pm-k,0]]


resB7:= (-I/Pi^2) OneLoop[p,ampB7] /. onshell

ansB7=PaVeReduce[resB7]


B7c1=Coefficient[ansB7,ME[1]]
B7c2=Coefficient[ansB7,ME[2]]
B7c3=Coefficient[ansB7,ME[3]]
B7c4=Coefficient[ansB7,ME[4]]
B7c5=Coefficient[ansB7,ME[5]]
B7c6=Coefficient[ansB7,ME[6]]

TestB7:=Simplify[ansB7 - B7c1* ME[1]- B7c2* ME[2]- B7c3* ME[3]- B7c4* ME[4]- B7c5* ME[5]- B7c6* ME[6]]


(* Diagram B8 *)

numB8:= Spinor[pp,0] . dm[mu] . prop[-pp+k,0] . dm[a] . dm[7] . prop[p+pm,0] . dm[a] . dm[7] . Spinor[pm,0] 


ampB8:=  Contract[ (1/2) numB8 PV[k,mu] ] \
FeynAmpDenominator[PropagatorDenominator[p,mW], \
PropagatorDenominator[p-k+pp+pm,mW], \
PropagatorDenominator[p+pm,me], \
PropagatorDenominator[pp-k,0]]


resB8:= (-I/Pi^2) OneLoop[p,ampB8] /. onshell

ansB8=PaVeReduce[resB8]


B8c1=Coefficient[ansB8,ME[1]]
B8c2=Coefficient[ansB8,ME[2]]
B8c3=Coefficient[ansB8,ME[3]]
B8c4=Coefficient[ansB8,ME[4]]
B8c5=Coefficient[ansB8,ME[5]]
B8c6=Coefficient[ansB8,ME[6]]

TestB8:=Simplify[ansB8 - B8c1* ME[1]- B8c2* ME[2]- B8c3* ME[3]- B8c4* ME[4]- B8c5* ME[5]- B8c6* ME[6]]


c1W=B2c1+B3c1+B7c1+B8c1
c2W=B2c2+B3c2+B7c2+B8c2
c3W=-(B2c3+B3c3+B7c3+B8c3)
c4W=-(B2c4+B3c4+B7c4+B8c4)
c5W=-(B2c5+B3c5+B7c5+B8c5)
c6W=-(B2c6+B3c6+B7c6+B8c6);



(* Get Gauge Non-Invariant part of W boxes*)

GNIBoxW:=Simplify[c2W - c4W sp[k,pm] - c6W sp[k,pp]]

kinxs={sp[pm,pp]->s/2,sp[k,pp]->(s-2 kpm -mH2)/2,sp[k,pm]->kpm,q2->s}


(* Gauge Invariance Tests *)

TestGIBoxW:= Simplify[Simplify[GNIBoxW*mW + Simplify[c3GHG (-1/2/mW) /. kin] /. kinxs] /. {mH^2->mH2,me->0}]


(* Final Coefficients *)


AL[mu_, rho_] :=  Sum[al[i] ((sp[k, pp] mt[mu, rho] - fv[k, rho] \
  fv[pp, mu]) Gp[i] + (sp[k, pm] mt[mu, rho] - fv[k, rho] fv[pm, mu]) \
  Gm[i]), {i, 1, 4}]

AR[mu_, rho_] :=  Sum[ar[i] ((sp[k, pp] mt[mu, rho] - fv[k, rho] \
  fv[pp, mu]) Gp[i] + (sp[k, pm] mt[mu, rho] - fv[k, rho] fv[pm, mu]) \
  Gm[i]), {i, 1, 4}]

ALc[mu_, rho_] :=  Sum[al[i] ((sp[k, pp] mt[mu, rho] - fv[k, rho] \
  fv[pp, mu]) Gpc[i] + (sp[k, pm] mt[mu, rho] - fv[k, rho] fv[pm, mu]) \
  Gmc[i]), {i, 1, 4}]

ARc[mu_, rho_] :=  Sum[ar[i] ((sp[k, pp] mt[mu, rho] - fv[k, rho] \
  fv[pp, mu]) Gpc[i] + (sp[k, pm] mt[mu, rho] - fv[k, rho] fv[pm, mu]) \
  Gmc[i]), {i, 1, 4}]


ar[3]=0

res:=Contract[ Tr[ds[pp] . dm[rho] . (AL[mu,rho] dm[7] + AR[mu,rho] dm[6]) . ds[pm] . dm[a] . (ALc[nu,a] dm[7] + ARc[nu,a] dm[6])] (- mt[mu,nu])]

simp={ar[1]->al[1],ar[3]->0,Gm[1]->Gp[1],Gmc[1]->Gpc[1],Gm[2]->Gp[2],Gmc[2]->Gpc[2]}

simp2={Gp[i_]->ReGp[i] + I ImGp[i],Gpc[i_]->ReGp[i] - I ImGp[i],Gm[i_]->ReGm[i] + I ImGm[i],Gmc[i_]->ReGm[i] - I ImGm[i]}

simp3={ReGp[1]->ReGp1,ReGp[2]->ReGp2,ReGp[3]->ReGp3,ReGp[4]->ReGp4,ReGm[3]->ReGm3,ReGm[4]->ReGm4,ImGp[1]->ImGp1,ImGp[2]->ImGp2,ImGp[3]->ImGp3,ImGp[4]->ImGp4,ImGm[3]->ImGm3,ImGm[4]->ImGm4}


res= res /. onshell;
res= res /. kinxs
res= Simplify[res /. simp]
res= Simplify[res/mW2^3 /. simp2]
res= Simplify[res /. simp3]

(* Separate the contributions from various Diagrams *)

TG=Simplify[res /. {al[2]->0,ar[2]->0,al[3]->0,ar[3]->0,al[4]->0,ar[4]->0}]
TZ=Simplify[res /. {al[1]->0,al[3]->0,ar[3]->0,al[4]->0,ar[4]->0}]
TGTZ=Simplify[res -TG -TZ /. {al[3]->0,ar[3]->0,al[4]->0,ar[4]->0}]
BW=Simplify[res /. {al[1]->0,al[2]->0,ar[2]->0,al[4]->0,ar[4]->0}]
BZ=Simplify[res /. {al[1]->0,al[2]->0,ar[2]->0,al[3]->0,ar[3]->0}]
BWBZ=Simplify[res - BW -BZ /. {al[1]->0,al[2]->0,ar[2]->0}]
TGZBWZ=Simplify[res- TG-TZ-TGTZ-BW-BZ-BWBZ]
TestSeparation:=Simplify[res-TG-TZ-TGTZ-BW-BZ-BWBZ-TGZBWZ]


Gp3=Simplify[mW^4 c6W /. kinxs]
Gm3=Simplify[mW^4 c4W /. kinxs]
Gp4=Simplify[mZ^4 c5Z/g6^2 /. kinxs]
Gm4=Simplify[mZ^4 c3Z/g6^2 /. kinxs]

Gp3=Coefficient[Gp3,me,0] /. me->eps
Gp4=Coefficient[Gp4,me,0] /. me->eps
Gm3=Coefficient[Gm3,me,0] /. me->eps
Gm4=Coefficient[Gm4,me,0] /. me->eps


SetOptions[$Output,PageWidth->62]

fermions={Qf->Qf[i],gv->gv[i],ga->ga[i],mf->mf[i]}

CFGHGaux= Simplify[Simplify[Simplify[CFGHG /.kin] /.kinxs] /. fermions]
CFZHGaux= Simplify[Simplify[Simplify[CFZHG /.kin] /.kinxs] /. fermions]
c1GHGaux= Simplify[Simplify[c1GHG /.kin] /.kinxs]
d1ZHGaux= Simplify[Simplify[d1ZHG /.kin] /.kinxs]

stmp=OpenWrite["functions.f",FormatType -> FortranForm]
Write[stmp,"TG=",TG]
Write[stmp,"TZ=",TZ]
Write[stmp,"TGTZ=",TGTZ]
Write[stmp,"BW=",BW]
Write[stmp,"BZ=",BZ]
Write[stmp,"BWBZ=",BWBZ]
Write[stmp,"TGZBWZ=",TGZBWZ]
Write[stmp,"Gp3=",Gp3]
Write[stmp,"Gm3=",Gm3]
Write[stmp,"Gp4=",Gp4]
Write[stmp,"Gm4=",Gm4]
Write[stmp,"CFGHG=",CFGHGaux]
Write[stmp,"CFZHG=",CFZHGaux]
Write[stmp,"c1GHG=",c1GHGaux]
Write[stmp,"d1ZHG=",d1ZHGaux]
Close[stmp]


