Nothing
#' @title Create empty spectrum
#'
#' @description \code{spectrum} return an
#' empty spectrum, i.e. a matrix with four
#' columns particle, kinetic energy (mid bin, in MeV.u),
#' bin width and number of particles in the mid
#' (absolute, not density)
#'
#' @param nrow Number of rows in the spectrum matrix
#'
#' @seealso \code{\link{dataSpectrum}}
#'
#' @family Data structures
#'
#' @author Steffen Greilich
spectrum <- function(nrow){
return(matrix(nrow = nrow,
ncol = 4,
dimnames = list(NULL,
c("particle.no",
"E.MeV.u",
"dE.MeV.u",
"N"))))
}
#' @title Class for particle spectra
#'
#' @description \code{dataSpectrum} is the main class
#' to hold particle spectra information, i.e. the number of
#' particles differential in charge and energy
#'
#' @slot depth.g.cm2 Depth at which the spectrum is located in g/cm2
#' @slot spectrum \code{\link{spectrum}} (i.e. matrix) structure containing the actual particle numbers
#'
#' @family Data structures
#'
#' @author Steffen Greilich
setClass( Class = "dataSpectrum",
slots = c( depth.g.cm2 = "numeric",
spectrum = "matrix"),
prototype = list( depth.g.cm2 = numeric(),
spectrum = spectrum(0)) )
#' @title Constructor for \code{\link{dataSpectrum}} class
#'
#' @description Constructs a dataSpectrum object from a given dataSPC object,
#' returning the spectrum at a specific depth
#'
#' @details Interpolation is done by linear mixture of the spectra which can lead
#' to funny shapes such as double peaks. This is, however, the only way without
#' using explicit transport calculation and also employed for example in TRiP98.
#'
#' @param SPC.data dataSPC object holding spectrum information for multiple sample depths
#' @param depth.g.cm2 Depth at which the spectrum is located in g/cm2. This can be any
#' depth, spectra will be interpolated between the adjacent sample depths.
#'
#' @family Data structures
#'
#' @author Steffen Greilich
dataSpectrum <- function(SPC.data, depth.g.cm2){
res <- lapply(1:length(depth.g.cm2),
function(i, s, d){
new("dataSpectrum",
depth.g.cm2 = d[i],
spectrum = spectrum.at.depth.g.cm2(s@spectra, d[i]))},
s = SPC.data,
d = depth.g.cm2)
}
spectrum.get.depth.g.cm2 <- function(spectrum){
return(spectrum@depth.g.cm2)
}
#' R function for interpolation of spc files
#' with depth, was in libamtrack and moved to HITXML Jul15, sgre
spectrum.at.depth.g.cm2 <- function(spectra, depth.g.cm2, interpolate = TRUE)
{
depth.step.of.spc <- sort(unique(spectra[,"depth.step"]))
depth.g.cm2.of.spc <- sort(unique(spectra[,"depth.g.cm2"]))
depth.step.interp <- approx( x = depth.g.cm2.of.spc,
y = depth.step.of.spc,
xout = depth.g.cm2,
rule = 2)$y
depth.step.int <- floor(depth.step.interp)
depth.step.frac <- depth.step.interp - depth.step.int
spectrum.upstream <- spectra[spectra[,"depth.step"] == depth.step.int,]
spectrum.interp <- spectrum(nrow(spectrum.upstream))
spectrum.interp[,"E.MeV.u"] <- spectrum.upstream[,"E.MeV.u"]
spectrum.interp[,"dE.MeV.u"] <- spectrum.upstream[,"dE.MeV.u"]
spectrum.interp[,"particle.no"] <- spectrum.upstream[,"particle.no"]
if(interpolate & depth.step.frac != 0){
spectrum.downstream <- spectra[spectra[,"depth.step"] == (depth.step.int+1),]
spectrum.interp[,"N"] <- (1 - depth.step.frac) * spectrum.upstream[,"N.per.primary"] + depth.step.frac * spectrum.downstream[,"N.per.primary"]
}else{
spectrum.interp[,"N"] <- spectrum.upstream[,"N.per.primary"]
}
return(spectrum.interp)
}
#' @title Total number of particles in a spectrum
#'
#' @description Return sum of particles
#'
#' @param x Object of class \code{\link{dataSpectrum}} or list of objects of this class
#' @param particle.no if given, only these particles will be counted
spectrum.total.n.particles <- function(x, particle.no = NULL){
if(is.null(particle.no)){
FUN <- function(xx){sum(xx@spectrum[,"N"])}
}else{
FUN <- function(xx){
ii <- xx@spectrum[,"particle.no"] == particle.no
return(sum(xx@spectrum[ii,"N"]))
}
}
if(class(x) == "dataSpectrum"){
return(FUN(x))
}else{
return(sapply(x, FUN))
}
}
#' @title Fluence weighted LET of particles in a spectrum (in water, using ICRU)
#'
#' @description Return fLET of particles
#'
#' @param x Object of class \code{\link{dataSpectrum}} or list of objects of this class
#' @param particle.no if given, only these particles will be counted
spectrum.fLET <- function(x, stopping.power.source.no = 3, material.no = 1, particle.no = NULL){
if(is.null(particle.no)){
ii <- rep(TRUE, nrow(x@spectrum))
}else{
ii <- x@spectrum[,"particle.no"] == particle.no
}
x@spectrum <- x@spectrum[ii,]
FUN <- function(xx){
AT.fluence.weighted.LET.MeV.cm2.g(E.MeV.u = xx@spectrum[,"E.MeV.u"],
particle.no = xx@spectrum[,"particle.no"],
fluence.cm2 = xx@spectrum[,"N"],
material.no = 1,
stopping.power.source.no = 3)$returnValue}
if(class(x) == "dataSpectrum"){
return(FUN(x))
}else{
return(sapply(x, FUN))
}
}
#' @title Dose weighted LET of particles in a spectrum (in water, using ICRU)
#'
#' @description Return dLET of particles
#'
#' @param x Object of class \code{\link{dataSpectrum}} or list of objects of this class
#' @param particle.no if given, only these particles will be counted
spectrum.dLET <- function(x, particle.no = NULL){
if(is.null(particle.no)){
ii <- rep(TRUE, nrow(x@spectrum))
}else{
ii <- x@spectrum[,"particle.no"] == particle.no
}
x@spectrum <- x@spectrum[ii,]
FUN <- function(xx){
AT.dose.weighted.LET.MeV.cm2.g(E.MeV.u = xx@spectrum[,"E.MeV.u"],
particle.no = xx@spectrum[,"particle.no"],
fluence.cm2 = xx@spectrum[,"N"],
material.no = 1,
stopping.power.source.no = 3)$returnValue}
if(class(x) == "dataSpectrum"){
return(FUN(x))
}else{
return(sapply(x, FUN))
}
}
#' @title Mass stopping power for all energy bins in a spectrum
#'
#' @description (List of) vector(s) of mass stopping power values (in the order of particle number and energies)
#' of the given spectrum/a.
#'
#' @details Uses \code{\link[libamtrack]{libamtrack}} stopping power function
#'
#' @param x Object of class \code{\link{dataSpectrum}} or list of objects of this class
#' @param stopping.power.source Descriptor for source of stopping power data (\code{\link[libamtrack]{stopping.power.source}})
#' @param target.material Descriptor for target material (\code{\link[libamtrack]{material.no}})
spectrum.Mass.Stopping.Power.MeV.cm2.g <- function(x, stopping.power.source, target.material){
FUN <- function(xx, s, t){ AT.Mass.Stopping.Power( s,
xx@spectrum[,"E.MeV.u"],
xx@spectrum[,"particle.no"],
AT.material.no.from.material.name(t))$stopping.power.MeV.cm2.g}
if(class(x) == "dataSpectrum"){
return( FUN(x,
s = stopping.power.source,
t = target.material))
}else{ return( lapply(x,
FUN,
s = stopping.power.source,
t = target.material))
}
}
#' @title Dose for a spectrum
#'
#' @description (Vector of) dose in Gy
#'
#' @details Uses \code{\link[libamtrack]{libamtrack}} stopping power function
#'
#' @param x Object of class \code{\link{dataSpectrum}} or list of objects of this class
#' @param stopping.power.source Descriptor for source of stopping power data (\code{\link[libamtrack]{stopping.power.source}})
#' @param target.material Descriptor for target material (\code{\link[libamtrack]{material.no}})
spectrum.dose.Gy <- function(x, stopping.power.source, target.material){
FUN <- function(xx, s, t, d){ m <- spectrum.Mass.Stopping.Power.MeV.cm2.g(xx, s, t)
sum(m * xx@spectrum[,"N"]) / d * 1.60217657e-10}
density.g.cm3 <- AT.get.materials.data(AT.material.no.from.material.name(target.material))$density.g.cm3
if(class(x) == "dataSpectrum"){
return(FUN(x,
s = stopping.power.source,
t = target.material,
d = density.g.cm3))
}else{
return(sapply(x,
FUN,
s = stopping.power.source,
t = target.material,
d = density.g.cm3))
}
}
#' Method plot
setMethod(f = "plot",
signature = c("dataSpectrum"),
definition = function(x) {
Z <- AT.Z.from.particle.no(x@spectrum[,"particle.no"])$Z
lattice::xyplot(x@spectrum[,"N"] ~ x@spectrum[,"E.MeV.u"],
type = "S",
grid = TRUE,
groups = Z,
ylab = "particles / Gy",
scales = list(y = list(log = 10)),
auto.key = list(title = "Z",
space = "top",
columns = max(Z),
lines = TRUE,
points = FALSE))
})
setMethod(f = "+",
signature = c("dataSpectrum", "dataSpectrum"),
definition = function(e1, e2) {
new.spectrum <- e1@spectrum
new.spectrum[,"N"] <- e1@spectrum[,"N"] + e2@spectrum[,"N"]
return(new("dataSpectrum",
depth.g.cm2 = e1@depth.g.cm2,
spectrum = new.spectrum))
})
setMethod(f = "*",
signature = c("dataSpectrum", "numeric"),
definition = function(e1, e2) {
new.spectrum <- e1@spectrum
new.spectrum[,"N"] <- new.spectrum[,"N"] * e2
new("dataSpectrum",
depth.g.cm2 = e1@depth.g.cm2,
spectrum = new.spectrum)
})
#' @title Converts dataSpectrum to data table
#'
#' @description Converts dataSpectrum to data table
#'
#' @param x Object of class \code{\link{dataSpectrum}} or list of objects of this class
get.data.table <- function(x){
FUN <- function(xx){
cbind(data.table(depth.g.cm2 = rep(xx@depth.g.cm2, nrow(xx@spectrum))),
as.data.table(xx@spectrum))}
if(class(x) == "dataSpectrum"){
return(FUN(x))
}else{
return(rbindlist(lapply(x,
FUN)))
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.