require (signal)
require (inline)

dt <- 0.01				# time in us
t <- seq(0, 10000, by=dt)	# simulation (10ms)

Cf <- 2e-12				# feedback capacitance (F)
Rf <- 100e6				# feedback resistor (Ohm)

Cd <- 10e-12      # detector capacitance (F)
Rp <- 100e6       # polarisation resistor (Ohm)
Irev <- 1e-9     # reverse current (A)

Sioa <- 1e-15    # //current noise gen of OpAmp (A/Hz^1/2)
Seoa <- 5e-9     # + voltage noise gen of OpAmp (V/Hz^1/2)
GBP <- 1e8        # gain bandwith product (Hz)
k <- 1.38e-23
T <- 300
e <- 1.6e-19

## s.current
## generates the detector current at random time
## counting.rate = counting rate (c/s)
## energy = particles energy (MeV) in Ge (w=3eV)
## variable.rise.time = TRUE -> random e/h 
s.current <- function(counting.rate, energy, variable.rise.time=FALSE) {

	print (paste ("average current/pulse =", energy * 1e6/3.6*1.6e-19 / (dt * 1e-6)))

	I <- rep(0, length(t))
    #intervalle de temps entre les particules
	t.pulse <- rexp (100000, counting.rate)*1e6
    #integration des temps
	t.pulse <- 10+cumsum(t.pulse)
	t.pulse <- t.pulse[t.pulse < t[length(t)]]
	if (variable.rise.time==FALSE) {
		I[t.pulse/dt] <- energy * 1e6/3.6*1.6e-19 / (dt * 1e-6)
	} 
	I
}

## s.CSP
## generates de voltage signal after Charge sensing preamplifier
s.CSP <- function (i) {
	my.CSP <- Arma(b = Rf, 
			a = c(1+Rf*Cf/(dt*1e-6), -Rf*Cf/(dt*1e-6)))
	filter (my.CSP, i)
}

n.CSP <- function () {
  my.par <- Arma (b = Rf,
      a = c(1+Rf*Cf/(dt*1e-6), -Rf*Cf/(dt*1e-6)))
  # parallel noise
  noise.par <- rnorm(length(t), 0, sqrt((2*e*Irev + 4*k*T/Rf + 4*k*T/Rp + Sioa^2) / (2 * dt * 1e-6)))  

  tau <- 1/(2 * pi * GBP) * Cd/Cf
  my.ser <- Arma (b= 1+Cd/Cf,
      a = c(1+tau/(dt*1e-6), -tau/(dt*1e-6))) 
  noise.ser <- rnorm(length(t), 0, sqrt(Seoa^2/ (2 * dt * 1e-6)))
  filter (my.par, noise.par) + filter (my.ser, noise.ser)
}

## s.fastout
## produces a "fast out" signal
## tau = time constant (?s)
s.fastout <- function (v.csp, tau) {
	my.HP <- Arma (b = c(tau / dt, -tau / dt),
			a = c(1+tau/dt, -tau/dt))
	dummy <- filter (my.HP, v.csp)
	my.LP <- Arma (b = 1,
			a = c(1+tau/dt, -tau/dt))
	filter (my.LP, dummy)
}

## s.sa
## produces th CR-RC4 signal without PZ compensation
## tau = time constant (?s)
s.SA.CRRC4 <- function (v.csp, tau) {
	my.HP <- Arma (b = c(tau / dt, -tau / dt),
			a = c(1+tau/dt, -tau/dt))
	dummy <- filter (my.HP, v.csp)
	my.LP <- Arma (b = 1,
			a = c(1+tau/dt, -tau/dt))
	dummy <- filter (my.LP, dummy)
	dummy <- filter (my.LP, dummy)
	dummy <- filter (my.LP, dummy)
	filter (my.LP, dummy)
} 

