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


(* 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+ *)

os1={sp[p1,p1]->me^2,sp[p2,p2]->me^2,sp[q1,q1]->MW^2,sp[q2,q2]->MW^2}
os2={sp[Q1,Q1]->-1,sp[Q2,Q2]->-1}
os3={sp[Q1,q1]->0,sp[Q2,q2]->0}
os4={sp[p1,p2]-> (s -2 me^2)/2}
os5={sp[q1,q2]-> (s -2 MW^2)/2}
os6={sp[p1,q1]-> (MW^2 +me^2-t)/2}
os7={sp[p2,q2]-> ( MW^2+me^2-t)/2}
os8={sp[p1,q2]-> (s+t- MW^2-me^2)/2}
os9={sp[p2,q1]-> (s+t- MW^2-me^2)/2}
onshell=Flatten[{os1,os2,os3,os4,os5,os6,os7,os8,os9}]


os1L={sp[Q1,p1]->s/4/MW (beta - betae cteta)}
os2L={sp[Q1,p2]->s/4/MW (beta + betae cteta)}
os3L={sp[Q1,q2]->s/2/MW beta}
os4L={sp[Q2,p1]->s/4/MW (beta +betae cteta)}
os5L={sp[Q2,p2]->s/4/MW (beta - betae cteta)}
os6L={sp[Q2,q1]->s/2/MW beta}

os1T1={sp[Q1,p1]->0}
os2T1={sp[Q1,p2]->0}
os3T1={sp[Q1,q2]->0}
os4T1={sp[Q2,p1]->0}
os5T1={sp[Q2,p2]->0}
os6T1={sp[Q2,q1]->0}

os1T2={sp[Q1,p1]->Sqrt[s]/2 betae (+steta)}
os2T2={sp[Q1,p2]->Sqrt[s]/2 betae (-steta)}
os3T2={sp[Q1,q2]->0}
os4T2={sp[Q2,p1]->Sqrt[s]/2 betae (-steta)}
os5T2={sp[Q2,p2]->Sqrt[s]/2 betae (+steta)}
os6T2={sp[Q2,q1]->0}

osLL={sp[Q1,Q2]->s/4/MW^2 (beta^2 +1)}
osT1T1={sp[Q1,Q2]->-1}
osT1T2={sp[Q1,Q2]->0}
osT2T1={sp[Q1,Q2]->0}
osT2T2={sp[Q1,Q2]->1}
osT1L={sp[Q1,Q2]->0}
osT2L={sp[Q1,Q2]->0}
osLT1={sp[Q1,Q2]->0}
osLT2={sp[Q1,Q2]->0}


PolLL=Flatten[{os1L,os2L,os3L,os4L,os5L,os6L,osLL}]
PolT1L=Flatten[{os1T1,os2T1,os3T1,os4L,os5L,os6L,osT1L}]
PolT2L=Flatten[{os1T2,os2T2,os3T2,os4L,os5L,os6L,osT2L}]
PolLT1=Flatten[{os1L,os2L,os3L,os4T1,os5T1,os6T1,osLT1}]
PolLT2=Flatten[{os1L,os2L,os3L,os4T2,os5T2,os6T2,osLT2}]
PolT1T1=Flatten[{os1T1,os2T1,os3T1,os4T1,os5T1,os6T1,osT1T1}]
PolT1T2=Flatten[{os1T1,os2T1,os3T1,os4T2,os5T2,os6T2,osT1T2}]
PolT2T1=Flatten[{os1T2,os2T2,os3T2,os4T1,os5T1,os6T1,osT2T1}]
PolT2T2=Flatten[{os1T2,os2T2,os3T2,os4T2,os5T2,os6T2,osT2T2}]

Pol=PolLL

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]}
subeps={Eps[Momentum[a_],Momentum[b_],Momentum[c_],Momentum[d_]]->0}
simpbetas={beta -> Sqrt[1 - 4 MW^2/s], betae -> Sqrt[1 - 4 me^2/s]}

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

