Creating GamBet SRC files with R

(This is ninth article in the series Using R for Gambet statistical analysis.[1]) The focus of this article is to learn how to generate particle distributions with R and how to export them in a form that can be used directly as input to GamBet. A second topic is how to run R scripts directly from the Windows command prompt or from batch files. With this capability, you can write scripts for automatic analyses of output from GamBet and other technical programs.

The parameters of primary particles for Monte Carlo simulations in GamBet are specified in source (SRC) files. Output escape files have the same format, so the output from one GamBet calculation can be used as the input to a subsequent one. For initial GamBet simulations, a source file representing specific distributions in position, velocity and energy in usually prepared. Although the GenDist utility can create several useful distributions, the possibilities with R are greatly expanded.

As a test case, we’ll generate an input electron beam for a 3D calculation with current CurrTotal = 2.5 A and average energy E0 = 20.0 keV. The beam has a circular cross section with a Gaussian distribution in x and y of width Xw = Yw = 1.4 mm. The average position is [X0,Y0,Z0] = [0.0 mm, 0.0 mm, 0.0 mm). The parallel electrons move in z. There is a Gaussian energy spread of Ew = 500 eV.

To begin, we set up a script that can be tested in the interactive environment of RStudio and then see how to convert it to an autonomous program that can be run from a batch file. The first set of commands clears the workspace, sets a working directory and defines parameters:

rm(list=objects())
WorkDir = "C:/USINGRFORGAMBETSTATISTICALANALYSIS/Examples/Section09/"
setwd(WorkDir)
CurrTotal = 2.5
X0 = 0.0
Y0 = 0.0
Z0 = 0.0
Xw = 1.4
Yw = 1.4
E0 = 20000.0
Ew = 500.0
Ux0 = 0.0
Uy0 = 0.0
Uz0 = 1.0
NPart= 100000

The previous article discussed the SRC file format. The strategy will be to set up a vector of length NPart for each quantity in the file data lines, combine them into a data frame and then use the write.table() command to create the file. These commands create the vector for the Type column:

Type = character(length=NPart)
for (n in 1:NPart) {
Type[n] = "E"
}

The first command creates a character vector called Type with NPart blank entries. The loop fills the vector with the character E.

The quantities Z, Ux, Uy, and Uz are numerical vectors of length NPart that contain identical values. For this, it is convenient to create and to fill the vectors with the seq() and c() commands:

Z = c(seq(from=Z0,to=Z0,length.out=NPart))
Ux = c(seq(from=Ux0,to=Ux0,length.out=NPart))
Uy = c(seq(from=Uy0,to=Uy0,length.out=NPart))
Uz = c(seq(from=Uz0,to=Uz0,length.out=NPart))

The position vectors X and Y and the Energy vector are created using the rnorm() function. It generates NPart values following a specified normal distribution:

X = rnorm(NPart,mean=X0,sd=Xw)
Y = rnorm(NPart,mean=Y0,sd=Yw)
Energy = rnorm(NPart,mean=E0,sd=Ew)

The set of vectors is assembled into a data frame. Note that the column names are the same as the vector names, Type, Energy, X, …:

SRCRaw = data.frame(Type,Energy,X,Y,Z,Ux,Uy,Uz)

We must take some precautions. With the normal distribution, there is a very small but non-zero probability of extreme values of the position and energy. They could result in electrons outside the GamBet solution volume or negative energy values. Both conditions would lead to a program error. We create a subset of the raw data such that no electron has r > 7.0 mm or Energy ? 0.0:

SRCFile = subset(SRCRaw,((X^2+Y^2)<=25.0*Xw^2) & (Energy > 0.0))

We need to add the current per electron to complete the SRCFile data frame. Note the use of NLength, the number of electrons in the modified data frame, which may be less than NPart:

NLength = length(SRCFile$Type)
dCurr = CurrTotal/NLength

A column vector of identical current values is constructed and appended to the SRCFile data frame with the cbind() command:

Curr = c(seq(from=dCurr,to=dCurr,length.out=NLength))
SRCFile = cbind(SRCFile,Curr)

We also define two plotting vectors for latter use:

RVector = sqrt(SRCFile$X^2 + SRCFile$Y^2)
EVector = SRCFile$Energy

We’re ready to write the file. Some variables are set up in preparation:

FNameOut = "TestSRCGeneration.SRC"
HLine1 = "* GamBet Particle Escape File (Field Precision)"
HLine2 = paste("* NPart:",NLength)
HLine3 = "* Type   Energy         X            Y            Z             Ux           Uy           Uz       Curr"
HLine4 = "* =========================================================================================================="

Note that the actual number of electrons was written to HLine2. The following commands open the file (over-writing any previous version) and write the header using the cat() command:

cat(HLine1,file=FNameOut,append=FALSE,fill=TRUE)
cat(HLine2,file=FNameOut,append=TRUE,fill=TRUE)
cat(HLine3,file=FNameOut,append=TRUE,fill=TRUE)
cat(HLine4,file=FNameOut,append=TRUE,fill=TRUE)

We could add the table in its current form, but it would be nice to have the fixed-width format illustrated in the previous article. This command adds three initial spaces to the Type column:

SRCFile$Type = paste("   ",SRCFile$Type,sep="")

These commands convert the numbers in the columns to character representations with width 12 in either scientific or standard notation:

