#---------------------------------------------------------------
# Inference of rates of signal and background Poisson processes
#
# GdA, March 2022
#----------------------------------------------------------------
# function pausa
pause <- function() {
  cat ("\n >> Give enter to continue\n")
  scan()
}

library(rjags)

model = "inf_r_bck_measured.bug"  # model file

data <- NULL   # R list containing data
data$X <- 1    # observed nr of counts from signal+background
data$T <- 1   # time of measurement signal+background
data$TB <- 10000   # time of measurement of background alone
data$XB <- 10000  # observed nr of counts from background alone

jm <- jags.model(model, data)  # define the model

update(jm, 100)    # "burn in": the chain runs but history not recorded, 
                   # just to get rid of initial positions
                   # (exaggerated in this case!)
chain <- coda.samples(jm, c("r","rB"), n.iter=10000)  # sampling 

print(summary(chain))
plot(chain)

pause()
# let us change the format of the chain into a 'data.frame'  
chain.df <- as.data.frame( as.mcmc(chain) )

layout(matrix(c(1:2),c(2,1)))
# plot a modo nostro
hist(chain.df$r, nc=100, prob=TRUE,
     col='yellow',
     xlab='r', ylab='f(r)', main='Inference on r, intensity of signal process')
hist(chain.df$rB, nc=100, prob=TRUE,
     col='red',
     xlab='rB', ylab='f(rB)', main='Inference rB, intensity of backd process')
layout(1)
      
pause()
plot(chain.df$rB, chain.df$r, xlab='r background', ylab='r signal', col='blue')