(* Remove["Global`*"] *)

xs=0

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_] := Pair[Momentum[p], Momentum[q]]
li[mu_] := LorentzIndex[mu]
prop[p_, m_] := ds[p] + m

pv[p_,mu_]:= PolarizationVector[p,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] 


(* Couplings *)

PL=dm[7]
PR=dm[6]

gR=sw^2
gL=-1/2+sw^2


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

vp1={Sqrt[s]/2,0,0,betae Sqrt[s]/2}
vp2={Sqrt[s]/2,0,0,-betae Sqrt[s]/2}
vq1={Sqrt[s]/2,Sqrt[s]/2 beta steta,0,Sqrt[s]/2 beta cteta}
vq2={Sqrt[s]/2,-Sqrt[s]/2 beta steta,0,-Sqrt[s]/2 beta cteta}
vQ1L={Sqrt[s]/2/MW beta,Sqrt[s]/2/MW  steta,0,Sqrt[s]/2/MW cteta}
vQ2L={Sqrt[s]/2/MW beta,-Sqrt[s]/2/MW  steta,0,-Sqrt[s]/2/MW cteta}
vQ1P={0,cteta/Sqrt[2], I/Sqrt[2] , -steta/Sqrt[2]}
vQ1M={0,cteta/Sqrt[2], -I/Sqrt[2] , -steta/Sqrt[2]}
vQ2P={0,-cteta/Sqrt[2], I/Sqrt[2] , steta/Sqrt[2]}
vQ2M={0,-cteta/Sqrt[2], -I/Sqrt[2] , steta/Sqrt[2]}
vQ1cL=vQ1L
vQ2cL=vQ2L
vQ1cP={0,cteta/Sqrt[2], -I/Sqrt[2] , -steta/Sqrt[2]}
vQ1cM={0,cteta/Sqrt[2], I/Sqrt[2] , -steta/Sqrt[2]}
vQ2cP={0,-cteta/Sqrt[2],-I/Sqrt[2] , steta/Sqrt[2]}
vQ2cM={0,-cteta/Sqrt[2], I/Sqrt[2] , steta/Sqrt[2]}


GetEps=Function[{a,b,c,d},mat={a,b,c,d}; Det[mat]]
GetDot=Function[{a,b},a[[1]]*b[[1]]-a[[2]]*b[[2]]-a[[3]]*b[[3]]-a[[4]]*b[[4]]]

vlist={p1,p2,q1,q2,Q1L,Q2L,Q1P,Q1M,Q2P,Q2M,Q1cL,Q2cL,Q1cP,Q1cM,Q2cP,Q2cM}
fvlist={vp1,vp2,vq1,vq2,vQ1L,vQ2L,vQ1P,vQ1M,vQ2P,vQ2M,vQ1cL,vQ2cL,vQ1cP,vQ1cM,vQ2cP,vQ2cM}
m2list={me^2,me^2,MW^2,MW^2,-1,-1,0,0,0,0,-1,-1,0,0,0,0}
onshell=Flatten[Table[If[i==j,Pair[Momentum[vlist[[i]]],Momentum[vlist[[j]]]]->m2list[[i]],Pair[Momentum[vlist[[i]]],Momentum[vlist[[j]]]]-> GetDot[fvlist[[i]],fvlist[[j]]]],{i,1,16},{j,i,16}]];
onshell = Simplify[onshell /. {steta^2->1-cteta^2,beta^2->1-4 MW^2/s}];


Q1=Polarization[q1];
Q2=Polarization[q2];
Q1c=Momentum[Polarization[q1, -I]];
Q2c=Momentum[Polarization[q2, -I]];

VList={p1,p2,q1,q2,Q1,Q1c,Q2,Q2c}
VListLL={p1,p2,q1,q2,Q1L,Q1cL,Q2L,Q2cL}
VListLP={p1,p2,q1,q2,Q1L,Q1cL,Q2P,Q2cP}
VListLM={p1,p2,q1,q2,Q1L,Q1cL,Q2M,Q2cM}
VListPL={p1,p2,q1,q2,Q1P,Q1cP,Q2L,Q2cL}
VListML={p1,p2,q1,q2,Q1M,Q1cM,Q2L,Q2cL}
VListPP={p1,p2,q1,q2,Q1P,Q1cP,Q2P,Q2cP}
VListPM={p1,p2,q1,q2,Q1P,Q1cP,Q2M,Q2cM}
VListMP={p1,p2,q1,q2,Q1M,Q1cM,Q2P,Q2cP}
VListMM={p1,p2,q1,q2,Q1M,Q1cM,Q2M,Q2cM}

