(*
Program to calculate the amplitudes squared for the process

  e^- e^+ -> W^- e^+ \nu_e

in the SM

For details see http://porthos.ist.utl.pt/CTQFT

Author: jorge.romao@ist.utl·pt

Last Revision: 26/04/2012
*)

Remove["Global`*"]

Batch=True

If[Batch,<<FC]

TimeNow=TimeUsed[];

(* Definitons *)

SumPolV=Function[{k,m,mu,nu},-MT[mu,nu]+FV[k,mu]*FV[k,nu]/m^2];
PL:=DiracGamma[7];
PR:=DiracGamma[6];

p5=p1+p2-p3-p4;

kin={SP[p1,p1]->0,SP[p2,p2]->0,SP[p3,p3]->0,SP[p4,p4]->mW^2,SP[p3,p4]->(-mW^2 -2 SP[p1,p2]+ 2 SP[p1,p3] + 2 SP[p1,p4]+ 2 SP[p2,p3]+2 SP[p2,p4])/2};



vecmom={SP[p1,p2]->p1p2,SP[p1,p3]->p1p3,SP[p1,p4]->p1p4,
SP[p2,p3]->p2p3,SP[p2,p4]->p2p4,SP[p3,p4]->p3p4,LC[][p1, p2, p3, p4]->Epsp1p2p3p4};

kinFCI={Pair[Momentum[p1], Momentum[p1]] -> 0, 
 Pair[Momentum[p2], Momentum[p2]] -> 0, 
 Pair[Momentum[p3], Momentum[p3]] -> 0, 
 Pair[Momentum[p4], Momentum[p4]] -> mW^2, 
 Pair[Momentum[p3], Momentum[p4]] -> 
  (-mW^2 - 2*Pair[Momentum[p1], Momentum[p2]] + 
           2*Pair[Momentum[p1], Momentum[p3]] + 
           2*Pair[Momentum[p1], Momentum[p4]] + 
           2*Pair[Momentum[p2], Momentum[p3]] + 
           2*Pair[Momentum[p2], Momentum[p4]])/2, gWWV[A]->gWWA,gWWV[Z]->gWWZ,gW->g2/Sqrt[2]};



Pol=Function[{a,b},a PL + b PR];

VWMWPA=Function[{al,be,ga,p1,p2,p3,V},
       gWWV[V] (MT[al,be] ((FourVector[p1,ga]) - (FourVector[p2,ga])) +
                MT[be,ga] ((FourVector[p2,al]) - (FourVector[p3,al])) +
                MT[ga,al] ((FourVector[p3,be]) - (FourVector[p1,be])))];


stmpAll=OpenWrite["Diags.f",FormatType->FortranForm,PageWidth->68];
stmpImag=OpenWrite["Imag.f",FormatType->FortranForm,PageWidth->68];


(* Propagators *)

SignVB=-1
SymmetryEven=+1
SymmetryOdd=-1

DenA1A=(-2 SP[p2,p3]) SignVB
DenA1El=(mW^2+2 SP[p4,p5])
DenA1 = Calc[DenA1A DenA1El SymmetryEven] /.kinFCI//Simplify

DenA2Z=(-mZ^2-2 SP[p2,p3]) SignVB
DenA2El=(mW^2+2 SP[p4,p5])
DenA2 =  Calc[DenA2Z DenA2El SymmetryEven] /.kinFCI//Simplify

DenA3Nu = (mW^2 -2 SP[p1,p4])
DenA3Z  = (-mZ^2 -2 SP[p2,p3]) SignVB
DenA3 = Calc[DenA3Nu DenA3Z SymmetryEven] /.kinFCI//Simplify

DenA4W = (-mW^2 -2 SP[p1,p5]) SignVB
DenA4A = (-2 SP[p2,p3]) SignVB
DenA4 = Calc[DenA4W DenA4A SymmetryEven] /.kinFCI//Simplify

