Clear[ga,gv];

coef=2*(1.16637*10^-5)^2*0.511*10^-3*0.3894*10^-27;(*  This coefficient if just 2*((Subscript[G, F])^2)*(Subscript[m, e])*(Subscript[E, ν]) in units of (1/GeV)(cm^2)  *)
coef=coef*10^42 ;(*  The 10^42 normalizes the number in units of the measured quantities in the experiment, see the problem sheet    *)
(**)
(*  The following are the cross sections computed: *)
σνμ=(coef/Pi)((gv+ga)^2+(((gv-ga)^2)/3));
σνbμ=(coef/Pi)((gv-ga)^2+(((gv+ga)^2)/3));
σνe=(coef/Pi)((1+gv+ga)^2+(((gv-ga)^2)/3));
σνbe=(coef/Pi)((gv-ga)^2+(((gv+ga+1)^2)/3));
(**)
(*  Experimental results: m is mean, δ is 1σ uncertainty in the measurement  *)
mσνμ=1.1;
δσνμ=0.6;
mσνbμ=2.2;
δσνbμ=1.0;
mσνbe=5.0;
δσνbe=1.4;
(**)
(*  The (χ^2) functions. The are build as ((( prediction - measurement ))^2)/((δunc)^2)  *)
χ2νμ=(((σνμ-mσνμ)^2)/(δσνμ^2));
χ2νbμ=(((σνbμ-mσνbμ)^2)/(δσνbμ^2));
χ2νbe=(((σνbe-mσνbe)^2)/(δσνbe^2));

Show[
RegionPlot[χ2νμ<1,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Red,Opacity[0.2]],BoundaryStyle->Red,Frame->True,FrameLabel->{"(Subscript[g, A])","(Subscript[g, V])"},FrameStyle->Directive[Black,FontSize->16],PlotRangePadding->None],
RegionPlot[χ2νbμ<1,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Blue,Opacity[0.2]],BoundaryStyle->Blue],
RegionPlot[χ2νbe<1,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Green,Opacity[0.2]],BoundaryStyle->Green]
]

(*  Now the total (χ^2) with the minimum substracted:  *)
χtotal=χ2νμ+χ2νbμ+χ2νbe;
min=NMinimize[χtotal,{ga,gv}];
χtotal=χtotal-min[[1]]//Expand;

(*  Now plot the 68%CL regions for the experiments, and the 68% and 95% for the combination:  *)
s12dof=2.3; (*  why? see below  *)
s22dof=5.991;
Show[
RegionPlot[χ2νμ<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Red,Opacity[0.2]],BoundaryStyle->Red,Frame->True,FrameLabel->{"(Subscript[g, A])","(Subscript[g, V])"},FrameStyle->Directive[Black,FontSize->16],PlotRangePadding->None],
RegionPlot[χ2νbμ<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Blue,Opacity[0.2]],BoundaryStyle->Blue],
RegionPlot[χ2νbe<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Green,Opacity[0.2]],BoundaryStyle->Green],
RegionPlot[χtotal<s22dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Gray,Opacity[0.1]],BoundaryStyle->{Black,Dashed}],
RegionPlot[χtotal<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Gray],BoundaryStyle->Black]
]

(*   Now use the Electroweak theory, or model, to write gv and ga in terms of a single parameters (Subscript[sinθ, W]):   *)
gv=((-1+4*sw2)/4);
ga=-1/4;
(*  Build the single-variable (χ^2):  *)
χ2νμmin=χ2νμ-NMinimize[χ2νμ,sw2][[1]];
χ2νbμmin=χ2νbμ-NMinimize[χ2νbμ,sw2][[1]];
χ2νbemin=χ2νbe-NMinimize[χ2νbe,sw2][[1]];

Show[
Plot[χ2νμmin,{sw2,0,1},PlotStyle->Red,PlotRange->{0,10}],
Plot[χ2νbμmin,{sw2,0,1},PlotStyle->Blue],
Plot[χ2νbemin,{sw2,0,1},PlotStyle->Green]
]

(*  Build the total (χ^2):  *)
χ2νtot=χ2νμ+χ2νbμ+χ2νbe-NMinimize[χ2νμ+χ2νbμ+χ2νbe,sw2][[1]];

