(* ------------------------------------------------------------
Input program for FeynCalc for the process H-> gamma gamma. The notation 
and other explanations can be found at the page:

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

Author: Jorge C. Romão
email: jorge.romao@ist.utl.pt
---------------------------------------------------------------*)

(* Definitons *)
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]


(* Diagrams F1 and F1a *)

numF1:=Tr[prop[k,mf] . prop[k-q1-q2,mf] . dm[nu] . prop[k-q1,mf] .dm[mu]]
numF1a:=Tr[prop[-k,mf] . dm[mu] . prop[-k+q1,mf] . dm[nu] . prop[-k+q1+q2,mf]]


ampF:=Contract[(numF1+numF1a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mf], \
PropagatorDenominator[k-q1-q2,mf],PropagatorDenominator[k-q1,mf]]


onshell={sp[q1,q1]->0, sp[q2,q2]->0,sp[p,p]->mH^2}
kin={sp[q1,q2]->mH^2/2}

resF:= (-I/Pi^2) OneLoop[k,ampF] /. onshell


ansF=PaVeReduce[resF]


aux =(Pair[Momentum[q1], Momentum[Polarization[q2, I]]]*
     Pair[Momentum[q2], Momentum[Polarization[q1, I]]] -
    Pair[Momentum[q1], Momentum[q2]]*Pair[Momentum[Polarization[q1, I]], 
      Momentum[Polarization[q2, I]]])

aux1=Pair[Momentum[q1], Momentum[Polarization[q2, I]]]*
 Pair[Momentum[q2], Momentum[Polarization[q1, I]]]

aux2=Pair[Momentum[Polarization[q1, I]],Momentum[Polarization[q2, I]]]


(* Diagrams G8 e G8a *)

numG8:= (1/2 mH^2) VScalar[k,-k+q1,mu] VScalar[k-q1,-k+q1+q2,nu]
numG8a:=(1/2 mH^2) VScalar[-k+q1,k,mu] VScalar[-k+q1+q2,k-q1,nu]


ampG8:=Contract[(numG8+numG8a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1-q2,mW],PropagatorDenominator[k-q1,mW]]


resG8:= (-I/Pi^2) OneLoop[k,ampG8] /. onshell


(* Diagrams G13 *)

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


ampG13:=Contract[numG13 PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k+q1+q2,mW]]


resG13:= (-I/Pi^2) OneLoop[k,ampG13] /. onshell


(* Diagrams G8+G8a+G13 *)

ansG8=PaVeReduce[resG8] 

ansG13=PaVeReduce[resG13] 

ansG8G13=ansG8+ansG13

auxG8G13 = -Coefficient[ansG8G13,aux] 

TestG8G13:= Simplify[ansG8G13 - (-aux) auxG8G13]

XG8G13:= auxG8G13 /. kin


(* Diagram G1+G1a *)

numG1:= (mW^2) WMWPA[-k+q1,k,-q1,b,a,mu] WMWPA[-k+q1+q2,k-q1,-q2,a,b,nu]
numG1a:= (mW^2) WMWPA[k,-k+q1,-q1,a,b,mu] WMWPA[k-q1,-k+q1+q2,-q2,b,a,nu]


ampG1:=Contract[(numG1+numG1a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1,mW],PropagatorDenominator[k-q1-q2,mW]]


resG1:= (-I/Pi^2) OneLoop[k,ampG1] /. onshell

ansG1=PaVeReduce[resG1] 

YG1:= ansG1 /. kin


(* Diagram G2+G2a *)

numG2:= (-mW^4) mt[mu,nu]
numG2a:= (-mW^4) mt[mu,nu]


ampG2:=Contract[(numG2+numG2a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1-q2,mW],PropagatorDenominator[k-q1,mW]]


resG2:= (-I/Pi^2) OneLoop[k,ampG2] /. onshell

ansG2=PaVeReduce[resG2] 

YG2:= ansG2 /. kin


(* Diagram G3+G3a  *)

