phy-4660/notes.R

103 lines
2.6 KiB
R
Raw Normal View History

2017-04-28 08:25:56 +00:00
plot(y,xlab="Channel",ylab="Count",lwd=4,pch="-")
lines(v[3]*exp(-1/2*(x-v[1])^2/v[2]^2),col=2,lwd=1.4)
abline(v = 7813.63,col=4,lwd=1)
axis(3,at=7813.63,"5.486 MeV")
dev.print(pdf,"../report/calibration.pdf")
Alpha Activity
plot(data_4_2,xlab="Energy [MeV]",ylab="Count [Alpha Particles]",pch=".")
minor.tick(nx=5,ny=2)
abline(v = 5.631,col=2,lwd=1)
axis(3,at=5.631,"5.6")
abline(v = 5.323,col=2,lwd=1)
axis(3,at=5.323,"5.3")
dev.print(pdf,"../report/activity.pdf")
For Gaussian Fits:
─────────────
x <- seq_along(y)
f <- function(par)
{
m <- par[1]
sd <- par[2]
k <- par[3]
rhat <- k * exp(-0.5 * ((x - m)/sd)^2)
sum((y - rhat)^2)
}
optim(c(15, 2, 1), f, method="BFGS", control=list(reltol=1e-9))
─────────────
I propose to use non-linear least squares for this analysis.
# First present the data in a data-frame
tab <- data.frame(x=seq_along(y), y=y)
#Apply function nls
(res <- nls( y ~ k*exp(-1/2*(x-mu)^2/sigma^2), start=c(mu=7500,sigma=50,k=1) , data = tab))
And from the output, I was able to obtain the following fitted "Gaussian curve":
v <- summary(res)$parameters[,"Estimate"]
plot(y~x, data=tab)
plot(function(x) v[3]*exp(-1/2*(x-v[1])^2/v[2]^2),col=2,add=T,xlim=range(tab$x) )
Mirror Gumbel Distribution:
(res <- nls( y ~ (1/beta)*exp((x-alpha)/beta - exp((x - alpha)/beta)), data=cal,start=c(alpha=7800,beta=1),control=list(minFactor=1/50000,maxiter=1000)))
v <- summary(res)$parameters[,"Estimate"]
plot(function(x) (1/v[2])*exp((x-v[1])/v[2] - exp((x - v[1])/v[2])),col=2,add=T,xlim=range(tab$x))
──────────────────────────────────────────────────────────────────────────
─────────────
File writing
write(x, file = "data",
ncolumns = if(is.character(x)) 1 else 5,
append = FALSE, sep = " ")
dev.print(pdf,'filename.pdf')
──────────────────────────────────────────────────────────────────────────
─────────────
Axis Formatting
# Add minor tick marks
library(Hmisc)
minor.tick(nx=n, ny=n, tick.ratio=n)
──────────────────────────────────────────────────────────────────────────
─────────────
Importing CSV