decay.rate <- 60
Navg <- 100
gain <- 1e6

dt <- 0.01	#timestep (ns)
t <- seq (0, by=dt, length.out=2^16)
cur <- rep (0, length(t))

t[length(t)]

N.exp <- rpois(1, Navg)

t.PE <- rexp(N.exp, 1/decay.rate)
s.PE <- rpois(N.exp, 4) / 4 * gain

for (i in 1:length(t.PE)) {
	index <- t.PE[i] / dt
	if (t.PE[i] < t[length(t)]) {
		cur[index] <- cur[index] + s.PE[i] * 1.6e-19 / (dt * 1e-9)
	}	
}

#plot (t, cur, type="l")

##shape of single PE

PM.resp <- t/1*exp(-t/1)
PM.resp <- PM.resp / sum(PM.resp)

plot (t, PM.resp, type = "l", xlim=c(0, 20))

V <- 50 * Re(fft(fft(cur) * fft(PM.resp), inverse=TRUE))/length(cur)
plot (t, V, type="l")

V <- V + rnorm(length(t), 0, 0.0003)