windows (5,3)
par(mar=c(4.2, 4.2, 1, 1))
require (signal)


## joli PU de 2.5mV, Rise Time=2ns, Fall Time=3ns
## t en ns
signal.PU <- function (t) (t/1)^2*exp(-(t/1)^1.1)/1.435065e-9*1.6e-19*1e6*50

## distri de PU gaussienne tronquée
distri.PU <- function (n) {
	A <- rep(-1, length.out=n)
	for (i in 1:n) 
		while (A[i] < 0) A[i] <- rnorm(1, 1/1.28, 1/1.28)		
	A
}

## distri lumière 
rlight.n <- function (NPE) {
	nl <- rbinom(1, NPE, 0.5)
	c (rexp (NPE-nl, 1/9), rexp(nl, 1/60)) 
}

rlight.g <- function (NPE) {
	nl <- rbinom(1, NPE, 0.2)
	c (rexp (NPE-nl, 1/9), rexp(nl, 1/60)) 
}

rlight.t <- function (NPE) runif (NPE,0, 100)

# on suppose explicitement que dt=0.1 (100ps)
simu.scinti <- function (NPE_, t, vrms=1e-3, type="BC400") {

	NPE <- rpois(1,NPE_)
	#NPE <- NPE_
	s <- rep(0, length(t))	#signal de sortie
	n0 <- which(t>=0)[1]	#t=0

	tPU <- seq (0, 10, by=0.1)
	sPU <- signal.PU(tPU)	# production su signal PU sur 10ns
	#tPU <- 1
	#sPU <- c(1)
	if (type == "n") t.PE <- rlight.n (NPE)
	if (type == "g") t.PE <- rlight.g (NPE)
	if (type == "t") t.PE <- rlight.t (NPE)
	A.PE <- distri.PU (NPE)

	for (i in 1:NPE) {
		ns <- n0 + t.PE[i]/0.1
		ne <- ns + length(tPU)-1
		if (ne < length(s)) s[ns:ne] <- s[ns:ne] + sPU * A.PE[i]
	}

	#noise <- rnorm(length(s), 0, vrms)
	#bf <- butter (1, 500/5000)
	#s <- s + filter (bf, noise)
	s
}

t <- seq(-10,200,by=0.1)

## signaux beaux
Ntot <- 50000
n.good <- simu.scinti (Ntot, t, type="n")
g.good <- simu.scinti (Ntot, t, type="g")
plot (c(0,200), c(0, max(g.good)), type="n")
lines (t, g.good, col="blue")
lines (t, n.good, col="red")

## affichage de quelques signaux
Ntot <- 1000
plot (c(0,200), c(0, 0.6), type="n", xlab="time (ns)", ylab="amplitude (V)")
grid()
for (i in 1:10) lines (t, simu.scinti (Ntot, t, type="g"), col="blue")
for (i in 1:10) lines (t, simu.scinti (Ntot, t, type="n"), col="red")
#for (i in 1:10) lines (t, simu.scinti (Ntot, t, type="t"), col="green")


discri.ng <- function(M, VI) {
	Q <- rep(0, dim(M)[1])
	for (i in 1:dim(M)[1]) Q[i] <- sum(VI * M[i,])
	Q
}

build.VA <- function (Ts, TS) {
	ns <- which(t>=Ts)[1]
	nS <- which(t>TS)[1]
	gate <- rep(0, length(t))
	gate[ns:nS] <- 1
	gate
}

## matrice de simulation (long)
Nsimu <- 1000
Ntot <- round(runif(Nsimu, 98, 102))
Mn <- matrix(0, nrow=Nsimu, ncol=length(simu.scinti (100, t, type="n")))
for (i in 1:Nsimu) Mn[i,] <- simu.scinti (Ntot[i], t, type="n")

Ntot <- round(runif(Nsimu, 98, 102))
Mg <- matrix(0, nrow=Nsimu, ncol=length(simu.scinti (100, t, type="g")))
for (i in 1:Nsimu) Mg[i,] <- simu.scinti (Ntot[i], t, type="g")

