(******************************** Program nunuWLWL.m *************************************
    This program evaluates the amplitudes and cross section for the process

                      nu + nubar -> W^+_L + W^-_L

 Author: Jorge C. Romão
 email: jorge.romao@ist.utl.pt
********************************* Program nunuWLWL.m ************************************)
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=dm[7]
PR=dm[6]


(* 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 /.ct]; Simplify[TmpAux2 /. beta->Sqrt[1-4 MW^2/s]]] ; 

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



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

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


mandel={sp[p1,p2]-> s/2,sp[q1,q2]-> (s -2 MW^2)/2, sp[p1,q1]-> (MW^2 -t)/2,sp[p2,q2]-> ( MW^2-t)/2, sp[p1,q2]-> (s+t- MW^2)/2,sp[p2,q1]-> (s+t- MW^2)/2}

highenergy={s->4 MW^2 x,t-> MW^2 (1 - 2 x (1- Sqrt[(1 -1/x)] cteta)),beta->Sqrt[1-1/x]}
weakangle={sw^2->1 -cw^2}

ct = Part[Solve[t == MW^2 - s/2*(1 - beta*cteta), cteta], 1]

ampsZ= - 1/2 PVL[Q1,mu] PVL[Q2,nu] Spinor[-p2] . dm[a] .  PL .  Spinor[p1] V[mu,nu,b,-q1,-q2,q1+q2] (-mt[a,b] + fv[p1+p2,a] fv[p1+p2,b]/MW^2)/(s-MW^2/cw^2)

MsZ = Simplify[DiracSimplify[Contract[ampsZ] /. onshell] /. mandel]


ampNeu= - 1/2 Spinor[-p2] . dm[nu] . PL . prop[p1-q1,0] . dm[mu] . PL . Spinor[p1] PVL[Q1,mu] PVL[Q2,nu] 1/t

MtNeu= Simplify[DiracSimplify[Contract[ampNeu] /. onshell] /. mandel]



diracsub1={DiracGamma[Momentum[Q1]]->1/MW DiracGamma[Momentum[q1]],DiracGamma[Momentum[Q2]]->1/MW DiracGamma[Momentum[q2]]};

diracsub2={Spinor[-Momentum[p2], 0, 1] . (DiracGamma[Momentum[q1]]/MW) . DiracGamma[7] . Spinor[Momentum[p1], 0, 1]-> 1/MW Spinor[-Momentum[p2], 0, 1] . DiracGamma[Momentum[q1]] . DiracGamma[7] . Spinor[Momentum[p1], 0, 1],Spinor[-Momentum[p2], 0, 1] . (DiracGamma[Momentum[q2]]/MW) . DiracGamma[7] . Spinor[Momentum[p1], 0, 1]-> 1/MW Spinor[-Momentum[p2], 0, 1] . DiracGamma[Momentum[q2]] . DiracGamma[7] . Spinor[Momentum[p1], 0, 1]};



diracsub3={DiracGamma[Momentum[q2]]->DiracGamma[Momentum[p1+p2-q1]]}


MsZaux1= MsZ /. diracsub1
MsZaux2= Simplify[MsZaux1 /. diracsub2]
MsZaux3= MsZaux2 /. highenergy
MsZaux4= MsZaux3 /. x->1/z
MsZaux5= Normal[Series[MsZaux4,{z,0,1}]] 

MtNeuaux1 = DiracSimplify[MtNeu /. diracsub1]
MtNeuaux2 = MtNeuaux1 /. onshell
MtNeuaux3 = MtNeuaux2 /. highenergy
MtNeuaux4 = MtNeuaux3 /. x->1/z
MtNeuaux5 = Normal[Series[MtNeuaux4,{z,0,1}]] 

Mtotal= Simplify[MsZaux5 + MtNeuaux5]

MHighEnergy:= Simplify[DiracSimplify[Mtotal  /. diracsub3] /. z->1/x]

Print["1) High Energy Behaviour of the Amplitudes: \n x=s/(4 MW^2) , \[CurlyPhi](-p2).(\[Gamma]\[CenterDot]q1).\[Gamma]^7.\[CurlyPhi](p1) goes like s"] 
Print["MHighEnergy = ",MHighEnergy]


(* Cross Section *)
Print["\n 2) Exact Amplitudes Squared: \n"]

