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_] := ds[p] + m

PVL[Q_,mu_]:= FourVector[Q,mu]


V[a_,b_,mu_,p_,k_,q_]:=mt[a,b] fv[p-k,mu] + mt[b,mu] fv[k-q,a] + \
                       mt[mu,a] fv[q-p,b] 


PL:=(1 - dm[5])/2
PR:=(1 + dm[5])/2


(* incoming: p1=W-, p2=W+ outgoing: q1=W-, q2=W+ *)

(* definitions *)
onshell1={sp[p1,p1]->MW^2,sp[p2,p2]->MW^2,sp[q1,q1]->MW^2,sp[q2,q2]->MW^2}
onshell2={sp[P1,p1]->0,sp[P2,p2]->0,sp[Q1,q1]->0,sp[Q2,q2]->0}
onshell3={sp[P1,p2]->s/2/MW beta,sp[P1,q1]->s/4/MW beta (1-cteta),sp[P1,q2]->s/4/MW beta (1+cteta),sp[P2,p1]->s/2/MW beta,sp[P2,q1]->s/4/MW beta (1+cteta),sp[P2,q2]->s/4/MW beta (1-cteta),sp[Q1,p1]->s/4/MW beta (1-cteta),sp[Q1,p2]->s/4/MW beta (1+cteta),sp[Q1,q2]->s/2/MW beta,sp[Q2,p1]->s/4/MW beta (1+cteta),sp[Q2,p2]->s/4/MW beta (1-cteta),sp[Q2,q1]->s/2/MW beta}
onshell4={sp[P1,P1]->-1,sp[P2,P2]->-1,sp[Q1,Q1]->-1,sp[Q2,Q2]->-1}
onshell5={sp[P1,P2]->s/4/MW^2 (beta^2 +1),sp[P1,Q1]->s/4/MW^2 (beta^2 -cteta),sp[P1,Q2]->s/4/MW^2 (beta^2 +cteta), sp[P2,Q1]->s/4/MW^2 (beta^2 +cteta),sp[P2,Q2]->s/4/MW^2 (beta^2 -cteta), sp[Q1,Q2]->s/4/MW^2 (beta^2 +1)}
onshell=Flatten[{onshell1,onshell2,onshell3,onshell4,onshell5}]


mandel={sp[p1,p2]-> (s -2 MW^2)/2,sp[q1,q2]-> (s -2 MW^2)/2, sp[p1,q1]-> (2 MW^2-t)/2,sp[p2,q2]-> (2 MW^2-t)/2, sp[p1,q2]-> (s+t-2 MW^2)/2,sp[p2,q1]-> (s+t-2 MW^2)/2}
highenergy={s->4 MW^2 x,t->2 MW^2 (1 - x (1- (1 -1/x) cteta)),beta->Sqrt[1-1/x]}
weakangle={sw->Sqrt[1 -cw^2]}
ct = Part[Solve[t == 2 MW^2 - s/2 (1 - beta^2 cteta), cteta], 1]

(* Functions *)
TakeLimit=Function[{exp}, aux1=exp /. highenergy; aux2= aux1 /. x->1/z; \
	   aux3=Normal[Series[aux2,{z,0,0}]]; aux4=Expand[aux3 /. z->1/x]]

CoefA=Function[{exp},Coefficient[TakeLimit[exp],x,2]]
CoefB=Function[{exp},Coefficient[TakeLimit[exp],x,1]]
CoefC=Function[{exp},Coefficient[TakeLimit[exp],x,0]]

WriteMandel=Function[{exp}, TmpAux1 = exp /. onshell; TmpAux2 = Simplify[TmpAux1] ; \
TmpAux3 = Simplify[TmpAux2 /. ct];   TmpAux4=Simplify[TmpAux3 /. beta->Sqrt[1-4 MW^2/s]]]

WriteTeX=Function[{exp}, TmpAux1 = WriteMandel[exp]; TmpAux2=TeXForm[TmpAux1]]



(* s-channel gamma *)
ampsG= sw^2 PVL[P1,a] PVL[P2,b] PVL[Q1,c] PVL[Q2,d] mt[mu,nu] \
       V[a,b,mu,p1,p2,-p1-p2] V[d,c,nu,-q2,-q1,p1+p2] /s

MsG = Simplify[Simplify[Contract[ampsG] /. onshell] /. mandel]
MsG = Simplify[MsG /. weakangle]
MsGA:=CoefA[MsG]
MsGB:=CoefB[MsG]
MsGC:=CoefC[MsG]

(* s-channel Z *)
ampsZ= cw^2 PVL[P1,a] PVL[P2,b] PVL[Q1,c] PVL[Q2,d] \
       V[a,b,mu,p1,p2,-p1-p2] V[d,c,nu,-q2,-q1,p1+p2] \
       (mt[mu,nu] - fv[p1+p2,mu] fv[p1+p2,nu]/MW^2/cw^2)/(s -MW^2/cw^2)