## essai de discri n/g
Qtot.gate <- build.VA (-5, 190)
Qana.gate <- build.VA (30, 190)

n.tot <- discri.ng (Mn,Qtot.gate)
n.ana <- discri.ng (Mn,Qana.gate)
g.tot <- discri.ng (Mg,Qtot.gate)
g.ana <- discri.ng (Mg,Qana.gate)

plot (c(0,20), c(0,.6), type = "n", xlab="TOT (A.U.)", ylab="FAST/TOT")
grid()
points (n.tot, n.ana/n.tot, col="red",pch=20)
points (g.tot, g.ana/g.tot, col="blue",pch=20)

##
## FOM 
##

n.proj <- n.ana[n.tot>8 & n.tot<10]/n.tot[n.tot>8 & n.tot<10]
g.proj <- g.ana[g.tot>8 & g.tot<10]/g.tot[g.tot>8 & g.tot<10]


n.h <- hist(n.proj, freq=FALSE, breaks=seq(0,0.6,by=0.02))
rug(n.proj,col="red")
lines (n.h$mids, dnorm(n.h$mids, mean(n.proj), sd(n.proj)), col="red")
g.h <- hist(g.proj, freq=FALSE, breaks=seq(0,0.6,by=0.02), same=TRUE)

plot (c(0,0.6), c(0,15), type="n", xlab="Qana/Qtot", ylab="density")
grid()
points (n.h$mids, n.h$density, col="red", pch="_")
points (g.h$mids, g.h$density, col="blue", pch="_")
lines (n.h$mids, dnorm(n.h$mids, mean(n.proj), sd(n.proj)), col="red")
lines (g.h$mids, dnorm(g.h$mids, mean(g.proj), sd(g.proj)), col="blue")
rug (n.proj, col="red")
rug (g.proj, col="blue")

(mean(n.proj)-mean(g.proj))/(sd(n.proj) + sd(g.proj))


##
## etude théorique
##

n.theo <- 0.5 * dexp(t, 1/9) + 0.5 * dexp(t, 1/60) 
g.theo <- 0.8 * dexp(t, 1/9) + 0.2 * dexp(t, 1/60) 
pu.nan <- signal.PU(t)
pu.theo <- c(pu.nan[!is.nan(pu.nan)], rep(0, which(!is.nan(pu.nan))[1]-1))
neut <- Re(fft(fft(n.theo) * fft(pu.theo), inverse=TRUE))
gamm <- Re(fft(fft(g.theo) * fft(pu.theo), inverse=TRUE))
neut <- neut / (sum(neut) * 0.1)
gamm <- gamm / (sum(gamm) * 0.1)

plot (c(0,200), c(0,max(gamm)), type="n", xlab="time (ns)", ylab="density probability")
grid ()
lines (t, neut, col="red", lwd=2)
lines (t, gamm, col="blue", lwd=2)

s.numb <- seq(-10, 200, by=2)
s.neut <- rep(0, length(s.numb))
s.gamm <- rep(0, length(s.numb))

resample <- function (pdf) {
	s <- rep(0, length(s.numb))
	for (i in 1:length(t)) s[6+t[i]/2] <- s[6+t[i]/2] + pdf[i]
	s		
}

s.neut <- resample (neut)
s.neut[1:5] <- 0
s.gamm <- resample (gamm)
s.gamm[1:5] <- 0
s.neut <- s.neut / sum(s.neut)
s.gamm <- s.gamm / sum(s.gamm)

plot (c(0,200), c(0,max(s.gamm)), type="n", xlab="time (ns)", ylab="density probability")
grid ()
points (s.numb, s.neut, col="red", lwd=3, pch=15)
points (s.numb, s.gamm, col="blue", lwd=3, pch=16)

V1 <- 8
Vrms <- 0
BL <- 0

signal <- function (A, T, S, N) {
	sum(A * S) / sum(T * S)
}