SRCFile$Energy = format(SRCFile$Energy,scientific=TRUE,digits=5,width=12)
SRCFile$X = format(SRCFile$X,scientific=TRUE,digits=5,width=12)
SRCFile$Y = format(SRCFile$Y,scientific=TRUE,digits=5,width=12)
SRCFile$Z = format(SRCFile$Z,scientific=TRUE,digits=5,width=12)
SRCFile$Ux = format(SRCFile$Ux,scientific=FALSE,digits=6,width=12)
SRCFile$Uy = format(SRCFile$Uy,scientific=FALSE,digits=6,width=12)
SRCFile$Uz = format(SRCFile$Uz,scientific=FALSE,digits=6,width=12)
SRCFile$Curr = format(SRCFile$Curr,scientific=TRUE,digits=6,width=12)

Note that RVector and EVector were defined when the data entries were still numbers. Finally, this command writes the SRC file:

write.table(SRCFile,file=FNameOut,sep=" ",append=TRUE,col.names=FALSE,quote=FALSE,row.names=FALSE)

Here is a sample of the result:

* GamBet Particle Escape File (Field Precision)
* NPart: 100000
* Type   Energy         X            Y            Z             Ux           Uy           Uz       Curr
* ==========================================================================================================
E   2.0448e+04  -6.8232e-01   8.4354e-01        0e+00            0            0            1      2.5e-05
E   1.9949e+04  -4.8475e-02   1.5114e+00        0e+00            0            0            1      2.5e-05
E   2.1422e+04  -6.9644e-01  -4.9890e-01        0e+00            0            0            1      2.5e-05
E   1.9291e+04  -1.9870e+00   5.5780e-01        0e+00            0            0            1      2.5e-05
E   2.0032e+04  -5.6534e-01  -2.0669e-01        0e+00            0            0            1      2.5e-05
E   2.0752e+04   1.2376e+00   7.0758e-01        0e+00            0            0            1      2.5e-05
...

The format isn’t precisely the way we would like. For example, despite the setting digits=6, the entries in Ux are 0 rather than 0.00000. This is a quirk of R. The program will not write more significant figures than the defining values, no matter what you tell it. Perhaps there is a solution, but I haven’t found it. GamBet and GenDist will recognize the above form. Figure 1 shows a GenDist scatter plot of model particles in the x-y plane.

Scatter plot in x-y produced by GenDist

Figure 1. Scatter plot in x-y produced from the SRC file by GenDist

The following commands create histograms to display the distribution of electrons in radius and energy (Figure 2):

hist(RVector,breaks=100)
hist(EVector,breaks=100)
Variation of particle density with radius and the energy spectrum

Figure 2. Variation of particle density with radius and the energy spectrum created with the hist() command.

To conclude, we’ll modify the script so it may be called from a batch file. Suppose we want to make a series of GamBet runs with beams of differing width. For each run, we regenerate the input SRC file with R to represent the new width, run GamBet and rename critical output files so they are not over-written. Here’s a sample of a Windows batch file showing the first operation

START /WAIT RScript Sect09DemoScript.R "C:\USINGRFORGAMBETSTATISTICALANALYSIS\Examples\Section09\" "1.4"
START /WAIT C:\fieldp_basic\gambet\gambet BatchDemo
RENAME BatchDemo.GLS BatchDemo01.GLS
RENAME BatchDemo.G3D BatchDemo01.G3D
...

The /WAIT option in the START command ensures that the SRC file is available before starting GamBet and that there will not be a file access conflict when GamBet is running. Consider the first command. RScript is a special form of R to run scripts. To avoid typing the full path to program every time, I modified the Windows path to include C:\Program Files\R\R-3.2.0\bin\. The next command-line parameter is the name of the script followed by two character values that are passed to the script. The first is the working directory for the data and second is the value of Xw.

Here’s how the script we have been discussing is modified to act as an autonomous program. The first set of commands becomes:

rm(list=objects())
args = commandArgs()
WorkDir = args[6]
setwd(WorkDir)

The second command stores the command line parameters in a character vector args. You may wonder, why start with args[6]? A listing of args gives the following information:

[1] C:\PROGRA~1\R\R-32~1.0\bin\i386\Rterm.exe
[2] --slave
[3] --no-restore
[4] --file=Sect09DemoScript.R
[5] --args
[6] C:\USINGRFORGAMBETSTATISTICALANALYSIS\Examples\Section09\
[7] 1.4

The first argument is the name of the running program, a typical convention in Windows. The next four are the values of options set in RTerm. The values of interest start at the sixth component. The next change is in the initialization section:

Xw = as.numeric(args[7])

Because command line parameters are strings, we force the value to be interpreted as a number to define Xw. Finally, suppose we want to inspect the histograms, even though the program is running in the background. The ending statements are modified to:

pdf(file="Test.pdf")
hist(RVector,breaks=100)
hist(EVector,breaks=100)
dev.off()

The results of all plot commands between pdf() and dev.off() are sent to the specified PDF file where they can be viewed after the run.

In conclusion, the R package has a vast set of available commands and options that could occupy several textbooks. In this set of short articles, I’ve tried to cover the fundamental core. One of my goals has been to clarify the sometimes arcane syntax of R so you have the background to explore additional functions. A book form of the articles with an index will be available soon on our website.

Footnotes

[1] The entire series is available as a PDF book at http://www.fieldp.com/rintroduction.html.

[2] Contact us : techinfo@fieldp.com.

[3] Field Precision home page: www.fieldp.com.

Comments are closed.