diracsubLT={DiracGamma[Momentum[Q1]]->1/MW DiracGamma[Momentum[q1]]};
diracsubTL={DiracGamma[Momentum[Q2]]->1/MW DiracGamma[Momentum[q2]]};
diracsubTT={};

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


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

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 PVL[Q1,mu] PVL[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= - PVL[Q1,mu] PVL[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 PVL[Q1,mu] PVL[Q2,nu] /t
MtNeu= DiracSimplify[Contract[ampNeu] /. onshell] 


LineH = Spinor[-p2,me] . (PL+PR). Spinor[p1,me]
ampH =  me/2 LineH PVL[Q1,mu] PVL[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 /. onshell]
MtNeuaux= DiracSimplify[MtNeuaux /. Pol]
MtNeuHEL = TakeLimit[MtNeuaux,1]

MsGaux= DiracSimplify[MsG /. diracsub1]
MsGaux= DiracSimplify[MsGaux /. diracsub2]
MsGaux= DiracSimplify[MsGaux /. onshell]
MsGaux= DiracSimplify[MsGaux /. Pol]
MsGHEL = TakeLimit[MsGaux,1]

MsZaux= DiracSimplify[MsZ /. diracsub1]
MsZaux= DiracSimplify[MsZaux /. diracsub2]
MsZaux= DiracSimplify[MsZaux /. onshell]
MsZaux= DiracSimplify[MsZaux /. Pol]
MsZHEL = TakeLimit[MsZaux,1]


MsHaux= DiracSimplify[MsH /. diracsub1]
MsHaux= DiracSimplify[MsHaux /. diracsub2]
MsHaux= DiracSimplify[MsHaux /. onshell]
MsHaux= DiracSimplify[MsHaux /. Pol]
MsHHEL = TakeLimit[MsHaux,1]


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
MsqNuNu = Contract[Tr[LineNuNu] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineNuG =  prop[p2,-me] . GammaNu . prop[p1,me] . GammaBarG
MsqNuG = Contract[Tr[LineNuG] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineNuZ =  prop[p2,-me] . GammaNu . prop[p1,me] . GammaBarZ
MsqNuZ = Contract[Tr[LineNuZ] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineGG =  prop[p2,-me] . GammaG . prop[p1,me] . GammaBarG
MsqGG =  Contract[Tr[LineGG] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineGNu =  prop[p2,-me] . GammaG . prop[p1,me] . GammaBarNu
MsqGNu = Contract[Tr[LineGNu] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineGZ =  prop[p2,-me] . GammaG . prop[p1,me] . GammaBarZ
MsqGZ =  Contract[Tr[LineGZ] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineZG =  prop[p2,-me] . GammaZ . prop[p1,me] . GammaBarG
MsqZG =  Contract[Tr[LineZG] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineZNu =  prop[p2,-me] . GammaZ . prop[p1,me] . GammaBarNu
MsqZNu = Contract[Tr[LineZNu] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineZZ =  prop[p2,-me] . GammaZ . prop[p1,me] . GammaBarZ
MsqZZ =  Contract[Tr[LineZZ] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineHH =  prop[p2,-me] . GammaH . prop[p1,me] . GammaBarH
MsqHH =  Contract[Tr[LineHH] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineNuH =  prop[p2,-me] . GammaNu . prop[p1,me] . GammaBarH
MsqNuH = Contract[Tr[LineNuH] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineGH =  prop[p2,-me] . GammaG . prop[p1,me] . GammaBarH
MsqGH = Contract[Tr[LineGH] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps

LineZH =  prop[p2,-me] . GammaZ . prop[p1,me] . GammaBarH
MsqZH = Contract[Tr[LineZH] PVL[Q1,mu] PVL[Q2,nu] PVL[Q1,mup] PVL[Q2,nup]] /. subeps


MsqNu = MsqNuNu /. onshell
MsqNuLL = MsqNu /. PolLL
MsqNuLL = Simplify[MsqNuLL*g^4/4 /. ct] /. g^4->32 MW^4 GF^2

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

MsqNuGZ = MsqNuNu + MsqGG + MsqZZ + 2 MsqNuG + 2 MsqNuZ + 2 MsqZG;
MsqNuGZLL = MsqNuGZ /. onshell;
MsqNuGZLL = MsqNuGZLL /. PolLL;
MsqNuGZLL = Simplify[MsqNuGZLL*g^4/4 /. ct] /. g^4->32 MW^4 GF^2 ;


MsqHiggs = MsqHH + 2 MsqNuH + 2 MsqGH + 2 MsqZH ;
MsqAll = MsqNuGZ + MsqHiggs ;

MsqTotalAux = MsqAll /. onshell
MsqTotalLL = MsqTotalAux /. PolLL;
MsqTotalLT1 = MsqTotalAux /. PolLT1;
MsqTotalLT2 = MsqTotalAux /. PolLT2;
MsqTotalT1L = MsqTotalAux /. PolT1L;
MsqTotalT2L = MsqTotalAux /. PolT2L;
MsqTotalT1T1 = MsqTotalAux /. PolT1T1;
MsqTotalT1T2 = MsqTotalAux /. PolT1T2;
MsqTotalT2T1 = MsqTotalAux /. PolT2T1;
MsqTotalT2T2 = MsqTotalAux /. PolT2T2;

MsqLL=Simplify[MsqTotalLL /. steta->Sqrt[1-cteta^2]];
MsqLL=Simplify[MsqLL /. ct];
MsqLL=Simplify[MsqLL*g^4/4 /. weakangle] /. g^4->32 MW^4 GF^2 ;

MsqTL=Simplify[MsqTotalT1L+ MsqTotalT2L + MsqTotalLT1 + MsqTotalLT2 ];
MsqTL=Simplify[MsqTL /. steta->Sqrt[1-cteta^2]];
MsqTL=Simplify[MsqTL /. ct];
MsqTL=Simplify[MsqTL*g^4/4 /. weakangle] /. g^4->32 MW^4 GF^2;

MsqTT=Simplify[MsqTotalT1T1+MsqTotalT1T2+MsqTotalT2T1+MsqTotalT2T2 ]
MsqTT=Simplify[MsqTT /. steta->Sqrt[1-cteta^2]];
MsqTT=Simplify[MsqTT /. ct];
MsqTT=Simplify[MsqTT*g^4/4 /. weakangle] /. g^4->32 MW^4 GF^2

MsqTotal = Simplify[MsqLL + MsqTL + MsqTT]

(* Write Fortran Code *)

xs=1




If[xs==1,
   SetOptions[$Output, PageWidth -> 70];
   stmp=OpenWrite["ansNu.f",FormatType -> FortranForm];
   aux=Simplify[MsqNuLL /. simpbetas];
   Write[stmp,"MsqNu = ", aux];
   Close[stmp];
   stmp=OpenWrite["ansNuGZ.f",FormatType -> FortranForm];
   aux=Simplify[MsqNuGZLL /. simpbetas];
   Write[stmp,"MsqNuGZ = ", aux];
   Close[stmp];
   stmp=OpenWrite["ansLL.f",FormatType -> FortranForm];
   aux=Simplify[MsqLL /. simpbetas];
   Write[stmp,"MsqLL = ", aux];
   Close[stmp];
   stmp=OpenWrite["ansTL.f",FormatType -> FortranForm];
   aux=Simplify[MsqTL /. simpbetas];
   Write[stmp,"MsqTL = ", aux];
   Close[stmp];
   stmp=OpenWrite["ansTT.f",FormatType -> FortranForm];
   aux=Simplify[MsqTT /. simpbetas];
   Write[stmp,"MsqTT = ", aux];
   Close[stmp];
   stmp=OpenWrite["ansTotal.f",FormatType -> FortranForm];
   aux=Simplify[MsqTotal /. simpbetas];
   Write[stmp,"MsqTotal = ", aux];
   Close[stmp];
   stmp=OpenWrite["ansH.f",FormatType -> FortranForm];
   aux=Simplify[MsqHLL /. simpbetas];
   Write[stmp,"MsqH = ", aux];
   Close[stmp], Print[" "]];




Result = Function[{rs, cteta, exp}, factor = 3.893795*10^8; 
  beta = Sqrt[1 - 4 MW^2/rs^2] /. num;
  betae = Sqrt[1 - 4 me^2/rs^2] /. num;
  t = me^2 + MW^2 - rs^2/2*(1 - beta*betae*cteta) /. num; 
  ola = exp*factor/32/Pi/rs^2*beta/betae /. num; ola /. s -> rs^2]
