[Alpha]em = 1/137
me = 0.511/1000
mmu = 102.8/1000
F[x_] := (5*x^4 - 14*x^3 + 39*x^2 - 38*x - 18*x^2*Log[x] +
8)/(12*(1 - x)^4)
G[x_] := (x^3 + 3*x - 6*x*Log[x] - 4)/(2*(1 - x)^3)
[CapitalGamma]totalmu = 3*10^-19
M4C = 200
[CapitalDelta]amumax = (26.1 + 8.0*2)*10^-10
[CapitalDelta]amumin = (26.1 - 8.0*2)*10^-10

c12L = 1
s14L = [Pi]/4

gL[Mu][Mu][s24L2_,
s12L_] := (c12L*Sqrt[s24L2] - Sqrt[1 - s24L2]*s12L*s14L)^2;
gL[Mu]e[s24L2_,
s12L_] := (s12L*Sqrt[s24L2] +
 c12L*Sqrt[1 - s24L2]*s14L)*(s12L*Sqrt[s24L2] -
 Sqrt[1 - s24L2]*s12L*s14L);
gLee[s24L2_,
s12L_] := (s12L*Sqrt[s24L2] + c12L*Sqrt[1 - s24L2]*s14L)^2;
gL[Mu]E[s24L2_,
s12L_] := (Sqrt[1 - s14L^2]*
 Sqrt[1 - s24L2])*(Sqrt[1 - s12L^2]*Sqrt[s24L2] -
 Sqrt[1 - s24L2]*s12L*s14L);
gLeE[s24L2_,
s12L_] := (Sqrt[1 - s14L^2]*Sqrt[1 - s24L2])*(s12L*Sqrt[s24L2] +
 Sqrt[1 - s12L^2]*Sqrt[1 - s24L2]*s14L);
gLEE[s24L2_, s12L_] := (1 - s14L^2)*(1 - s24L2);

c12R = 1
s14R = 0

gR[Mu][Mu][s24R2_,
s12R_] := (c12R*Sqrt[s24R2] - Sqrt[1 - s24R2]*s12R*s14R)^2;
gR[Mu]e[s24R2_,
s12R_] := (s12R*Sqrt[s24R2] +
 c12R*Sqrt[1 - s24R2]*s14R)*(s12R*Sqrt[s24R2] -
 Sqrt[1 - s24R2]*s12R*s14R);
gRee[s24R2_,
s12R_] := (s12R*Sqrt[s24R2] + c12R*Sqrt[1 - s24R2]*s14R)^2;
gR[Mu]E[s24R2_,
s12R_] := (Sqrt[1 - s14R^2]*
 Sqrt[1 - s24R2])*(Sqrt[1 - s12R^2]*Sqrt[s24R2] -
 Sqrt[1 - s24R2]*s12R*s14R);
gReE[s24R2_,
s12R_] := (Sqrt[1 - s14R^2]*Sqrt[1 - s24R2])*(s12R*Sqrt[s24R2] +
 Sqrt[1 - s12R^2]*Sqrt[1 - s24R2]*s14R);
gREE[s24R2_, s12R_] := (1 - s14R^2)*(1 - s24R2);

[CapitalDelta]a[Mu]Zprime[MZprime_, s24L2_, s12L_, s24R2_,
s12R_] := -(mmu^2/(8 [Pi]^2 MZprime^2))*(
((gL[Mu][Mu][s24L2, s12L])^2 + (gR[Mu][Mu][s24R2, s12R])^2)*
 F[mmu^2/MZprime^2]
+ ((gL[Mu]e[s24L2, s12L])^2 + (gR[Mu]e[s24R2, s12R])^2)*
 F[me^2/MZprime^2]
+ ((gL[Mu]E[s24L2, s12L])^2 + (gR[Mu]E[s24R2, s12R])^2)*
 F[M4C^2/MZprime^2]
+ mmu/mmu*
 Re[gL[Mu][Mu][s24L2,
    s12L]*(gR[Mu][Mu][s24R2, s12R])[Conjugate]]*
 G[mmu^2/MZprime^2]
+ me/mmu*
 Re[gL[Mu]E[s24L2, s12L]*(gR[Mu]E[s24R2, s12R])[Conjugate]]*
 G[me^2/MZprime^2]
+ M4C/mmu*
 Re[gL[Mu]E[s24L2,
    s12L]*(gR[Mu][Mu][s24R2, s12R])[Conjugate]]*
 G[1001^2/MZprime^2])

ContourPlot[[CapitalDelta]a[Mu]Zprime[MZprime, s24L2, 10^-3, 10^-3,
10^-3], {MZprime, 50, 1000}, {s24L2, 0, 1},
ScalingFunctions -> {"Log", "Log"},
FrameLabel -> {Style["!(*SubscriptBox[(M), (Z')])[GeV]",
FontSize -> 16],
Style["!(*SuperscriptBox[(sin), (2)])!(*SubscriptBox[(
[Theta]), (24  L)])", FontSize -> 16]},
BaseStyle -> {FontFamily -> "Times", FontSize -> 14},
RegionFunction ->
Function[{MZprime, s24L2, 10^-3, 10^-3,
10^-3} , [CapitalDelta]a[Mu]Zprime[MZprime, s24L2, 10^-3, 10^-3,
  10^-3] <= [CapitalDelta]amumin],
PlotRange -> {{50, 1000}, {10^-3, 1.0}, All},
ContourShading -> {Opacity[0.4, Blue]}, Contours -> {10},
MaxRecursion -> 5]