noise <- function (A, T, S, N) {
	num <- sum(A * S) * V1 * N 
	den <- sum(T * S) * V1 * N 
	v.quanta <- sum (((A*den-T*num)/den^2)^2 * V1^2 * S * 1.25 * N)
	v.noise <- sum (((A*den-T*num)/den^2)^2 * Vrms^2)
	v.baseline <- ((sum(A)*den-sum(T)*num)/den^2)^2 * 0.5 * BL^2

	sqrt(v.quanta+v.noise+v.baseline)
}

A <- c(rep(0,20), rep(1, 106-20))
T <- rep(1, length(A))

signal(A, T, s.neut, 100)
signal(A, T, s.gamm, 100)

noise (A, T, s.neut, 100)
noise (A, T, s.gamm, 100)

FOM <- function (A, T, N) {
	delta <- abs (signal(A, T, s.neut, N)-signal(A, T, s.gamm, N))
	delta / (noise (A, T, s.neut, N) + noise (A, T, s.gamm, N))
}

FOM (A, T, 100)

##
## FOM optimisation
##

cost <- function (params, ...) {
	-FOM(params, ...)
}

V1 <- 8
Vrms <- 0.45
BL <- 0

A <- c(rep(0,20), rep(1, 106-20))
T <- rep(1, length(A))
A_ <- nlm (cost, A, T, 100)

FOM (A_$estimate, T, 100)

plot (c(0,200), c(-0.5,1.2), type="n", xlab="time (ns)", ylab="ai coefficients")
grid ()
points(s.numb, A_$estimate)

##
## FOM vs N
##

V1 <- 8
Vmrs <- 0.45
BL <- 0

n <- c(10, 20, 50, 100, 200, 500, 1000)
plot (c(10, 1000), c(0, 7), type="n", log="x", xlab="number of PEs", ylab="FOM")
grid()
for (i in 1:length(n)) points (n[i], FOM (A_$estimate, T, n[i]), pch=16, col="blue")
n <- 10:1000
lines (n, 2*sqrt(n)/sqrt(100), lty=3)


##
## FOM vs Vrms 
##

N <- 100
V1 <- 8
Vrms<-0
BL<-0

plot (c(0, 2), c(1,2), type="n", xlab="electronics noise (mVrms)", ylab="FOM")
grid()
Vrms.list <- seq(0, 2, by=0.1)
FOM.list <- rep(0, length(Vrms.list))
FOMopt.list <- rep(0, length(Vrms.list))
for (i in 1:length(Vrms.list)) {
	Vrms <- Vrms.list[i]
	FOMopt.list[i] <- FOM(A_$estimate, T, N)
	FOM.list[i] <- FOM(A, T, N)
}
lines (Vrms.list, FOMopt.list, pch=16, lwd=3, col="red")
lines (Vrms.list, FOM.list, pch=16, lwd=3, col="green")


##
## FOM vs BL 
##

V1 <-8
N <- 100
Vrms<-0
BL<-0

plot (c(0, 2), c(1,2), type="n", xlab="baseline variations (mVpktopk)", ylab="FOM")
grid()
BL.list <- seq(0, 2, by=0.1)
FOM.list <- rep(0, length(BL.list))
FOMopt.list <- rep(0, length(BL.list))
for (i in 1:length(BL.list)) {
	BL<- BL.list[i]/2
	FOMopt.list[i] <- FOM(A_$estimate, T, N)
	FOM.list[i] <- FOM(A, T, N)
}
lines (BL.list, FOMopt.list, pch=16, lwd=3, col="red")
lines (BL.list, FOM.list, pch=16, lwd=3, col="green")


##
## FOM vs gain 
##

V1 <- 8
N <- 100
Vrms<-0.45
BL<-0