DenA5W = (-mW^2 -2 SP[p1,p5]) SignVB
DenA5Z = (-mZ^2 -2 SP[p2,p3]) SignVB
DenA5 = Calc[DenA5W DenA5Z SymmetryEven] /.kinFCI//Simplify

DenA6W = (-mW^2 -2 SP[p1,p5]) SignVB
DenA6Nu = (mW^2+2 SP[p3,p4])
DenA6 = Calc[DenA6W DenA6Nu SymmetryEven] /.kinFCI//Simplify


DenB1A = (2 SP[p1,p2]) SignVB
DenB1W = ( (-mW^2 +2 SP[p3,p5]) + I (wW mW) ) SignVB
DenB1 = Calc[DenB1A DenB1W SymmetryOdd] /.kinFCI//Simplify

DenB2A = (2 SP[p1,p2]) SignVB
DenB2El = (mW^2 + 2 SP[p4,p5])
DenB2 = Calc[DenB2A DenB2El SymmetryOdd] /.kinFCI//Simplify


DenB3Z = ( (-mZ^2 + 2 SP[p1,p2]) + I (wZ mZ) ) SignVB
DenB3W = ( (-mW^2 +2 SP[p3,p5]) + I (wW mW) ) SignVB
DenB3 = Calc[DenB3Z DenB3W SymmetryOdd] /.kinFCI//Simplify


DenB4Z = ( (-mZ^2 + 2 SP[p1,p2]) + I (wZ mZ) ) SignVB
DenB4El = (mW^2 + 2 SP[p4,p5])
DenB4 = Calc[DenB4Z DenB4El SymmetryOdd] /.kinFCI//Simplify

DenB5Z = ( (-mZ^2 + 2 SP[p1,p2]) + I (wZ mZ) ) SignVB
DenB5Nu = (mW^2 + 2 SP[p3,p4])
DenB5 = Calc[DenB5Z DenB5Nu SymmetryOdd] /.kinFCI//Simplify

DenB6Nu = (mW^2 -2 SP[p1,p4]) 
DenB6W = ( (-mW^2 +2 SP[p3,p5]) + I (wW mW) ) SignVB
DenB6 = Calc[DenB6Nu DenB6W SymmetryOdd] /.kinFCI//Simplify

DenA={DenA1,DenA2,DenA3,DenA4,DenA5,DenA6};
DenB={DenB1,DenB2,DenB3,DenB4,DenB5,DenB6};
ReDenB=ComplexExpand[Re[DenB]];
ImDenB=ComplexExpand[Im[DenB]];
DenBc=ReDenB - I ImDenB; 


(* Definiton of fermion lines *)


Gamma1A1 = GA[mu] . (gW PL) . (GS[p4]+GS[p5]) . (ee GA[id1]) ; 
Gamma2A1 = (ee GA[id1]) ;
Gamma1A1c = (ee GA[id2])  . (GS[p4]+GS[p5]) . GA[mup] .(gW PL);
Gamma2A1c = (ee GA[id2])  ;

Gamma1A2 = GA[mu] .(gW PL) . (GS[p4]+GS[p5]) . GA[id3] . (gZL PL + gZR PR); 
Gamma2A2 = GA[id3] . ( gZL PL + gZR PR);
Gamma1A2c = GA[id4] . (gZL PL + gZR PR) . (GS[p4]+GS[p5]) . GA[mup] .(gW PL);
Gamma2A2c = GA[id4] . ( gZL PL + gZR PR);

Gamma1A3 =GA[id5] .(gZnuL PL + gZnuR PR) . (GS[p1]-GS[p4]). GA[mu] . (gW PL);
Gamma2A3 = GA[id5] . ( gZL PL + gZR PR); 
Gamma1A3c= GA[mup] . (gW PL).(GS[p1]-GS[p4]). GA[id6] .(gZnuL PL + gZnuR PR);
Gamma2A3c =GA[id6] . ( gZL PL + gZR PR); 