numG3:= (1/2 mW^2) VScalar[-k,q1+q2,a] mt[b,mu]\
WMWPA[-k+q1+q2,k-q1,-q2,a,b,nu] 

numG3a:= (-1/2 mW^2) VScalar[-k,q1+q2,a] mt[b,mu]\
WMWPA[k-q1,-k+q1+q2,-q2,b,a,nu] 


ampG3:=Contract[(numG3+numG3a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1,mW],PropagatorDenominator[k-q1-q2,mW]]


resG3:= (-I/Pi^2) OneLoop[k,ampG3] /. onshell

ansG3=PaVeReduce[resG3] 

YG3:= ansG3 /. kin

(* Diagram G4+G4a  *)

numG4:= (-1/2 mW^2) VScalar[k-q1-q2,q1+q2,a] mt[b,nu] \
WMWPA[-k+q1,k,-q1,b,a,mu] 

numG4a:= (1/2 mW^2) VScalar[k-q1-q2,q1+q2,a] mt[b,nu]\
WMWPA[k,-k+q1,-q1,a,b,mu] 


ampG4:=Contract[(numG4+numG4a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1,mW],PropagatorDenominator[k-q1-q2,mW]]


resG4:= (-I/Pi^2) OneLoop[k,ampG4] /. onshell

ansG4=PaVeReduce[resG4] 

YG4:= ansG4 /. kin


(* Diagram G5+G5a  *)

numG5:= (-1/2 mW^2 ) VScalar[-k,q1+q2,a] VScalar[k,-k+q1,mu] mt[a,nu] 

numG5a:= (1/2 mW^2 ) VScalar[-k,q1+q2,a] VScalar[-k+q1,k,mu] mt[a,nu] 


ampG5:=Contract[(numG5+numG5a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1,mW],PropagatorDenominator[k-q1-q2,mW]]


resG5:= (-I/Pi^2) OneLoop[k,ampG5] /. onshell

ansG5=PaVeReduce[resG5] 

YG5:= ansG5 /. kin


(* Diagram G6+G6a  *)

numG6:= (1/2 mW^2) VScalar[k-q1-q2,q1+q2,a] VScalar[k-q1,-k+q1+q2,nu] mt[a,mu] 

numG6a:= (-1/2 mW^2 )  VScalar[k-q1-q2,q1+q2,a] VScalar[-k+q1+q2,k-q1,nu] \
mt[a,mu] 



ampG6:=Contract[(numG6+numG6a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1,mW],PropagatorDenominator[k-q1-q2,mW]]


resG6:= (-I/Pi^2) OneLoop[k,ampG6] /. onshell

ansG6=PaVeReduce[resG6] 

YG6:= ansG6 /. kin


(* Diagram G7+G7a  OK *)

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

numG7a:= numG7 



ampG7:=Contract[(numG7+numG7a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1,mW],PropagatorDenominator[k-q1-q2,mW]]


resG7:= (-I/Pi^2) OneLoop[k,ampG7] /. onshell

ansG7=PaVeReduce[resG7] 

YG7:= ansG7 /. kin


(* Diagram G9+G9a  *)

numG9:= 2 ( - 1/2 mW^2) fv[k-q1,mu] fv[k-q1-q2,nu]

numG9a:= 2 ( - 1/2 mW^2) fv[-k,mu] fv[-k+q1,nu]


ampG9:=Contract[(numG9+numG9a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1,mW],PropagatorDenominator[k-q1-q2,mW]]


resG9:= (-I/Pi^2) OneLoop[k,ampG9] /. onshell

ansG9=PaVeReduce[resG9] 

YG9:= ansG9 /. kin



(* Diagrams G10 and G10a *)

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


ampG10:=Contract[(numG10 +numG10a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q2,mW]]


resG10:= (-I/Pi^2) OneLoop[k,ampG10] /. onshell


ansG10=PaVeReduce[resG10] 

YG10:= ansG10 /. kin


(* Diagrams G11 and G11a *)

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


ampG11:=Contract[(numG11 +numG11a) PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k-q1,mW]]