plot (c(0, 20), c(1,2), type="n", xlab="PE charge (mV.ns)", ylab="FOM")
grid()
V1.list <- seq(0, 20, by=0.1)
FOM.list <- rep(0, length(V1.list))
FOMopt.list <- rep(0, length(V1.list))
for (i in 1:length(V1.list)) {
	V1<- V1.list[i]
	FOMopt.list[i] <- FOM(A_$estimate, T, N)
	FOM.list[i] <- FOM(A, T, N)
}
lines (V1.list, FOMopt.list, pch=16, lwd=3, col="red", lty=3)
lines (V1.list, FOM.list, pch=16, lwd=3, col="green", lty=3)
V1 <- 8



##
## delayed gate
##

plot (c(0,100), c(0.5,1), type="n", xlab="time (ns)", ylab="FOM / FOMopt")
grid()
for (i in 10:90) {
	A <- c(rep(0,i), rep(1, 106-i))
	points(s.numb[i], FOM(A, T, 100)/FOM (A_$estimate, T, 100), pch=16, col="red")
}

A <- c(rep(0,21), rep(1, 106-21))
FOM(A, T, 100)/FOM (A_$estimate, T, 100)
plot (c(0,200), c(0,1), type="n", xlab="time (ns)", ylab="ai coefficients")
grid ()
points(s.numb, A, col="red")



##
## check
##

random.pulse <- function (s, N) {
	#rpois (length(s), s*N) * V1 + rnorm(length(s), 0, Vrms)
	support <- rpois(length(s), s * N)
	
	randomPU <- function (n) sum(rpois(n, 4)/4)
	sapply(support, randomPU) * V1 + rnorm(length(s), 0, Vrms) + BL*sin(10000*runif(1,0,1))
}

MC <- function (A, T, s, Nrange) {
	Nsimu <- 1000
	A.g <- rep(0, Nsimu)
	T.g <- A.g
	Q.g <- A.g
	for (i in 1:Nsimu) {
		n <- runif(1, Nrange[1], Nrange[2])
		signal <- random.pulse (s, n)
		A.g[i] <- sum(signal * A)
		T.g[i] <- sum(signal * T)
		Q.g[i] <- sum(signal) 
	}
	data.frame (A=A.g, T=T.g, Q=Q.g)
}

n.dist <- MC (A, T, s.neut, c(99,101))
g.dist <- MC (A, T, s.gamm, c(99,101))

n.dist <- MC (A_$estimate, T, s.neut, c(99,101))
g.dist <- MC (A_$estimate, T, s.gamm, c(99,101))

plot (c(0, max(g.dist$Q)), c(.5,1), type="n")
points(n.dist$Q, n.dist$A/n.dist$T, col="blue")
points(g.dist$Q, g.dist$A/g.dist$T, col="red")

(mean(n.dist$A/n.dist$T)-mean(g.dist$A/g.dist$T)) / (sd(n.dist$A/n.dist$T) + sd(g.dist$A/g.dist$T))

##
## figures
##

V1 <- 0
Vrms <- 0.45
BL <- 0

plot (c(0,200), c(-2,2), type = "n", xlab="sample", ylab="noise (mV)")
grid ()
points (s.numb, random.pulse (s.neut, 100), col="green")
points (s.numb, random.pulse (s.neut, 100), col="red")
points (s.numb, random.pulse (s.neut, 100), col="blue")

V1 <- 0
Vrms <- 0
BL <- 0.5

plot (c(0,200), c(-2,2), type = "n", xlab="sample", ylab="baseline (mV)")
grid ()
points (s.numb, random.pulse (s.neut, 100), col="green")
points (s.numb, random.pulse (s.neut, 100), col="red")
points (s.numb, random.pulse (s.neut, 100), col="blue")

V1 <- 8
Vrms <- 0
BL <- 0

plot (c(0,200), c(0,100), type = "n", xlab="sample", ylab="signal (mV)")
grid ()
points (s.numb, random.pulse (s.neut, 100), col="green")
points (s.numb, random.pulse (s.neut, 100), col="red")
points (s.numb, random.pulse (s.neut, 100), col="blue")