Gamma1A4 = GA[al] . (gW PL)  VWMWPA[al,mu,be,p1-p5,-p4,p2-p3,A];
Gamma2A4 = (ee GA[be]);
Gamma1A4c = GA[alp] . (gW PL)  VWMWPA[alp,mup,bep,p1-p5,-p4,p2-p3,A];
Gamma2A4c = (ee GA[bep]);

Gamma1A5 = GA[al] . (gW PL) VWMWPA[al,mu,be,p1-p5,-p4,p2-p3,Z];
Gamma2A5 = GA[be] . ( gZL PL + gZR PR);
Gamma1A5c= GA[alp] . (gW PL) VWMWPA[alp,mup,bep,p1-p5,-p4,p2-p3,Z];
Gamma2A5c= GA[bep] . ( gZL PL + gZR PR);

Gamma1A6 = GA[id7] . (gW PL);
Gamma2A6 = GA[id7] . (gW PL) . (-GS[p3]-GS[p4]) . GA[mu] . (gW PL);
Gamma1A6c = GA[id8] . (gW PL);
Gamma2A6c = GA[mup] . (gW PL) . (-GS[p3]-GS[p4]) . GA[id8] . (gW PL)  ;

Gamma1B1 = (ee GA[id9]);
Gamma2B1 = GA[id10] . (gW PL) VWMWPA[id10,mu,id9,-p3-p5,-p4,p1+p2,A];
Gamma1B1c = (ee GA[id11]);
Gamma2B1c = GA[id12] . (gW PL) VWMWPA[id12,mup,id11,-p3-p5,-p4,p1+p2,A];

Gamma1B2 = (ee GA[id13]);
Gamma2B2 = GA[mu] . (gW PL) . (GS[p4]+GS[p5]) . (ee GA[id13]);
Gamma1B2c = (ee GA[id14]);
Gamma2B2c = (ee GA[id14]) . (GS[p4]+GS[p5]) . GA[mup] . (gW PL)  ;

Gamma1B3 = GA[id15] . (gZL PL + gZR PR);
Gamma2B3 = GA[id16] . (gW PL)  VWMWPA[id16,mu,id15,-p3-p5,-p4,p1+p2,Z];
Gamma1B3c = GA[id17] . (gZL PL + gZR PR);
Gamma2B3c = GA[id18] . (gW PL) VWMWPA[id18,mup,id17,-p3-p5,-p4,p1+p2,Z];

Gamma1B4 = GA[id19] . (gZL PL + gZR PR);
Gamma2B4 = GA[mu] . (gW PL) . (GS[p4]+GS[p5]) . GA[id19] . (gZL PL + gZR PR);
Gamma1B4c = GA[id20] . (gZL PL + gZR PR);
Gamma2B4c = GA[id20] . (gZL PL + gZR PR) . (GS[p4]+GS[p5]) . GA[mup] . (gW PL);

Gamma1B5 = GA[id21] . (gZL PL + gZR PR);
Gamma2B5 = GA[id21] . (gZnuL PL + gZnuR PR).(-GS[p3]-GS[p4]) . GA[mu] . (gW PL);
Gamma1B5c = GA[id22] . (gZL PL + gZR PR);
Gamma2B5c = GA[mup] . (gW PL).(-GS[p3]-GS[p4]) . GA[id22].(gZnuL PL + gZmuR PR);

Gamma1B6 = GA[id23] . (gW PL) . (GS[p1]-GS[p4]) . GA[mu] . (gW PL);
Gamma2B6 = GA[id23] . (gW PL);
Gamma1B6c = GA[mup] . (gW PL) . (GS[p1]-GS[p4]) . GA[id24] . (gW PL);
Gamma2B6c = GA[id24] . (gW PL);

