<< ~/SERVER/FYZIKA/MATHEMATICA/PARTIALWAVE/partial_waves.m

(* in p=proj=1 scenario Gama3 is contribution of rho -> bb. Gama3=Gama2 in the limit of negligible top and bottom masses *)
(* in proj=0, b_1 = 0 scenario replace the Gama3 line with Gama3 = 0 *)
Gama1 = 4 Gpi^2/(192 Pi) MRho^3/vev^2;
Gama2 = MRho/(4 Pi) (GVt^2 + GAt^2);
Gama3 = MRho/(4 Pi) (GVb^2 + GAb^2);
GamaRho = Gama1 + Gama2 + Gama3;

(* tt -> tt helicity amplitudes *)

Mpmpm = (-2*s*(GAt + GVt)^2*Cos[theta/2]^2)/(-MRho^2 + s + I MRho GamaRho) + (2*s*(GAt + GVt)^2*Cos[theta/2]^2)/(-MRho^2 - (s*(1 - Cos[theta]))/2);
Mmpmp = (-2*s*(GAt - GVt)^2*Cos[theta/2]^2)/(-MRho^2 + s + I MRho GamaRho) + (2*s*(GAt - GVt)^2*Cos[theta/2]^2)/(-MRho^2 - (s*(1 - Cos[theta]))/2);
Mpmmp = 2*s*(-GAt^2 + GVt^2)*Sin[theta/2]^2/(-MRho^2 + s + I MRho GamaRho); 
Mmppm = 2*s*(-GAt^2 + GVt^2)*Sin[theta/2]^2/(-MRho^2 + s + I MRho GamaRho) ;
Mpppp = 2*s*(-GAt^2 + GVt^2)/(-MRho^2 - (s*(1 - Cos[theta]))/2);
Mmmmm = 2*s*(-GAt^2 + GVt^2)/(-MRho^2 - (s*(1 - Cos[theta]))/2);

a0pppp = Simplify[PWA0[0,1/2,1/2,1/2,1/2,Mpppp]];
a0mmmm = Simplify[PWA0[0,-1/2,-1/2,-1/2,-1/2,Mmmmm]];
a1pmpm = Simplify[PWA0[1,1/2,-1/2,1/2,-1/2,Mpmpm]];
a1mpmp = Simplify[PWA0[1,-1/2,1/2,-1/2,1/2,Mmpmp]];
a1pmmp = Simplify[PWA0[1,1/2,-1/2,-1/2,1/2,Mpmmp]];
a1mppm = Simplify[PWA0[1,-1/2,1/2,1/2,-1/2,Mmppm]];
a1pppp = Simplify[PWA0[1,1/2,1/2,1/2,1/2,Mpppp]];
a1mmmm = Simplify[PWA0[1,-1/2,-1/2,-1/2,-1/2,Mmmmm]];

(* plot partial waves as function of s and GV for GA = 0 *)

GVt = -BR*G''/8. - BL*G''/8.;
GAt = -BR*G''/8. + BL*G''/8.;
GVb = +p^2*BR*G''/8. + BL*G''/8.;
GAb = +p^2*BR*G''/8. - BL*G''/8.;
GVtb = -p*BR*G''/(4.*Sqrt[2]) - BL*G''/(4.*Sqrt[2]); 
GAtb = -p*BR*G''/(4.*Sqrt[2]) + BL*G''/(4.*Sqrt[2]);
Gpi = MRho/(vev*G'');
s = ECM^2;
MRho = 1.0;
mt = 0.175;
vev = 0.246;
(*ECM = 2.5;*)
G'' = 100.;
p = 1.;

Aa0pppp = Sqrt[Re[a0pppp]^2 + Im[a0pppp]^2];
Aa0mmmm = Sqrt[Re[a0mmmm]^2 + Im[a0mmmm]^2];
Aa1pmpm = Sqrt[Re[a1pmpm]^2 + Im[a1pmpm]^2];
Aa1mpmp = Sqrt[Re[a1mpmp]^2 + Im[a1mpmp]^2];
Aa1pmmp = Sqrt[Re[a1pmmp]^2 + Im[a1pmmp]^2];
Aa1mppm = Sqrt[Re[a1mppm]^2 + Im[a1mppm]^2];
Aa1pppp = Sqrt[Re[a1pppp]^2 + Im[a1pppp]^2];
Aa1mmmm = Sqrt[Re[a1mmmm]^2 + Im[a1mmmm]^2];

(*Plot3D[Aa1pmpm,{ECM,0.5,3},{GA,0,3}, PlotRange -> {0,4}]*)

ECM = 2.5;
p1=ContourPlot[Aa0pppp,{BL,-1,1},{BR,-1,1}, PlotRange -> {1,100}, Contours ->{1}, PlotPoints -> {100,100}]

(*ContourPlot[Aa0pppp,{BL,0,3},{BR,0,3},PlotRange -> {0,100},Contours->{1},PlotPoints -> {100,100}, ColorFunction->(If[#1>1.,Red,White]&)]*)

(*ContourPlot[{Evaluate[Aa0pppp] ==1,Evaluate[Aa1pmpm]==1,Evaluate[Aa1mpmp]==1,Evaluate[Aa1pmmp]==1,Evaluate[Aa1pppp]==1},{BL,-3,3},{BR,-3,3}, PlotRange -> {0,1}, Contours ->{1}, PlotPoints -> {200,200}]*)

(* Blank plot means zero partial wave *)

p2=ContourPlot[Evaluate[Aa1pmpm],{BL,-1,1},{BR,-1,1}, PlotRange -> {1,100}, Contours ->{1}, PlotPoints -> {100,100}, AxesLabel -> {"ECM (TeV)", "GVt = GAt"}, FrameLabel -> {"", ""} ]

p3=ContourPlot[Evaluate[Aa1mpmp],{BL,-1,1},{BR,-1,1}, PlotRange -> {1,100}, Contours ->{1}, PlotPoints -> {100,100}, AxesLabel -> {"ECM (TeV)", "GVt = GAt"}, FrameLabel -> {"", ""} ]

p4=ContourPlot[Evaluate[Aa1pmmp],{BL,-10,10},{BR,-10,10}, PlotRange -> {1,100}, Contours ->{1}, PlotPoints -> {100,100}, AxesLabel -> {"ECM (TeV)", "GVt = GAt"}, FrameLabel -> {"", ""} ]

p5=ContourPlot[Evaluate[Aa1pppp],{BL,-1,1},{BR,-1,1}, PlotRange -> {1,100}, Contours ->{1}, PlotPoints -> {100,100}, AxesLabel -> {"ECM (TeV)", "GVt = GAt"}, FrameLabel -> {"", ""} ]

Show[p1,p2,p3,p4,p5]