M2List={me^2,me^2,MW^2,MW^2}
ms1=Flatten[Table[sp[VList[[i]],VList[[i]]]->M2List[[i]],{i,1,4}]]
ms2=Flatten[Table[sp[vlist[[i]],vlist[[j]]]-> GetDot[fvlist[[i]],fvlist[[j]]],{i,1,4},{j,i+1,4}]]
mshell=ms1


PolLL = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListLL[[i]],VListLL[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolLL = DeleteCases[DeleteCases[PolLL,0->0],-1->-1];
PolLP = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListLP[[i]],VListLP[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolLP = DeleteCases[DeleteCases[PolLP,0->0],-1->-1];
PolLM = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListLM[[i]],VListLM[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolLM = DeleteCases[DeleteCases[PolLM,0->0],-1->-1];
PolPL = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListPL[[i]],VListPL[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolPL = DeleteCases[DeleteCases[PolPL,0->0],-1->-1];
PolML = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListML[[i]],VListML[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolML = DeleteCases[DeleteCases[PolML,0->0],-1->-1];
PolPP = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListPP[[i]],VListPP[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolPP = DeleteCases[DeleteCases[PolPP,0->0],-1->-1];
PolPM = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListPM[[i]],VListPM[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolPM = DeleteCases[DeleteCases[PolPM,0->0],-1->-1];
PolMP = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListMP[[i]],VListMP[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolMP = DeleteCases[DeleteCases[PolMP,0->0],-1->-1];
PolMM = Flatten[{mshell,Table[sp[VList[[i]],VList[[j]]]-> (sp[VListMM[[i]],VListMM[[j]]] /. onshell),{i,1,8},{j,i+1,8}]}];
PolMM = DeleteCases[DeleteCases[PolMM,0->0],-1->-1];

fvListLL={vp1,vp2,vq1,vq2,vQ1L,vQ1cL,vQ2L,vQ2cL}
fvListLP={vp1,vp2,vq1,vq2,vQ1L,vQ1cL,vQ2P,vQ2cP}
fvListLM={vp1,vp2,vq1,vq2,vQ1L,vQ1cL,vQ2M,vQ2cM}
fvListPL={vp1,vp2,vq1,vq2,vQ1P,vQ1cP,vQ2L,vQ2cL}
fvListML={vp1,vp2,vq1,vq2,vQ1M,vQ1cM,vQ2L,vQ2cL}
fvListPP={vp1,vp2,vq1,vq2,vQ1P,vQ1cP,vQ2P,vQ2cP}
fvListPM={vp1,vp2,vq1,vq2,vQ1P,vQ1cP,vQ2M,vQ2cM}
fvListMP={vp1,vp2,vq1,vq2,vQ1M,vQ1cM,vQ2P,vQ2cP}
fvListMM={vp1,vp2,vq1,vq2,vQ1M,vQ1cM,vQ2M,vQ2cM}

subepsLL = DeleteCases[Flatten[Table[If[i != j && i != k && i != m && j != k\
	&& j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
      Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListLL[[i]], \
    fvListLL[[j]], fvListLL[[k]], fvListLL[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8}, {k, 1, 8}, {m, 1, 8}]],_Integer];

subepsLP = DeleteCases[Flatten[Table[If[i != j && i != k && i != m && j != k \
         && j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
    Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListLP[[i]], \
    fvListLP[[j]], fvListLP[[k]], fvListLP[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8},{k, 1, 8}, {m, 1, 8}]],_Integer];

subepsLM = DeleteCases[Flatten[Table[If[i != j && i != k && i != m && j != k \
         && j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
      Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListLM[[i]], \
    fvListLM[[j]], fvListLM[[k]], fvListLM[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8},{k, 1, 8}, {m, 1, 8}]],_Integer];

subepsPL = DeleteCases[Flatten[Table[If[i != j && i != k && i != m && j != k \
         && j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
       Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListPL[[i]], \
    fvListPL[[j]], fvListPL[[k]], fvListPL[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8},{k, 1, 8}, {m, 1, 8}]],_Integer];

subepsML = DeleteCases[Flatten[Table[If[i != j && i != k && i != m && j != k \
         && j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
      Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListML[[i]], \
    fvListML[[j]], fvListML[[k]], fvListML[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8},{k, 1, 8}, {m, 1, 8}]],_Integer];

subepsPP = DeleteCases[Flatten[Table[If[i != j && i != k && i != m && j != k && j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
    Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListPP[[i]], \
    fvListPP[[j]], fvListPP[[k]], fvListPP[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8},{k, 1, 8}, {m, 1, 8}]],_Integer];
subepsPM =DeleteCases[ Flatten[Table[If[i != j && i != k && i != m && j != k \
         && j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
    Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListPM[[i]], \
    fvListPM[[j]], fvListPM[[k]], fvListPM[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8}, {k, 1, 8}, {m, 1, 8}]],_Integer];

subepsMP = DeleteCases[Flatten[Table[If[i != j && i != k && i != m && j != k \
         && j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
    Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListMP[[i]], \
    fvListMP[[j]], fvListMP[[k]], fvListMP[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8},{k, 1, 8}, {m, 1, 8}]],_Integer];

subepsMM = DeleteCases[Flatten[Table[If[i != j && i != k && i != m && j != k \
         && j != m && k != m,Eps[Momentum[VList[[i]]], Momentum[VList[[j]]], \
    Momentum[VList[[k]]], Momentum[VList[[m]]]] -> (GetEps[fvListMM[[i]], \
    fvListMM[[j]], fvListMM[[k]], fvListMM[[m]]] /. onshell),1], {i, 1, 8}, \
    {j, 1, 8}, {k, 1, 8}, {m, 1, 8}]],_Integer];

Pol={PolLL,PolLM,PolLP,PolML,PolPL,PolPP,PolPM,PolMP,PolMM}
subeps={subepsLL,subepsLM,subepsLP,subepsML,subepsPL,subepsPP,subepsPM,subepsMP,subepsMM}
simpbetas={beta -> Sqrt[1 - 4 MW^2/s], betae -> Sqrt[1 - 4 me^2/s]}
weakangle={sw->Sqrt[1 -cw^2]}



he1={s->4 MW^2 x}
he2={t-> MW^2 + me^2 - 2 MW^2 x (1- Sqrt[(1 -1/x)] Sqrt[1-me^2/MW^2/x] cteta)}
he3={beta->Sqrt[1-1/x]}
he4={betae->Sqrt[1-me^2/x/MW^2]}

highenergy=Flatten[{he1,he2,he3,he4}]

weakangle={sw->Sqrt[1 -cw^2]}
simpbetas={beta -> Sqrt[1 - 4 MW^2/s], betae -> Sqrt[1 - 4 me^2/s]}

ds1LL={DiracGamma[Momentum[Polarization[q1,I]]]->1/MW DiracGamma[Momentum[q1]]};
ds2LL={DiracGamma[Momentum[Polarization[q2,I]]]->1/MW DiracGamma[Momentum[q2]]};
diracsubLL=Flatten[{ds1LL,ds2LL}];

diracsubLT={DiracGamma[Momentum[Polarization[q1, I]]]->1/MW DiracGamma[Momentum[q1]]};
diracsubTL={DiracGamma[Momentum[Polarization[q2, I]]]->1/MW DiracGamma[Momentum[q2]]};
diracsubTT={};

diracsub1=diracsubLL
diracsub2={DiracGamma[Momentum[q2]]->DiracGamma[Momentum[p1+p2-q1]]}

PolType=Pol[[1]]

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

ct = Part[Solve[t == me^2+ MW^2 - s/2 (1 - beta*betae cteta), cteta], 1]
st2 ={steta^2->1-cteta^2}
beta2={beta^2->1-4 MW^2/s}

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

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



(* Amplitudes *)
LineG= Spinor[-p2,me] . dm[a] . (PL + PR) . Spinor[p1,me] 
ampsG= sw^2 LineG pv[q1,mu] pv[q2,nu] V[nu,mu,b,-q2,-q1,p1+p2] (-mt[a,b])/s
MsG = DiracSimplify[Contract[ampsG] /. onshell] 

LineZ = Spinor[-p2,me] . dm[a] . (gL PL + gR PR).  Spinor[p1,me]
PropZ = (-mt[a,b] +  cw^2 fv[p1+p2,a] fv[p1+p2,b]/MW^2)/(s-MW^2/cw^2)
ampsZ= - pv[q1,mu] pv[q2,nu] LineZ V[nu,mu,b,-q2,-q1,p1+p2] PropZ
MsZ = DiracSimplify[Contract[ampsZ] /. onshell] 


LineNeu=Spinor[-p2,me] . dm[nu] . PL . prop[p1-q1,0] . dm[mu] . PL . Spinor[p1,me]
ampNeu= - 1/2 LineNeu pv[q1,mu] pv[q2,nu] /t
MtNeu= DiracSimplify[Contract[ampNeu] /. onshell] 


LineH = Spinor[-p2,me] . (PL+PR). Spinor[p1,me]
ampH =  me/2 LineH pv[q1,mu] pv[q2,nu] mt[mu,nu]/(s-mH^2)
MsH= DiracSimplify[Contract[ampH] /. onshell]


(* High Energy Limit *)


MtNeuaux= DiracSimplify[MtNeu /. diracsub1]
MtNeuaux= DiracSimplify[MtNeuaux /. diracsub2]
MtNeuaux= DiracSimplify[MtNeuaux /. PolType]
MtNeuaux= DiracSimplify[MtNeuaux /. diracsub1]
MtNeuHEL = TakeLimit[MtNeuaux,0]

MsGaux= DiracSimplify[MsG /. diracsub1]
MsGaux= DiracSimplify[MsGaux /. diracsub2]
MsGaux= DiracSimplify[MsGaux /. PolType]
MsGaux= DiracSimplify[MsGaux /. diracsub1]
MsGHEL = TakeLimit[MsGaux,0]

MsZaux= DiracSimplify[MsZ /. diracsub1]
MsZaux= DiracSimplify[MsZaux /. diracsub2]
MsZaux= DiracSimplify[MsZaux /. PolType]
MsZaux= DiracSimplify[MsZaux /. diracsub1]
MsZHEL = TakeLimit[MsZaux,0]


MsHaux= DiracSimplify[MsH /. diracsub1]
MsHaux= DiracSimplify[MsHaux /. diracsub2]
MsHaux= DiracSimplify[MsHaux /. PolType]
MsHHEL = TakeLimit[MsHaux,0]


MFinalHEL:= Simplify[MtNeuHEL+MsGHEL+MsZHEL+MsHHEL]

  (*
(* cross section *)
GammaNu =  - 1/(2 t) dm[nu] . PL . prop[p1-q1,0] . dm[mu] . PL 
GammaZ1 = - dm[b] . (gL PL + gR PR) V[nu,mu,b,-q2,-q1,p1+p2] 
GammaZ2 = (-me) (gL PL + gR PR) V[nu,mu,b,-q2,-q1,p1+p2] cw^2 fv[p1+p2,b]/MW^2 
GammaZ3 = (me) (gL PR + gR PL) V[nu,mu,b,-q2,-q1,p1+p2]  cw^2 fv[p1+p2,b]/MW^2 
GammaZ  = - 1/(s-MW^2/cw^2) (GammaZ1+GammaZ2+GammaZ3)
GammaG  =  sw^2/s  dm[a] . (PL + PR)  V[nu,mu,b,-q2,-q1,p1+p2] (-mt[a,b])
GammaH  =  me/2 (PL+PR) mt[mu,nu]/(s-mH^2)
GammaBarNu =  - 1/(2 t) dm[mup] . PL . prop[p1-q1,0] . dm[nup] . PL 
GammaBarZ1 = - dm[bp] . (gL PL + gR PR) V[nup,mup,bp,-q2,-q1,p1+p2] 
GammaBarZ2 = (-me) (gL PL + gR PR) V[nup,mup,bp,-q2,-q1,p1+p2] cw^2 fv[p1+p2,bp]/MW^2 
GammaBarZ3 = (me) (gL PR + gR PL) V[nup,mup,bp,-q2,-q1,p1+p2]  cw^2 fv[p1+p2,bp]/MW^2 
GammaBarZ  = - 1/(s-MW^2/cw^2) (GammaBarZ1+GammaBarZ2+GammaBarZ3)
GammaBarG  =  sw^2/s  dm[ap] . (PL + PR)  V[nup,mup,bp,-q2,-q1,p1+p2] (-mt[ap,bp])
GammaBarH  =  me/2 (PL+PR) mt[mup,nup]/(s-mH^2)



(* Neutrino + gamma + Z cross section *)
LineNuNu =  prop[p2,-me] . GammaNu . prop[p1,me] . GammaBarNu
MsquareNuNu=Contract[Tr[LineNuNu] pv[q1,mu] pv[q2,nu] Conjugate[pv[q1,mup]] Conjugate[pv[q2,nup]]] /. mshell;

MsqNuNu = Table[(MsquareNuNu /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqNuNu = (MsqNuNu /. st2) /. beta2 //Simplify
Clear[MsquareNuNu]

LineNuG =  prop[p2,-me] . GammaNu . prop[p1,me] . GammaBarG 
MsquareNuG = Contract[Tr[LineNuG] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;

MsqNuG = Table[(MsquareNuG /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqNuG = (MsqNuG /. st2) /. beta2 //Simplify
Clear[MsquareNuG]


LineGNu =  prop[p2,-me] . GammaG . prop[p1,me] . GammaBarNu
MsquareGNu := Contract[Tr[LineGNu] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqGNuAux:=Table[(MsquareGNu /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}];
TestGNuNuG=Function[{},Msq=MsquareGNu; MsqGNuAux=Table[(Msq /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}];  Clear[Msq,MsquareGNu]; MsqGNu = (MsqGNuAux /. st2) /. beta2 //Simplify; Table[MsqNuG[[i]]-MsqGNu[[i]],{i,1,9}]]


LineNuZ =  prop[p2,-me] . GammaNu . prop[p1,me] . GammaBarZ
MsquareNuZ = Contract[Tr[LineNuZ] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;

MsqNuZ = Table[(MsquareNuZ /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqNuZ = (MsqNuZ /. st2) /. beta2 //Simplify
Clear[MsquareNuZ]


LineZNu =  prop[p2,-me] . GammaZ . prop[p1,me] . GammaBarNu
MsquareZNu := Contract[Tr[LineZNu] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqZNuAux:=Table[(MsquareZNu /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}];
TestZNuNuZ=Function[{},Msq=MsquareZNu; MsqZNuAux=Table[(Msq /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}];  Clear[Msq,MsquareZNu]; MsqZNu = (MsqZNuAux /. st2) /. beta2 //Simplify; Table[MsqNuZ[[i]]-MsqZNu[[i]],{i,1,9}]]


LineGG =  prop[p2,-me] . GammaG . prop[p1,me] . GammaBarG
MsquareGG =  Contract[Tr[LineGG] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqGG = Table[(MsquareGG /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqGG = (MsqGG /. st2) /. beta2 //Simplify
Clear[MsquareGG]


LineZZ =  prop[p2,-me] . GammaZ . prop[p1,me] . GammaBarZ
MsquareZZ =  Contract[Tr[LineZZ] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqZZ = Table[(MsquareZZ /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqZZ = (MsqZZ /. st2) /. beta2 //Simplify
Clear[MsquareZZ]

 
LineGZ =  prop[p2,-me] . GammaG . prop[p1,me] . GammaBarZ
MsquareGZ =  Contract[Tr[LineGZ] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqGZ = Table[(MsquareGZ /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqGZ = (MsqGZ /. st2) /. beta2 //Simplify
Clear[MsquareGZ]

LineZG =  prop[p2,-me] . GammaZ . prop[p1,me] . GammaBarG
MsquareZG:=  Contract[Tr[LineZG] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqZGAux:=Table[(MsquareZG /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}];
TestZGGZ=Function[{},Msq=MsquareZG; Clear[MsquareZG]; MsqZGAux=Table[(Msq /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]; Clear[Msq,MsquareZG]; MsqZG = (MsqZGAux /. st2) /. beta2 //Simplify; Table[MsqGZ[[i]]-MsqZG[[i]],{i,1,9}]]



LineHH =  prop[p2,-me] . GammaH . prop[p1,me] . GammaBarH
MsquareHH =  Contract[Tr[LineHH] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqHH = Table[(MsquareHH /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqHH = (MsqHH /. st2) /. beta2 //Simplify
Clear[MsquareHH]

LineNuH =  prop[p2,-me] . GammaNu . prop[p1,me] . GammaBarH
MsquareNuH = Contract[Tr[LineNuH] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqNuH = Table[(MsquareNuH /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqNuH = (MsqNuH /. st2) /. beta2 //Simplify
Clear[MsquareNuH]


LineHNu =  prop[p2,-me] . GammaH . prop[p1,me] . GammaBarNu
MsquareHNu:= Contract[Tr[LineHNu] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqHNuAux:=Table[(MsquareHNu /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}];
TestHNuNuH=Function[{},Msq=MsquareHNu; Clear[MsquareHNu]; MsqHNuAux=Table[(Msq /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]; Clear[Msq,MsquareHNu]; MsqHNu = (MsqHNuAux /. st2) /. beta2 //Simplify; Table[MsqNuH[[i]]-MsqHNu[[i]],{i,1,9}]]

LineGH =  prop[p2,-me] . GammaG . prop[p1,me] . GammaBarH
MsquareGH = Contract[Tr[LineGH] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqGH = Table[(MsquareGH /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqGH = (MsqGH /. st2) /. beta2 //Simplify
Clear[MsquareGH]


LineHG =  prop[p2,-me] . GammaH . prop[p1,me] . GammaBarG
MsquareHG:= Contract[Tr[LineHG] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqHGAux:=Table[(MsquareHG /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}];
TestHGGH=Function[{},Msq=MsquareHG; Clear[MsquareHG]; MsqHGAux=Table[(Msq /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]; Clear[Msq,MsquareHG]; MsqHG = (MsqHGAux /. st2) /. beta2 //Simplify; Table[MsqGH[[i]]-MsqHG[[i]],{i,1,9}]]

LineZH =  prop[p2,-me] . GammaZ . prop[p1,me] . GammaBarH
MsquareZH = Contract[Tr[LineZH] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqZH = Table[(MsquareZH /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]
MsqZH = (MsqZH /. st2) /. beta2 //Simplify
Clear[MsquareZH]


LineHZ =  prop[p2,-me] . GammaH . prop[p1,me] . GammaBarZ
MsquareHZ:= Contract[Tr[LineHZ] pv[q1,mu] pv[q2,nu]  Conjugate[pv[q1,mup]]  Conjugate[pv[q2,nup]]] /. mshell;
MsqHZAux:=Table[(MsquareHZ /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}];
TestHZZH=Function[{},Msq=MsquareHZ; Clear[MsquareHZ]; MsqHZAux=Table[(Msq /. subeps[[i]]) /. Pol[[i]] // Simplify,{i,1,9}]; Clear[Msq,MsquareHZ]; MsqHZ = (MsqHZAux /. st2) /. beta2 //Simplify; Table[MsqZH[[i]]-MsqHZ[[i]],{i,1,9}]]


								     
stmp=OpenWrite["status.1"];
Write[stmp,"Ponto 1"];
Close[stmp]


MsqNuLL = Simplify[MsqNuNu[[1]]*g^4/4 /. ct] /. g^4->32 MW^4 GF^2
MsqNuLL = Simplify[MsqNuLL /. simpbetas];

MsqHLL = Simplify[MsqHH[[1]]*g^4/4 /. ct ] /. g^4->32 MW^4 GF^2 ;
MsqHLL= Simplify[MsqHLL /. simpbetas];

MsqNuGZLL=MsqNuNu[[1]] + MsqGG[[1]] + MsqZZ[[1]] + \
          2*MsqNuG[[1]] + 2*MsqNuZ[[1]] + 2*MsqGZ[[1]];
MsqNuGZLL = Simplify[MsqNuGZLL*g^4/4 /. ct] /. g^4->32 MW^4 GF^2 ;
MsqNuGZLL =Simplify[MsqNuGZLL /. simpbetas];

MsqHiggsLL = MsqHH[[1]] + 2*MsqNuH[[1]] + 2*MsqGH[[1]] + 2*MsqZH[[1]] ;
MsqHiggsLL = Simplify[MsqHiggsLL*g^4/4 /. ct] /. g^4->32 MW^4 GF^2 ;
MsqHiggsLL = Simplify[MsqHiggsLL /. simpbetas];

MsqLL = MsqNuGZLL + MsqHiggsLL ;

stmp=OpenWrite["status.2"];
Write[stmp,"Ponto 2"];
Close[stmp]

MsqLL=Simplify[MsqLL /. steta->Sqrt[1-cteta^2]];
MsqLL=Simplify[MsqLL /. ct];
MsqLL=Simplify[MsqLL /. simpbetas];

If[xs==1,
   SetOptions[$Output, PageWidth -> 60];
   stmp=OpenWrite["ansNu.f",FormatType -> FortranForm];
   Write[stmp,"MsqNu      =      ", MsqNuLL];
   Close[stmp];
   stmp=OpenWrite["ansNuGZ.f",FormatType -> FortranForm];
   Write[stmp,"MsqNuGZ      =      ", MsqNuGZLL];
   Close[stmp];
   stmp=OpenWrite["ansLL.f",FormatType -> FortranForm];
   Write[stmp,"MsqLL      =      ", MsqLL];
   Close[stmp];
   stmp=OpenWrite["ansH.f",FormatType -> FortranForm];
   Write[stmp,"MsqH      =      ", MsqHLL];
   Close[stmp], Print[" "]];


stmp=OpenWrite["status.3"];
Write[stmp,"Ponto 3"];
Close[stmp]

(* For the Transverse polarization put me=0, betae->1 *)
(* Recall Pol={PolLL,PolLM,PolLP,PolML,PolPL,PolPP,PolPM,PolMP,PolMM} *)
(*               1     2     3     4     5     6     7     8     9    *)



MsqTL=Sum[MsqNuNu[[i]] + MsqGG[[i]] + MsqZZ[[i]] + \
          2*MsqNuG[[i]] + 2*MsqNuZ[[i]] + 2*MsqGZ[[i]] +\
          MsqHH[[i]] + 2*MsqNuH[[i]] + 2*MsqGH[[i]] + 2*MsqZH[[i]] \
          /. {me->0,betae->1} //Simplify,{i,2,5}] ;

MsqTL=Simplify[MsqTL /. steta->Sqrt[1-cteta^2]];
MsqTL=Simplify[MsqTL /. ct];
MsqTL=Simplify[MsqTL /. {me->0,betae->1}];
MsqTL=Simplify[MsqTL*g^4/4 /. weakangle] /. g^4->32 MW^4 GF^2;
MsqTL=Simplify[MsqTL /. beta2];

MsqTT=Sum[MsqNuNu[[i]] + MsqGG[[i]] + MsqZZ[[i]] + \
          2*MsqNuG[[i]] + 2*MsqNuZ[[i]] + 2*MsqGZ[[i]] +\
          MsqHH[[i]] + 2*MsqNuH[[i]] + 2*MsqGH[[i]] + 2*MsqZH[[i]] \
          /. {me->0,betae->1} //Simplify,{i,6,9}] ;

MsqTT=Simplify[MsqTT /. steta->Sqrt[1-cteta^2]];
MsqTT=Simplify[MsqTT /. ct];
MsqTT=Simplify[MsqTT /. {me->0,betae->1}];
MsqTT=Simplify[MsqTT*g^4/4 /. weakangle] /. g^4->32 MW^4 GF^2
MsqTT=Simplify[MsqTT /. beta2];

MsqTotal = Simplify[MsqLL + MsqTL + MsqTT]

(* Write Fortran Code *)

xs=1

stmp=OpenWrite["status.4"];
Write[stmp,"Ponto 4"];
Close[stmp]


If[xs==1,
   SetOptions[$Output, PageWidth -> 70];
   stmp=OpenWrite["ansTL.f",FormatType -> FortranForm];
   Write[stmp,"MsqTL = ", MsqTL];
   Close[stmp];
   stmp=OpenWrite["ansTT.f",FormatType -> FortranForm];
   Write[stmp,"MsqTT = ", MsqTT];
   Close[stmp];
   stmp=OpenWrite["ansTotal.f",FormatType -> FortranForm];
   Write[stmp,"MsqTotal = ", MsqTotal];
   Close[stmp], Print[" "]];


	  *)
