################################
library(neutRino)
library(inline)
library(Rcpp)
library(RcppArmadillo)
library(tikzDevice)
library(data.table)
library(ggplot2)
library(weights)
library(stats)
library(plyr)
###############################
Nu.ErecoAsEle <- function(Elepton, Plepton, cos.theta) {
Mp2 <- 1
Mn2 <- 1
me2 <- 1
Eb <- 27 #MeV
Ereco <- (Mp2 - (Mn - Eb)^2 - me2 + 2(Mn - Eb) * Elepton)
Ereco <- Ereco/2/(Mn - Eb - Elepton + Plepton * cos.theta)
return(Ereco)
}
Nu.ErecoAsMuon <- function(Elepton, Plepton, cos.theta) {
Mp2 <- 1
Mn2 <- 1
ml2 <- 1
Eb <- 27 # MeV
Ereco <- (Mp2 - (Mn - Eb)^2 - ml2 + 2(Mn - Eb) * Elepton)
Ereco <- Ereco/2/(Mn - Eb - Elepton + Plepton * cos.theta)
return(Ereco)
}
LogLike.Histo <- function(hobs, hexp) {
LL <- hobs - hexp - hobs * log(hobs/hexp)
LL[!is.finite(LL)] <- 0
LL <- sum(LL, na.rm = TRUE)
return(-2 * LL)
}
Chi2.Histo <- function(hobs, herr, hexp) {
chi2 <- 0.5 * sum((hobs - hexp)^2/(herr^2))
return(chi2)
}
####
set.seed(1234)
####
Stats <- function(x) {
Mean <- mean(x, na.rm = TRUE)
Median <- median(x, na.rm = TRUE)
SD <- sd(x, na.rm = TRUE)
Min <- min(x, na.rm = TRUE)
Max <- max(x, na.rm = TRUE)
xx <- x[x > Median]
Qr <- quantile(xx, 0.68)
xx <- x[x <= Median]
Ql <- quantile(xx, 1 - 0.68)
return(c(Mean = Mean, SD = SD, Min = Min, Max = Max, Median = Median, Ql = Ql, Qr = Qr))
}
stat_fluct <- function(n, opt = "Poisson") {
if (n < 0) {
return(0)
}
dn <- 0
if (opt == "Gaus")
dn <- rnorm(1, mean = 0, sd = sqrt(n)) else dn <- rpois(1, lambda = sqrt(n))
if (n + dn > 0)
return(n + dn) else return(0)
}
DataRead.AsFrame <- function(file = file.name, check.firstCol = T, ...) {
data <- read.delim(file = file.name, ...)
if (check.firstCol) {
# remove first colums if all NA's
nevents <- nrow(data)
col.1 <- data[, 1]
nna <- sum(is.na(col.1))
if (nna == nevents) {
data <- data[, 2:ncol(data)]
}
}
return(data)
}
####
file.name <- "data/nue_events_list.txt"
nue_data <- DataRead.AsFrame(file = file.name, check.firstCol = T, sep = " ", header = F)
nevents <- nrow(nue_data)
# add event osc prob (init at 1)
nue_data <- cbind(nue_data, rep(1, nevents))
# Names of event variables
var.names <- c("event", "iflav", "isel", "iint", "ifv", "Ereco", "Etrue", "pathLen", "weight", "IsNuE", "IsBkg", "IsGamma", "OscProb")
nvars <- length(var.names)
names(nue_data) <- var.names
## Ereco energy binning units: MeV
bins.ereco <- c(200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200,
2350, 2500, 2700, 3000, 3300, 3500, 4000, 4400, 5000, 6000, 10000)
bins.etrue <- seq(from = 0, to = 20000, length.out = 200)
bins.egamma <- c(200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1700, 1900, 2200, 2500, 2800, 4000,
10000)
##
labl.ereco <- "Ereco [MeV]"
labl.etrue <- "Etrue [MeV]"
labl.entries <- "Entries"
hErecoNuE <- wtd.hist(main = "NuE Selection", xlab = labl.ereco, ylab = labl.entries, x = nue_data[nue_data$IsNuE == 1, ]$Ereco,
weight = nue_data[nue_data$IsNuE == 1, ]$weight, freq = T, breaks = bins.ereco)
hErecoBkg <- wtd.hist(main = "Background", xlab = labl.ereco, ylab = labl.entries, x = nue_data[nue_data$IsBkg == 1, ]$Ereco, weight = nue_data[nue_data$IsBkg ==
1, ]$weight, freq = T, breaks = bins.ereco)
hErecoGamma <- wtd.hist(main = "Gamma Selection", xlab = labl.ereco, ylab = labl.entries, x = nue_data[nue_data$IsGamma == 1, ]$Ereco,
weight = nue_data[nue_data$IsGamma == 1, ]$weight, freq = T, breaks = bins.ereco)
hEtrueNuE <- wtd.hist(main = "NuE Selection", xlab = labl.etrue, ylab = labl.entries, x = nue_data[nue_data$IsNuE == 1, ]$Etrue,
weight = nue_data[nue_data$IsNuE == 1, ]$weight, freq = T, breaks = bins.etrue)
# apply stat fluc around the bin content value get the transpost in order ot have a column for each bin
opt <- "Gaus"
hist.fluct <- t(replicate(10000, sapply(hErecoNuE$counts, stat_fluct, opt = opt)))
hist.fluct.summary <- apply(hist.fluct, 2, Stats)
hist.fluct.Median <- apply(hist.fluct.summary, 2, function(x) as.list(x)$Median)
hist.fluct.Ql <- apply(hist.fluct.summary, 2, function(x) as.list(x)$Ql)
hist.fluct.Qr <- apply(hist.fluct.summary, 2, function(x) as.list(x)$Qr)
hist.fluct.summary
plot(hErecoNuE$mids, hist.fluct.Qr, type = "l", ylim = c(0, max(hist.fluct.Qr) * 1.1), ylab = "Entries", xlab = "Ereco [Mev]")
lines(hErecoNuE$mids, hist.fluct.Ql)
lines(hErecoNuE$mids, hist.fluct.Median, col = "red")
## Apply Oscillation
elist <- seq(0.1, 5, 0.01)
EneLen <- cbind(nue_data$Etrue/1000, nue_data$pathLen/1000)
pars <- c(sin2 = 0.7, dm2 = 50)
# pars <- c(sin2=1,dm2=2)
# baseline <- 0.280 # [km], nominal max
baseline <- 0.247 # [km], average
probs <- apply(EneLen, 1, function(x) Sterile.ProbNuEToNuE(x[1], x[2], pars = pars))
# probsLconst <- apply(EneLen,1,function(x) Sterile.ProbNuEToNuE(x[1],baseline,pars=pars))
probsLconst <- Sterile.ProbNuEToNuE(elist, baseline, pars = pars)
plot <- plot(EneLen[, 1], probs, type = "p", cex = 0.01, xlim = c(0, 5), ylab = "Probability", xlab = "True Energy [GeV]")
par(new = T)
plot <- plot(elist, probsLconst, type = "l", cex = 0.01, xlim = c(0, 5), ylab = "", xlab = "", main = "nue disappearance prob")
# Update Probs
# go here####
nue_data$OscProb <- probs
##
hErecoNuEOsc <- wtd.hist(main = "NuE Selection", xlab = labl.ereco, ylab = labl.entries, x = nue_data[nue_data$IsNuE == 1, ]$Ereco,
weight = nue_data[nue_data$IsNuE == 1, ]$OscProb * nue_data[nue_data$IsNuE == 1, ]$weight, freq = T, breaks = bins.ereco)
ratioEreco <- hErecoNuEOsc$counts/hErecoNuE$counts
hEtrueNuEOsc <- wtd.hist(main = "NuE Selection", xlab = labl.etrue, ylab = labl.entries, x = nue_data[nue_data$IsNuE == 1, ]$Etrue,
weight = nue_data[nue_data$IsNuE == 1, ]$OscProb * nue_data[nue_data$IsNuE == 1, ]$weight, freq = T, breaks = bins.etrue)
ratioEtrue <- hEtrueNuEOsc$counts/hEtrueNuE$counts
plot(hErecoNuE$mids, ratioEreco, type = "S")
plot(hEtrueNuE$mids, ratioEtrue, type = "S")
#
hObs.EtrueNuE <- hEtrueNuEOsc
base <- seq(0.21, 0.28, length.out = 200)
sin2 <- seq(0.5, 1, length.out = 100)
vobs <- hEtrueNuE$counts * Sterile.ProbNuEToNuE(hEtrueNuE$mids, 0.23, pars = pars)
nLL <- as.numeric()
for (i in 1:length(sin2)) {
p <- c(sin2[i], pars[2])
v <- hEtrueNuE$counts * Sterile.ProbNuEToNuE(hEtrueNuE$mids, 0.23, pars = p)
# nLL[i] <- LogLike.Histo(hObs.EtrueNuE$counts,v)
nLL[i] <- LogLike.Histo(vobs, v)
}
plot(sin2, nLL)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.