rm(list=ls())

dt <- 0.1					# ps
t <- seq(0, by=dt, length.out=2^18)	# time in ns
s <- rep(0, length(t))			# signal

Dt <- rexp(10000, 1e6)			# time intervals (in s)
T.det <- cumsum(Dt)*1e9+1000		# detection time (in ns)
T.det <- T.det[which(T.det < t[length(t)])]

PE.resp <- function () rnorm (1, 1, 0.05)

Nsec <- function () rpois(1, 0.4)

for (i in 1:length(T.det)) {
	secondaries <- Nsec ()
	index <- round(T.det[i]/dt)
	for (j in 1:(1+secondaries)) s[index] <- s[index]+ PE.resp()
}

plot (t, s, type="l")

cur <- s * 1.6e-19 * 1e6 / (dt * 1e-9) 

el.resp <- 1/3e-9*exp(-t/3)

det.resp <- 1/50e-9*exp(-t/50) 

RESP <- fft(det.resp) * fft(el.resp)
resp <- Re(fft(RESP, inverse = TRUE)) / length(t)
resp <- resp/sum(resp)

s.f <- 50 * 100 * Re(fft(fft(cur) * fft(resp), inverse = TRUE)) / length(t)
plot (t, s.f, type="l", xlab="time (ns)", ylab="amplitude (V)")

plot (t, s.f, type="l", xlab="time (ns)", ylab="amplitude (V)", xlim=c(14000, 18000))
grid()