MsZ = Simplify[Simplify[Contract[ampsZ] /. onshell] /. mandel]
MsZ = Simplify[MsZ /. weakangle]
MsZA:=CoefA[MsZ]
MsZB:=CoefB[MsZ]
MsZC:=CoefC[MsZ]


(* 4 Boson Interaction *)
amp4V = PVL[P1,a] PVL[P2,b] PVL[Q1,c] PVL[Q2,d] \
        (2 mt[a,d] mt[b,c] - mt[a,b] mt[d,c] - mt[a,c] mt[d,b])

M4V = Simplify[Simplify[Contract[amp4V] /. onshell] /. mandel]
M4V = Simplify[M4V /. weakangle]
M4VA:=CoefA[M4V]
M4VB:=CoefB[M4V]
M4VC:=CoefC[M4V]

(* t- channel gamma *)
amptG =  sw^2 PVL[P1,a] PVL[P2,b] PVL[Q1,c] PVL[Q2,d] mt[mu,nu] \
       V[a,c,mu,p1,-q1,q1-p1] V[d,b,nu,-q2,p2,q2-p2] / t

MtG = Simplify[Simplify[Contract[amptG] /. onshell] /. mandel]
MtG = Simplify[MtG /. weakangle]
MtGA:=CoefA[MtG]
MtGB:=CoefB[MtG]
MtGC:=CoefC[MtG]

(* t-channel Z *)
amptZ =  cw^2 PVL[P1,a] PVL[P2,b] PVL[Q1,c] PVL[Q2,d] \
       V[a,c,mu,p1,-q1,q1-p1] V[d,b,nu,-q2,p2,q2-p2] \
       (mt[mu,nu] - fv[p1-q1,mu] fv[p1-q1,nu]/MW^2/cw^2)/(t -MW^2/cw^2)

MtZ = Simplify[Simplify[Contract[amptZ] /. onshell] /. mandel]
MtZ = Simplify[MtZ /. weakangle]
MtZA:=CoefA[MtZ]
MtZB:=CoefB[MtZ]
MtZC:=CoefC[MtZ]
         
(* s-channel Higgs *)
ampsH = - MW^2 PVL[P1,a] PVL[P2,b] PVL[Q1,c] PVL[Q2,d] \
          mt[a,b] mt[c,d] /(s-MH^2)

MsH = Simplify[Simplify[Contract[ampsH] /. onshell] /. mandel]
MsH = Simplify[MsH /. weakangle]
MsHA:=CoefA[MsH]
MsHB:=CoefB[MsH]
MsHC:=CoefC[MsH]

(* t-channel Higgs *)
amptH = - MW^2 PVL[P1,a] PVL[P2,b] PVL[Q1,c] PVL[Q2,d] \
          mt[a,c] mt[b,d] /(t-MH^2)

MtH = Simplify[Simplify[Contract[amptH] /. onshell] /. mandel]
MtH = Simplify[MtH /. weakangle]
MtHA:=CoefA[MtH]
MtHB:=CoefB[MtH]
MtHC:=CoefC[MtH]


testA:=Simplify[MsGA+MsZA+MtGA+MtZA+M4VA+MsHA+MtHA /. weakangle]
testB:=Simplify[MsGB+MsZB+MtGB+MtZB+M4VB+MsHB+MtHB /. weakangle]


(* Output Fortran *)
xs=0

If[xs==1,
   stmp=OpenWrite["MsG.f",FormatType -> FortranForm];
   Write[stmp,g^2 MsG /. g^2->8 MW^2 GF/Sqrt[2]];
   Close[stmp];
   stmp=OpenWrite["MsZ.f",FormatType -> FortranForm];
   Write[stmp,g^2 MsZ /. g^2->8 MW^2 GF/Sqrt[2]];
   Close[stmp];
   stmp=OpenWrite["MtG.f",FormatType -> FortranForm];
   Write[stmp,g^2 MtG /. g^2->8 MW^2 GF/Sqrt[2]];
   Close[stmp];
   stmp=OpenWrite["MtZ.f",FormatType -> FortranForm];
   Write[stmp,g^2 MtZ /. g^2->8 MW^2 GF/Sqrt[2]];
   Close[stmp];
   stmp=OpenWrite["M4V.f",FormatType -> FortranForm];
   Write[stmp,g^2 M4V /. g^2->8 MW^2 GF/Sqrt[2]];
   Close[stmp];
   stmp=OpenWrite["MsH.f",FormatType -> FortranForm];
   Write[stmp,g^2 MsH /. g^2->8 MW^2 GF/Sqrt[2]];
   Close[stmp];
   stmp=OpenWrite["MtH.f",FormatType -> FortranForm];
   Write[stmp,g^2 MtH /. g^2->8 MW^2 GF/Sqrt[2]];
   Close[stmp];
   ,Print[" "]];




