-- Muon pair production data from SLAC PEP

-- Compute R-squared for QED 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) "nanobarn" "GeV"^2

"Observed cross section per bin"

"Predicted cross section"

alpha = 7.2973525693 10^(-3)
c = 299792458.0 meter / second
h = 6.62607015 10^(-34) joule second
hbar = h / float(2 pi)
eV = 1.602176634 10^(-19) joule

pi = float(pi) -- use numerical value of pi

C1 = pi alpha^2 / 2
C2 = (hbar c)^2
C3 = 10^37 "nanobarn" / meter^2
C4 = 10^9 eV / "GeV"

C = C1 C2 C3 / C4^2

yhat = zero(N)

for(k,1,N,
 yhat[k] = C (1 + x[k]^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