Gamma1A={Gamma1A1,Gamma1A2,Gamma1A3,Gamma1A4,Gamma1A5,Gamma1A6};
Gamma2A={Gamma2A1,Gamma2A2,Gamma2A3,Gamma2A4,Gamma2A5,Gamma2A6};
Gamma1Ac={Gamma1A1c,Gamma1A2c,Gamma1A3c,Gamma1A4c,Gamma1A5c,Gamma1A6c};
Gamma2Ac={Gamma2A1c,Gamma2A2c,Gamma2A3c,Gamma2A4c,Gamma2A5c,Gamma2A6c};

Gamma1B={Gamma1B1,Gamma1B2,Gamma1B3,Gamma1B4,Gamma1B5,Gamma1B6};
Gamma2B={Gamma2B1,Gamma2B2,Gamma2B3,Gamma2B4,Gamma2B5,Gamma2B6};
Gamma1Bc={Gamma1B1c,Gamma1B2c,Gamma1B3c,Gamma1B4c,Gamma1B5c,Gamma1B6c};
Gamma2Bc={Gamma2B1c,Gamma2B2c,Gamma2B3c,Gamma2B4c,Gamma2B5c,Gamma2B6c};

If[False,

(* Evaluate All Traces And Save *)


Trace1A1Ac=Table[Tr[GS[p5].Gamma1A[[i]]. Pol[a,b] . GS[p1].Gamma1Ac[[j]]],{i,1,6},{j,1,6}]/. kinFCI//Expand;
Trace2A2Ac=Table[Tr[Pol[c,d] . GS[p2].Gamma2A[[i]].GS[p3].Gamma2Ac[[j]]],{i,1,6},{j,1,6}]/. kinFCI//Expand;

file=Close[OpenWrite["TracesAA"]];
Save[file,{Trace1A1Ac,Trace2A2Ac}];


Trace1B1Bc=Table[Tr[Pol[c,d] . GS[p2].Gamma1B[[i]]. Pol[a,b] . GS[p1].Gamma1Bc[[j]]],{i,1,6},{j,1,6}]/. kinFCI//Expand;
Trace2B2Bc=Table[Tr[GS[p5].Gamma2B[[i]].GS[p3].Gamma2Bc[[j]]],{i,1,6},{j,1,6}]/. kinFCI//Expand;

file=Close[OpenWrite["TracesBB"]];
Save[file,{Trace1B1Bc,Trace2B2Bc}];



TraceABc=Table[Tr[GS[p5].Gamma1A[[i]]. Pol[a,b] . GS[p1].Gamma1Bc[[j]]. Pol[c,d] . GS[p2].Gamma2A[[i]].GS[p3].Gamma2Bc[[j]]],{i,1,6},{j,1,6}]/. kinFCI//Expand;
TraceBAc=Table[Tr[Pol[c,d] . GS[p2].Gamma1B[[i]]. Pol[a,b] . GS[p1].Gamma1Ac[[j]].GS[p5].Gamma2B[[i]].GS[p3].Gamma2Ac[[j]]],{i,1,6},{j,1,6}]/. kinFCI//Expand;


file=Close[OpenWrite["TracesABAndBA"]];
Save[file,{TraceABc,TraceBAc}];

]