s.CPZ<-function(v.csp,tau.zero){
    my.CPZ<-Arma(b=c(1+dt/tau.zero,-1),a=c(1,-1))
    filter (my.CPZ, v.csp)
}

## s.trigger
## produces a trigger signal 0->1 = trigger
## s = signal
## trigger.level = trigger level
## time.interval = timebetween pulses (?s)

Ccode <- "
int i = 0;
int triggered = 0, ToT;

do {
	if (triggered == 0) {
        if (s[i] > *trigger_level && s[i] - s[i-1] > 0.0) {
			triggered = 1;
			ToT = i;
			trig[i] = 1;
		} else trig[i] = 0;
	}
	else {
		if (i - ToT > *time_interval)	triggered = 0;
		trig[i] = 0;
	}
	i++;
} while (i < *n);
"

C.trigger <- cfunction (signature (n="integer", s="numeric", 
	trigger_level="numeric", time_interval="integer", trig="numeric"), 
	Ccode, language="C", convention=".C")


s.trigger <- function (s, trigger.level, time.interval) {
	trigger <- rep(0, length(t))	
	S.out <- C.trigger (length(t), s, trigger.level, 
		round(time.interval / dt), trigger)
	S.out$trig
}

## oscilloscope
## s = signal
## trig = trigger
## t.before, t.after = time window (us)
## v.min, v.max = amplitude window (V)
## N = number of pulses to show (0 = all)

oscilloscope <- function (s, trig, t.before=1, t.after=30, v.min=-10, v.max=10, N=0) {
	
	if (N == 0) { 
		pulse.index <- which(trig == 1)
	}
	else {
		pulse.index <- which(trig == 1)[1:N]
	}

	plot (c(-t.before, t.after), c(v.min, v.max), type = "n",
		xlab = "time (us)", ylab="amplitude (V)")
	grid()

	t.scope <- seq(-t.before, t.after, by=dt)
	for (i in 1:length(pulse.index)) {
		start <- t.before / dt
		stop <- t.after / dt
		lines (t.scope, s[pulse.index[i]+ (-start):stop])
	}

}


## s.peak_hold
## extract maxima of the signal
## s = signal
## trigger.level = trigger level
## trig start -> >1.1 trigger.level
## trig stop -> <0.9 trigger.level

Ccode <- "
int i = 0;
int triggered = 0;
double smax=0.0;
*n_detect = 0;

do {
	if (triggered == 0) {
        	if (s[i] > *trigger_level * 1.1) {
			smax = s[i];
			triggered = 1;
		}
	}
	else {
		if (s[i] > smax) smax = s[i];
		if (s[i] < *trigger_level * 0.9) {
			triggered = 0;
			if (*n_detect < 10000) {
				peaks[*n_detect] = smax;
				(*n_detect)++;
			}
		}
	}
	i++;
} while (i < *n);


"

C.peak_hold <- cfunction (signature (n="integer", s="numeric", 
	trigger_level="numeric", n_detect="integer", peaks="numeric"), 
	Ccode, language="C", convention=".C")


s.peak_hold <- function (s, trigger.level) {
	peaks <- rep(0, 10000)
	n_detect <- 0	
	S.out <- C.peak_hold (length(t), s, trigger.level, 
		n_detect, peaks)
	S.out$peaks[1:S.out$n_detect]
}

## spectro
## produces the histogram of input data
##

spectro <- function (data) {
	b.min <- 0.0041
	b.max <- 0.00432
	b <- seq(b.min, b.max, by=(b.max - b.min)/100)	

	data.plot <- data[data>b.min & data<b.max]
	hist (data.plot, freq=FALSE, breaks=seq(b.min,b.max,length.out=100))
	lines (b, dnorm(b, mean(data.plot), sd(data.plot)), col="red")	

	print (paste("peak avg = ", mean(data.plot)))
	print (paste("peak std = ", sd(data.plot)))
	print (paste("res = ", sd(data.plot) / mean(data.plot) * 1000))
		
}