rm(list=ls())

L <- 100	#cm
Th <- 10	#cm
E <- 10	#MeV
m <- 940	#MeV/c2
c <- 30	#cm/ns

v <- function (E) sqrt(2*E/m)
v(10)

T <- function (E) L/(v(E)*30)
T(10)

det.Th.contrib <- function (E) m*L/c^2/T(E)^2*Th/sqrt(12)
stop.contrib <- function (E) m*L^2/c^2/T(E)^3*0.2
start.contrib <- function (E) m*L^2/c^2/T(E)^3*0.170
tot.contrib <- function(E) sqrt(det.Th.contrib(E)^2 +
	stop.contrib(E)^2 + start.contrib(E)^2)

det.Th.contrib(10)
stop.contrib(10)
start.contrib(10)
tot.contrib(10)

##
## same proble, with MC
##

N <- 10000
start.p <- rnorm (N, 0, 0.17)
stop.p <- rnorm (N, T(10), 0.2)
L.p.expo <- function () {
	i.p <- rep(0, N)
	for (i in 1:N) {
		repeat {
			dL <- rexp(1, -log(0.8)/Th)
			if (dL < Th) break
		}			
		i.p[i] <- L + dL - Th/2 	
	}
	i.p
}
L.p <- L.p.expo ()

E <- 0.5 * m * (L.p/(stop.p-start.p)/c)^2
hist(E,freq=FALSE)

E.n <- seq(8, 12, by=0.1)
lines (E.n, dnorm(E.n, 10, 0.630), col="red", lwd=2)


##
## contribution at various energies
##

E <- seq(5,40, by=0.1)

plot (c(5,40), c(0,3), type="n",
	xlab="neutron energy (MeV)", ylab="uncertainty (MeV RMS)")
grid()
lines (E, tot.contrib(E), lwd=2)
lines (E, det.Th.contrib(E), col="red", lwd=2)
lines (E, start.contrib(E), col="blue", lwd=2)
lines (E, stop.contrib(E), col="green", lwd=2)

det.Th.contrib(10)/start.contrib(10)
L <- 200
det.Th.contrib(10)/start.contrib(10)

L <- 100
tot.contrib(10)
L <- 200
tot.contrib(10)