If[False,

Get["TracesAA"];

Get["TracesBB"];

Get["TracesABAndBA"];

If[True,
MsqAAc=Table[Contract[
   Trace1A1Ac[[i, j]] Trace2A2Ac[[i, j]] SumPolV[p4, mW, mu, mup]] /. 
  kinFCI //Expand//Simplify, {i, 1, 6}, {j, 1, 6}];

file=Close[OpenWrite["MsqAA"]];
Save[file,{MsqAAc}];
]

If[True,

FunAux=Function[{i,j}, aux1=Contract[Trace2B2Bc[[i, j]] SumPolV[p4, mW, mu, mup]] /. kinFCI // Expand; Contract[Trace1B1Bc[[i, j]] aux1] /. kinFCI // Expand];

MsqBBc=Table[FunAux[i,j], {i, 1, 6}, {j, 1, 6}];

file=Close[OpenWrite["MsqBB"]];
Save[file,{MsqBBc}];
]

If[True,
MsqABc=Table[Contract[
   TraceABc[[i, j]] SumPolV[p4, mW, mu, mup]] /. kinFCI //Expand//Simplify, {i, 1, 6}, {j, 1, 6}];

file=Close[OpenWrite["MsqAB"]];
Save[file,{MsqABc}];
]

If[True,
MsqBAc=Table[Contract[
   TraceBAc[[i, j]] SumPolV[p4, mW, mu, mup]] /. 
  kinFCI //Expand//Simplify, {i, 1, 6}, {j, 1, 6}];

file=Close[OpenWrite["MsqBA"]];
Save[file,{MsqBAc}];
]

]


If[True,

Get["MsqAA"];
Get["MsqBB"];
Get["MsqAB"];
Get["MsqBA"];


Do[ 

If[i==j,
aux=MsqAAc[[i,i]]/DenA[[i]]/DenA[[i]]/. kinFCI//Expand//Simplify,
aux=(MsqAAc[[i,j]]+MsqAAc[[j,i]])/DenA[[i]]/DenA[[j]]/. kinFCI//Expand//Simplify];

Reaux = ComplexExpand[Re[aux]] //FCE//Simplify; 


Write[stmpAll,"  "];
Write[stmpAll,"* AA",i,j];
Write[stmpAll,"dsig = dsig + ", Reaux /. vecmom];

Imaux=ComplexExpand[Im[aux]]//FCE//Simplify; 
Write[stmpImag,"  "];
Write[stmpImag,"* AA",i,j];
Write[stmpImag,"dsig = dsig + ", Imaux /. vecmom];,{i,1,6},{j,i,6}];

Write[stmpAll,"  "];
Write[stmpAll,"* End AA"];
Write[stmpImag,"  "];
Write[stmpImag,"* End AA"];



Do[ 

If[i==j,
aux=MsqBBc[[i,i]]/DenB[[i]]/DenBc[[i]]/. kinFCI//Expand//Simplify,
aux=(MsqBBc[[i,j]]/DenB[[i]]/DenBc[[j]]+MsqBBc[[j,i]]/DenB[[j]]/DenBc[[i]])/. kinFCI//Expand//Simplify];

Reaux = ComplexExpand[Re[aux]] //FCE//Simplify; 

Write[stmpAll,"  "];
Write[stmpAll,"* BB",i,j];
Write[stmpAll,"dsig = dsig + ", Reaux /. vecmom];

Imaux=ComplexExpand[Im[aux]]//FCE//Simplify; 
Write[stmpImag,"  "];
Write[stmpImag,"* BB",i,j];
Write[stmpImag,"dsig = dsig + ", Imaux /. vecmom];,{i,1,6},{j,i,6}];

Write[stmpAll,"  "];
Write[stmpAll,"* End BB"];
Write[stmpImag,"  "];
Write[stmpImag,"* End BB"];


Do[ 

aux=(MsqABc[[i,j]]/DenBc[[j]] + MsqBAc[[j,i]]/DenB[[j]])/DenA[[i]]/. kinFCI//Expand//Simplify;

Reaux = ComplexExpand[Re[aux]]//FCE//Simplify;

Write[stmpAll,"  "];
Write[stmpAll,"* AB",i,j];
Write[stmpAll,"dsig = dsig + ", Reaux /. vecmom];

Imaux=ComplexExpand[Im[aux]]//FCE//Simplify; 
Write[stmpImag,"  "];
Write[stmpImag,"* AB",i,j];
Write[stmpImag,"dsig = dsig + ", Imaux /. vecmom];,{i,1,6},{j,1,6}];

Write[stmpAll,"  "];
Write[stmpAll,"* End AB"];
Write[stmpImag,"  "];
Write[stmpImag,"* End AB"];

]



If[Batch,Close[stmpAll]];
If[Batch,Close[stmpImag]];



Print["Time = ",TimeUsed[]-TimeNow];

