-- Muon pair production data from SLAC PEP

-- Compute R-squared for electroweak model

-- www.hepdata.net/record/ins216031 (Table 1, 29.0 GeV)

N = 20 -- number of observations

-- x is cos(theta)

x = (
-0.925,
-0.85,
-0.75,
-0.65,
-0.55,
-0.45,
-0.35,
-0.25,
-0.15,
-0.05,
0.05,
0.15,
0.25,
0.35,
0.45,
0.55,
0.65,
0.75,
0.85,
0.925)

-- y is (2E)^2 times differential cross section (2E = 29 GeV)

y = (
67.08,
58.67,
54.66,
51.72,
43.7,
41.12,
39.71,
35.34,
33.35,
34.69,
34.05,
34.48,
34.66,
35.23,
35.6,
40.13,
42.56,
46.37,
49.28,
55.7)

y = y "nanobarn" "GeV"^2

"Observed cross section per bin"

"Predicted cross section"

pi = float(pi)
alpha = 7.2973525693 10^(-3)
hbar = 6.582119569 10^(-25) "GeV" "second"
c = 299792458.0 "meter" / "second"
C = 10^37 "nanobarn" / "meter"^2 -- convert square meters to nanobarns

gA = -0.5
gV = -0.0348
mZ = 91.17 "GeV"
G1 = 1.166 10^(-5) "GeV"^(-2) -- Fermi coupling constant

s = (29.0 "GeV")^2

F = alpha^2/4 (1 + gV^2/(sqrt(2) pi) mZ^2/(s - mZ^2) s G1/alpha +
    (gA^2 + gV^2)^2/(8 pi^2) (mZ^2/(s - mZ^2))^2 (s G1/alpha)^2)

G = alpha^2/4 (sqrt(2) gA^2/pi mZ^2/(s - mZ^2) s G1/alpha +
    gA^2 gV^2/pi^2 (mZ^2/(s - mZ^2))^2 (s G1/alpha)^2)

yhat = zero(N)

for(k,1,N,
  yhat[k] = C 2 pi (F (1 + x[k]^2) + G x[k]) (hbar c)^2
)

yhat

"Coefficient of determination (R squared)"

ybar = sum(y) / N

RSS = sum((y - yhat)^2) -- residual sum of squares
TSS = sum((y - ybar)^2) -- total sum of squares

1 - RSS / TSS