line1= ds[p2] . dm[nu] . PL . ds[p1-q1] . dm[mu] . PL . ds[p1] . dm[mup] . PL . ds[p1-q1] . dm[nup] . PL PVL[Q1,mu] PVL[Q1,mup] PVL[Q2,nu] PVL[Q2,nup]


cst1 = g^4/4/t^2
ans1= cst1 Contract[Tr[line1]] /. onshell
Msqt= Simplify[ans1 /. mandel] 

Print["|M_t|^2 = ",Msqt /. ct /. beta -> Sqrt[1 - 4 MW^2/s] // Simplify]

line2= ds[p2] . dm[a] . PL . ds[p1] . dm[b] . PL V[mu,nu,a,-q1,-q2,p1+p2] \
  V[mup,nup,b,-q1,-q2,p1+p2]  PVL[Q1,mu] PVL[Q1,mup] PVL[Q2,nu] PVL[Q2,nup] 

cst2= g^4/4/(s-MZ^2)^2
ans2= cst2 Contract[Tr[line2]] /. onshell
Msqs= Simplify[ans2 /. mandel] 

Print["|M_s|^2 = ",Msqs /. ct /. beta -> Sqrt[1 - 4 MW^2/s] // Simplify]

subeps={Eps[Momentum[a_],Momentum[b_],Momentum[c_],Momentum[d_]]->0}

line12=ds[p2] . dm[nu] . PL . (ds[p1]-ds[q1]) . dm[mu] . PL . ds[p1] . dm[a] . PL PVL[Q1,mu] PVL[Q1,mup] PVL[Q2,nu] PVL[Q2,nup] V[mup,nup,a,-q1,-q2,p1+p2]

line21=ds[p2] . dm[a] . PL . ds[p1] . dm[mu] . PL . (ds[p1]-ds[q1]) . dm[nu] . PL PVL[Q1,mu] PVL[Q1,mup] PVL[Q2,nu] PVL[Q2,nup] V[mup,nup,a,-q1,-q2,p1+p2]

cst12= -g^4/4/t/(s-MZ^2)
ans12=cst12 Contract[Tr[line12]+Tr[line21]] /. onshell
ans12= Simplify[ans12 /. mandel]
Msqstts= Simplify[ans12 /. subeps]

Print["M_s M_t^* + M_s^* M_t = ",Msqstts /. ct /. beta -> Sqrt[1 - 4 MW^2/s] // Simplify]



MsqTotal=Simplify[Msqt+Msqs+Msqstts];

(* High Energy Limit *)
Print["\n 3) High Energy Limit of Amplitudes Squared: \n"]

MsqtHEL=TakeLimit[Msqt] /. {cteta^2 g^4 x^2 -> g^4 x^2 ( 1 - steta^2), 
    cteta^2 g^4 x -> g^4 x (1 - steta^2)} // Simplify // Expand
Print["|M_t|^2 = ",MsqtHEL]

MsqsHEL=TakeLimit[Msqs] /. {cteta^2 g^4 x^2 -> g^4 x^2 ( 1 - steta^2), 
    cteta^2 g^4 x -> g^4 x (1 - steta^2)} // Simplify // Expand
Print["|M_s|^2 = ",MsqsHEL]

MsqsttsHEL=TakeLimit[Msqstts] /. {cteta^2 g^4 x^2 -> g^4 x^2 ( 1 - steta^2), 
    cteta^2 g^4 x -> g^4 x (1 - steta^2)} // Simplify // Expand
Print["2 Re(M_s M_t^*) = ",MsqsttsHEL]

MsqTotalHEL=MsqtHEL + MsqsHEL + MsqsttsHEL

MsqTotalHEL = MsqTotalHEL /. cteta^2->1-steta^2 //Simplify

Print["|M_total|^2 = ",MsqTotalHEL]

xs=1

If[xs==1,
   stmp=OpenWrite["anstotal.f",FormatType -> FortranForm];
   Write[stmp,MsqTotal /. g^4->32 MW^4 GF^2];
   Close[stmp];
   stmp=OpenWrite["ansneu.f",FormatType -> FortranForm];
   Write[stmp,Msqt /. g^4->32 MW^4 GF^2];
   Close[stmp];
   stmp=OpenWrite["ansZ.f",FormatType -> FortranForm];
   Write[stmp,Msqs /. g^4->32 MW^4 GF^2];
   Close[stmp];
   stmp=OpenWrite["ansneuZ.f",FormatType -> FortranForm];
   Write[stmp,Msqstts /. g^4->32 MW^4 GF^2];
   Close[stmp],Print[" "]];