resG11:= (-I/Pi^2) OneLoop[k,ampG11] /. onshell


ansG11=PaVeReduce[resG11] 

YG11:= ansG11 /. kin


(* Diagrams G12 and G12a *)

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

ampG12:=Contract[numG12 PV[q1,mu] PV[q2,nu] ] \
FeynAmpDenominator[PropagatorDenominator[k,mW], \
PropagatorDenominator[k+q1+q2,mW]]


resG12:= (-I/Pi^2) OneLoop[k,ampG12] /. onshell


ansG12=PaVeReduce[resG12] 

YG12:= ansG12 /. kin


(* Test that the divergences cancell out *)

div={B0[0,mW^2,mW^2]->Div,B0[mH^2,mW^2,mW^2]->Div}

ans=ansG1+ansG2+ansG3+ansG4+ansG5+ansG6+ansG7+ansG9+ansG10+ansG11+ansG12 
ansG=ans+ansG8G13 /. kin

TestDiv:= Simplify[Coefficient[ansG /.div ,Div]]

TestDivDiag=Function[exp,test=exp /. kin; test=Simplify[Coefficient[test /.div ,Div]]; test=Coefficient[test,aux2]; Coefficient[test,mW,2]]

(* Test Gauge Invariance *)

c2 = Simplify[Coefficient[ansG, aux2]/sp[q1,q2] /. kin];
c1 = Simplify[Coefficient[ansG, aux1] /. kin];

ansG= c2 

c2Diag=Function[exp,Simplify[Coefficient[exp, aux2]/sp[q1,q2] /. kin]]
c1Diag=Function[exp,Simplify[Coefficient[exp, aux1] /. kin] ] ;

TestGI:= Simplify[c1+c2]

(* Final results *)

ansF=-Coefficient[ansF,aux] mf (-1/2 ) /. kin

c0F = Coefficient[ansF,C0[0, 0, mH^2, mf^2, mf^2, mf^2],0]
c1F = Coefficient[ansF,C0[0, 0, mH^2, mf^2, mf^2, mf^2],1]
c0F = c0F /. mH->Sqrt[4 mf^2/tau]
c1F = c1F /. mH->Sqrt[4 mf^2/tau]
XF= Simplify[c0F + c1F C0[0, 0, mH^2, mf^2, mf^2, mf^2]]

c0G = Coefficient[ansG,C0[0, 0, mH^2, mW^2, mW^2, mW^2],0]
c1G = Coefficient[ansG,C0[0, 0, mH^2, mW^2, mW^2, mW^2],1]
c0G = c0G /. mH->Sqrt[4 mW^2/tau]
c1G = c1G /. mH->Sqrt[4 mW^2/tau]
XG= Simplify[c0G + c1G C0[0, 0, mH^2, mW^2, mW^2, mW^2]]


(* Comparison with NPB267(1986)509 *)

subJ[m_] := {J1 -> -m^2 C0[0, 0, mH^2, m^2, m^2, m^2], 
  J2 -> tau/4 (-1/2 - m^2 C0[0, 0, mH^2, m^2, m^2, m^2])}


XFBR=Simplify[-4 ( J1 -4 J2) /. subJ[mf]]
XGBR= Simplify[ 4 (4 J1 -( 6 + 4/tau) J2 )  /. subJ[mW]]

TestXFBPR:=Simplify[XF-XFBR]
TestXGBPR:=Simplify[XG-XGBR]

(* Comparison with Higgs Hunters Guide *)

subftau[m_] := {ftau[m] -> -2 m^2 C0[0, 0, mH^2, m^2, m^2, m^2] /tau}


XFHHG=Simplify[-2 tau (1 + (1-tau) ftau[mf]) /. subftau[mf]]
XGHHG= Simplify[ 2 + 3 tau + 3 tau (2-tau) ftau[mW]  /. subftau[mW]]

TestXFHHG:=Simplify[XF-XFHHG]
TestXGHHG:=Simplify[XG-XGHHG]

