pause <- function() { cat ("\n >> give Enter to continue\n"); scan() }

fun <- function(x, c) {
    f <- rep(0,length(x))
    for(i in 1:length(c))  f <- f + c[i]*x^(i-1)
    return(f)
}

c <- c(10,10,-0.75)

# tanto per dare un'occhiata alla funzione 
# x <- seq(0,10,len=11)
# y <- fun(x,c)
#plot(x,y, ty='l')
#pause()

# definiamo i bin
x.min <- 0
x.max <- 10
n.bin <- 10
Dx <- (x.max -x.min)/n.bin
x.c <- seq(x.min+Dx/2, x.max-Dx/2, len=n.bin)   # bin centers
p <- fun(x.c, c)*Dx
p <- p/sum(p)

# generiamo gli eventi
N <- 1000
n.evts <- rmultinom(1, N, p)[,1]
plot(x.c, n.evts, pch=19, col='blue',
     xlim=c(x.min, x.max), ylim=c(0,max(n.evts)*1.05),
     xaxs = "i", yaxs = "i")
grid()
for(i in 1:length(x.c)){
    lines(c(x.c[i]-Dx/2, x.c[i]+Dx/2), rep(n.evts[i],2), col='blue', cex=1)
}
points(x.c,p*N, ty='l', col='cyan')
#x <- seq(x.min,x.max,len=101)
#points(x, fun(x,c)*N, ty='l', col='cyan')