(*  Check where its minima is (the point that fits best the data) and the uncertainty:  *)
NMinimize[χ2νtot,sw2][[2,1,2]]
Select[NSolve[χ2νtot==1,sw2],Im[#[[1,2]]]==0&&Re[#[[1,2]]]>0&]

Show[
Plot[χ2νtot,{sw2,0,0.6},PlotStyle->Black,PlotRange->{0,10}],
Plot[{1,4,9},{x,0,1},PlotStyle->Directive[Gray,Dashed]]
]

Include ee to μμ data

Clear[ga,gv]

parameters={GF->1.16637*10^-5,eem->(Sqrt[4Pi αem]),αem->1/137};

Aparam=1-((8(Sqrt[2])GF)/(eem^2))s gv^2+32 GF^2 s^2 (ga^2+gv^2)^2//.parameters
Bparam=-((16(Sqrt[2])GF s)/(eem^2))ga^2+256GF^2s^2ga^2gv^2//.parameters

σtot=((4Pi αem^2)/3)(1/s)*0.3894*10^9*Aparam//.parameters;
Aasym=0.3445*(Bparam/Aparam)//.parameters;

χtot=((((σtot/.s->30^2)-96)^2)/((4)^2))+((((σtot/.s->35^2)-71)^2)/((6)^2));
χasym=((((Aasym/.s->30^2)-(-0.03))^2)/((0.1)^2))+((((Aasym/.s->35^2)-0.006)^2)/((0.073)^2));

minasym=NMinimize[χasym,{ga,gv}];
χasym=χasym-minasym[[1]];

mintot=NMinimize[χtot,{ga,gv}];
χtot=χtot-mintot[[1]];

s12dof=2.3; 
Show[
RegionPlot[χtot<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Red,Opacity[0.2]],BoundaryStyle->Red,Frame->True,FrameLabel->{"(Subscript[g, A])","(Subscript[g, V])"},FrameStyle->Directive[Black,FontSize->16],PlotRangePadding->None],
RegionPlot[χasym<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Blue,Opacity[0.2]],BoundaryStyle->Blue]
]

χtotal=χ2νμ+χ2νbμ+χ2νbe+χtot+χasym;
min=NMinimize[χtotal,{ga,gv}][[1]];
χtotal=χtotal-min;

(*  Now plot the 68%CL regions for the experiments, and the 68% and 95% for the combination:  *)
s12dof=2.3; (*  why? see below  *)
s22dof=5.991;
Show[
RegionPlot[χ2νμ<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Red,Opacity[0.2]],BoundaryStyle->Red,Frame->True,FrameLabel->{"(Subscript[g, A])","(Subscript[g, V])"},FrameStyle->Directive[Black,FontSize->16],PlotRangePadding->None],
RegionPlot[χ2νbμ<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Blue,Opacity[0.2]],BoundaryStyle->Blue],
RegionPlot[χ2νbe<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Green,Opacity[0.2]],BoundaryStyle->Green],
RegionPlot[χtot<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Yellow,Opacity[0.2]],BoundaryStyle->Yellow],
RegionPlot[χasym<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Orange,Opacity[0.2]],BoundaryStyle->Orange],
RegionPlot[χtotal<s22dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Gray,Opacity[0.1]],BoundaryStyle->{Black,Dashed}],
RegionPlot[χtotal<s12dof,{ga,-1,1},{gv,-1,1},PlotStyle->Directive[Gray],BoundaryStyle->Black]
]

gv=((-1+4*sw2)/4);
ga=-1/4;
χtotal=χ2νμ+χ2νbμ+χ2νbe+χtot+χasym;
min=NMinimize[χtotal,sw2][[1]];
χtotal=χtotal-min;

Show[
Plot[χtotal,{sw2,0,0.6},PlotStyle->Black,PlotRange->{0,10}],
Plot[χ2νtot,{sw2,0,0.6},PlotStyle->Gray,PlotRange->{0,10}],
Plot[{1,4,9},{x,0,1},PlotStyle->Directive[Gray,Dashed]]
]