-- Draw probability density for Bhabha scattering pi = float(pi) -- use numerical value of pi N = 4 -- number of bins a = pi / N M = (cos(theta)^2 + 3)^2 / (cos(theta) - 1)^2 f = M sin(theta) -- unnormalized probability density -- integrate f from a to b I(a,b) = do( dtheta = (b - a) / 100, sum(k,1,100,eval(f, theta, a + (k - 0.5) dtheta)) dtheta ) "Probability density" C = I(a,pi) -- normalization constant f = M sin(theta) / C f xrange = (a,pi) yrange = (0,4) draw(f,theta) "Probability distribution" f = M sin(theta) P = zero(N) for(k,2,N, theta2 = k pi/N, theta1 = theta2 - pi/N, P[k] = I(theta1,theta2) ) P = P / sum(P) -- normalize P Run