mirror of
https://asciireactor.com/otho/phy-4660.git
synced 2024-11-25 06:05:09 +00:00
103 lines
2.6 KiB
R
103 lines
2.6 KiB
R
|
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
|
||
|
|
|