Revised December 17
, 1997
THE MCLIB LIBRARY:
MONTE
CARLO SIMULATION OF NEUTRON SCATTERING INSTRUMENTS
Philip A. Seeger
239 Loma del Escolar, Los Alamos, NM 87544, USA
E

mail:
PASeeger@aol
.
com
ABSTRACT
This report describes the philosophy and structure of MCLIB, a Fortran library of Monte Carlo
subroutines which has been developed to test designs of neutron scattering instruments. Reports
from the proceedings of two meetings [1,2] are incorporated. A p
air of programs (
LQDGEOM
and
MC_RUN
) which use the library is shown as an example; a second example using a focusing mirror
to increase intensity in a small

angle instrument is also shown. Subroutine abstracts are
included in an appendix; the latest relea
se of the source code and documentation may be
obtained by anonymous ftp from
ftp://azoth.lansce.lanl.gov/pub/mclib/
. Work is also continuing
on a more friendly web

based user interface, currently located at
http://bayberry.lanl.gov/lansce/.
User input i
s requested for additional features to be added to the library.
1. Introduction
Monte Carlo is a method to integrate over a large number of variables. Random numbers are
used to select a value for each variable, and the integrand is evaluated. The pro
cess is repeated
a large number of times and the resulting values are averaged. For a neutron transport problem,
we first select a neutron from the source distribution, and project it through the instrument using
either deterministic or probabilistic algo
rithms to describe its interaction whenever it hits
something. If it hits a detector, we tally it in a histogram representing where and when it was
detected. This is intended to simulate the process of running an actual experiment (but it is
much
slower)
. Monte Carlo is a useful supplement to analytical treatment of an instrument, in
particular to check and demonstrate “non

intuitive” focusing arrangements, but should never be
used as a substitute for thinking. (I am grateful to Jack Carpenter for remin
ding us of this
limitation of Monte Carlo.)
The present MCLIB library has been derived from codes written by Mike Johnson at the
Rutherford Laboratory [3]. Significant additions and revisions were made by this author in
1984 [4], and the entire code was
rewritten in a structured form in 1994 [5]. Whenever the code
has been applied to new problems, additions have been made. Several applications of the code
have been presented at instrument design workshops [6]. Though not covered in this document,
the u
ser interface being developed in collaboration with Thierry Thelliez and Luke Daemen [7]
will have a major impact on the usefulness of the package. At this date, a prototype of the
interface may be visited at
http://bayberry.lanl.gov/lansce/Welcome.html
.
This will eventually
allow users to design instruments by defining the locations and properties of a variety of beam
elements.
The process is carried out in two stages. First, either the web page or a specific program must
be used to generate a descript
ion of the geometry of the instrument being simulated. For
2
MCLIB
example, the program
LQDGEOM
may be used to define a small

angle scattering instrument with
pinhole collimation, up to three choppers, and an on

axis 2

dimensional position sensitive
detector. Es
sentially all of the user interactions occur in this stage. The output is a geometry
file containing the complete problem definition. That file is then passed to a second

stage
program, for example
MC_RUN
which can be downloaded from
ftp://azoth.lansce.l
anl.gov/pub/mclib/
along with the Monte Carlo library and documentation),
which transports neutrons and tallies results in histograms. Principal outputs are a file with a
statistical summary, and a data file with histograms of the spectrum and detector.
A third stage,
which is not part of the Monte Carlo process, is to perform whatever data reduction is
appropriate to the experiment being simulated. Some measure of the information content (or a
“figure of merit”) is then used to evaluate the design of th
e instrument. Sample programs which
read and plot the output files are provided.
Features of MCLIB that are different from other Monte Carlo libraries include
Simplified transmission through materials. Rather than compute microscopic interaction
in a
simple (amorphous unpolarized) region, attenuation of the transmitted neutron is
calculated.
Optics at surfaces. When a neutron reaches a surface, the (complex) index of refraction is
computed to decide whether the neutron will reflect or refract.
Tim
e

dependent devices. There are region types to describe moving devices such as
choppers or a gravity focuser.
Acceleration of gravity is included in transport.
Scattering functions. Each kind of scattering sample is a region type. The scattering
algorithm may be deterministic (reflectometry), probabilistic (hard

sphere scatterer), or a
combination (Bragg reflection into a Debye

Scherrer cone).
2. Geometry Des
cription by Surfaces and Regions
The geometry of a system is described by surfaces and regions. A
surface
is defined by a
general 3

dimensional quadratic equation of the form
A
x
2
+
B
x +
C
y
2
+
D
y +
E
z
2
+
F
z +
G
+
P
xy +
Q
yz
+
R
zx
= 0
(1)
with 10 coefficients, plus a roughness (figure

error) parameter (
BETA
). The surface divides
3

dimensional space into two parts, which are called the + and
sides of the surface depending
on whether the left

hand side of eq. (1) evaluates to a positive or a negative value. For example,
a plane perpendicular to the z

axis at z = 1 can be expressed by the equation
z
1 = 0 ,
i.e.,
F
= 1 and
G
=
1 (a
ll other coefficients zero). Then all points with z < 1 are on the
side
and all points with z > 1 are on the + side of the surface. Higher

order surfaces such as toroids,
which can not be described by eq. (1), must instead be defined as parameters of a
special region
(
e.g.,
toroidal mirror, type 14). The scaling of eq. (1) is arbitrary, but we tend to evaluate
non

quadratic surfaces as m (coefficients
B
,
D
, and
F
dimensionless and
G
in m) and quadratic
surfaces as m
2
(coefficients
A
,
C
,
E
,
P
,
Q
, and
R
dimensionless, coefficients
B
,
D
, and
F
in m,
and
G
in m
2
). The parameter
BETA
is the length of a randomly oriented 3

dimensional vector
that is added to the unit vector normal to the mathematical surface to determine the surface
orientation when a particle interacts. For a perfect smooth surface,
BETA
= 0; for 0 <
BETA
< 1,
MCLIB
3
BETA
i
s the sine of the maximum angular deviation of the surface normal from smooth. If
BETA
< 0 (or
BETA
>> 1), the surface is completely random.
Note that the coordinate system being used is left handed: the instrument axis is the positive
z

direction, the x

axis is horizontal and positive to the right, and the y

axis is vertical with the
acceleration of gravity in the negative y

direction.
The geometric shape of each
region
is defined by its relationship to all of the defined surfaces.
A posit
ive or negative integer is placed in the region definition if every point in the region is on
the + or

side of the corresponding surface, and surfaces which do not bound the region are set
to zero. Special characteristics of the boundary are given by th
e value of the integer: ±1 for an
ordinary surface with roughness
BETA
and the possibility of refraction or critical reflection; ±2
for total reflection; ±3 for diffuse scattering (independent of
BETA
); ±4 for total absorption; ±5
for cases requiring special action (such as a coordinate transformation) whenever a particle
enters or leaves the region; and ±6 if the neutron history is to be split after crossing the surface.
Generally, no surface may be
used as a boundary of a region if any part of that surface is inside
the region, because then some points in the region would be on the + side and some on the
side
of that surface. Concave (reentrant) shapes are allowed when using quadratic surfaces as
in
fig.
1a, but the shape in fig. 1b requires that two regions be defined. Because definitions of this
sort become very complex when backscattering angles and multiple detectors are used, a special
form is provided for embedding one region in another. Fo
r example, the enclosing region
(“scattering chamber”) may be a large sphere. A detector embedded in the scattering chamber is
defined in the usual manner, but every one of its surfaces is also made a surface of the scattering
chamber with opposite sign a
nd with 10 added to the surface type. (Exceptions are if the
surface is already an actual boundary of the scattering chamber, or has already been used with
the opposite sign as now required. In the latter case it is necessary to define separate copies of
Figure 1.
Concave regions.
Equations of the numb
ered surfaces are
1:
x
=
0
2:
x
4
=
0
3:
y
=
0
4:
y
3
=
0
5:
x
2
6 x
+ y
2
6 y
+
14
=
0
or (x
3)
2
+ ( y
3)
2
= r
2
= 4
6:
x
1.5
=
0
7:
y
1
=
0
The shaded area in a) can be expressed unequivocally as a single region:
+1
1
+1
1
+1
0
0
To avoid ambiguity, the area in b) requires two regions, divided as shown by surface 7 (or alternatively
by surface 6).
+1
1
+1
0
0
0
1
+1
0
0
1
0
1
+1
The u
nshaded area in either case requires 5 regions for a complete unambiguous definition. (One
choice of region boundaries is shown by the faint lines.)
4
MCLIB
the surface for the two embedded regions that use it.) The presence of the 10s digit tells the
code to consider that surface as a
possible
exit from the enclosing region, but a particle will
remain within the enclosing region if no valid embedded region
is found on the other side.
3. Region Types
An
element
of the instrument consists of one or more
regions
. Every region has a
NAME
and a
pointer (
INDEX
). For a void drift region,
INDEX
= 0, and for every region which is not a void
INDEX
points to a loca
tion in a
REAL*4
parameter block which contains the region type number,
possibly followed by additional parameters. Future development of the library should be
accomplished by defining new region types and implementing the corresponding algorithms for
how
a neutron interacts in such regions (see section 9 below). Defined type numbers are listed
here, and the definitions of parameters are given in Appendix A. Types that have been added
since the ICANS 13 proceedings [1] are marked “
*
”.
Simple material typ
es:
type 0 =
total absorber; no parameters
type 1 =
amorphous unpolarized material; 4 parameters
type 2 =
aluminum, including Bragg edges; no parameters
type 3 =
hydrogenous, including multiple scattering
type 4 =
supermirror represented by trapezoidal reflectivity; 4 parameters
type 5 =
beryllium at 100K, including Bragg edges; no parameters
type 6 =
single

crystal filter, Freund formalism; 3 parameters
Complex regions:
type 10 =
multi

aperture collimator
type 11
=
multi

slit collimator, vertical blades; 3 or 5 parameters
type 12 =
multi

slit collimator, horizontal blades; 3 or 5 parameters
type 13 =
crystal monochromator; 10 parameters
*
type 14 =
toroidal mirror; 10 parameters
Time

dependent regions:
types 20.n
= chopper (disk or blade); 6 parameters
20.0 or 20.2 for motion in x

direction, 20.1 or 20.3 for y

direction
20.2 or 20.3 is counter

rotating (fully closed when edges at 0)
type 21 =
Fermi chopper (not yet implemented)
type 22 =
gravity focuser; 5 par
ameters
type 23 =
removable beamstop; no parameters
Scattering samples:
types 30.n = small

angle scattering sample; 2 parameters
30.0 scatters at constant Q; 30.1 is hard spheres of fixed radius
(formerly type 31)
type 32.n =
isotropic scatterer with
fixed energy change; 2 or 9 parameters
* 32.0 scatters into 4
, 32.1 has limited solid angle
type 34 =
inelastic scattering kernel; no parameters;
NAME
is ‘[path]filename’ of
S
(
) file
in MCNP Type I format
type 35 =
scattering from
layered reflectometry sample; 1 + 4
N
parameters
type 36 =
scattering from isotropic polycrystalline powder; 6 parameters + 2
table length
MCLIB
5
Detectors (zero

and one

dimensional):
type 40 =
detector; 9 parameters
type 41 =
vertical linear detector; 14 parameters
type 42 =
horizontal linear detector; 14 parameters
type 44 =
longitudinal linear detector; 14 parameters
Detectors (two

dimensional):
type 43.nm =
2

D detector in various coordinate systems; 22 parameters
43.00,
43.10, 43.20, rectilinear coordinates, respectively (X, Y), (Z, Y), (X, Z)
*
43.01, 43.11, 43.21, plane polar coordinates (
,
), respectively about axes
parallel to Z

, X

, and Y

axes
(Note: replace types 43.4, 43.5, and 43.6)
*
43.02, 43.12, 43.22, cyl
indrical coordinates, respectively (Z,
), (X,
), (Y,
)
*
43.03, 43.13, 43.23, spherical coordinates (
,
), with polar axis parallel to Z

, X

,
or Y

axes respectively
Scattering chamber:
type 50 =
scattering chamber, void

filled. No parameters, but
other regions may be embedded,
indicated by surface types with 10s digit on.
Sources:
types 90.n =
source size and phase space to be sampled
(Note: required!)
; 14

18 parameters
.1 = first aperture rectangular instead of circular, .2 = second aperture r
ectangular,
.4 = either/both aperture(s) offset vertically (Note: options are additive)
type 91 =
source energy distribution table and line shape parameters; 12 parameters plus
length of table
*
type 95 =
source file, direct

access binary; no parameters,
but NAME must be '[path]filename'
for the file. The file must be opened with logical unit number 95.
4. Program Structures
The relationships of the structures are shown schematically in fig. 2. The library source code is
available in Fortran 77 with VAX structure extensions (F77/VAX), and also in Fortran 90 (F90).
Unfortunately the F90 syntax did not adopt the relatively wid
espread F77/VAX
standard, most
distressingly by using “%” instead of “
.
” as the member pointer in a structure! In F77/VAX,
each surface is a
RECORD
of type
/SURFACE/
, and elements are referenced as (
e.g.)
SURFACE.G
.
In F90, the surfaces are defined as derived
TYPE
(SURFACE)
, and element references are of the
form
SURFACE%G
. Similarly, a geometric region is a
RECORD
of type
/REGION/
or a derived
TYPE
(REGION)
; the structure in either case is a vector
IGEOM
of 2

by
te integers, of length equal to
the maximum allowed number of surfaces. An additional structure,
MC_GEOM
, contains the
numbers of surfaces and regions in the problem,
NSURF
and
NREG
, and arrays of surface and
region structures. Information about the cont
ents of regions is contained in a structure called
MC_ELEMENT
,
which includes
NAME
and
INDEX
arrays, the parameter block
PARAM,
and the
pointer
NEXTINDX
to the next available location in
PARAM
.
The final structure is
PARTICLE
, which includes the position,
velocity, time of flight, mass (1 for
a neutron), charge (0 for a neutron), statistical weight, and vector polarization (not yet
implemented in the code). A purely “analog” Monte Carlo traces each individual neutron until
it is either lost or detected.
MCLIB uses “weighted” neutrons, and in many of the processes the
statistical weight is multiplied by the probability of survival instead of using a random number to
6
MCLIB
decide whether to terminate the history (“Russian Roulette”). This is especially beneficia
l when
scattering probability is small, as in subcritical reflection. To track
more long

wavelength
neutrons (which in general have larger scattering probability), the source distribution used is
2
I
(
) instead of
I
(
) and the initial weight is proportio
nal to 1/
2
. The tallied results are then
the sum of detected neutron weights. The relative error in each bin, however, depends on the
number of histories recorded.
The F77/VAX structure definitions are found in the file
MC_GEOM.INC
in Appendix A. The
differences required for strict F90 can be seen in file
MCGEOM90.INC
, also in Appendix A. As
compiled, the dimensions of the arrays are
MAXS
= 152,
MAXR
= 150, and
MAXP
= 1620.
5. Subroutine Library
Complete descriptions, calling sequences, revision history, and externals of all of the library
subroutines are given in Appendix B; the most up

to

date listing of the subroutine abstracts and
of the source codes may be found at
ftp://azoth.lansce.lanl.gov
/mclib/
. The following listing
divides the subroutines into (somewhat arbitrary) categories, and includes the latest modification
date and an abbreviated description. Note that if the operating system does not include the
function
RAN(ISEED)
, then such a
REAL*4
function must be provided as part of the library.
Vector analysis and transport
ANGLI
03 Feb 94
angle of incidence
ANGTORUS
04 Aug 96
angle of incidence on torus
DIST
09 Apr 96
distance to a surface
DISTORUS
08 Dec 96
dis
tance to a toroidal shell
SURFACE
MC_ELEMENT
A, B, C, D, E, F, G, P, Q, R, beta
PARAM
NAME
INDEX
pointer 1
TYPE
REGION
1
Problem (1st region) name
pointer 1
1
2
.
.
.
NSURF
2
Name of element in region 2
pointer 2
x
x
x
3
Name of element in region 3
0
void
x = { 0, ± [ {1,2,3,4,5,6} + {0,10} ] }
4
Name of element in region 4
pointer 3
pointer 2
TYPE
.
MC_GEOM
.
pointer 3
TYPE
NSURF
.
NREG
NREG
Last region name
pointer n
1
A, B, C, D, E, F, G, P, Q, R, beta
NEXTINDX
.
TYPE
2
A, B, C, D, E, F, G, P, Q, R, beta
.
.
NSURF
A, B, C, D, E, F, G, P, Q, R, beta
pointer n
TYPE
1
x
x
x
2
x
x
x
.
NEXTINDX
.
NREG
x
x
x
PARTICLE
X, Y, Z, VX, VY, VZ, TOF, M, Q, WT, PX, PY, PZ
Figure
2.
Structure
s used in MCLIB.
MCLIB
7
DTOEX
06 Mar 95
distance to nearest boundary
ELSCAT
26 Aug 95
do elastic scattering
LIGHTRFL
19 Mar 85
light reflection probability
LMONOCRM
04 Dec 97
reflection by crystal monochromator
LREFLCT
08 Sep 97
neutron reflection
probability
MOSAIC
27 Jan 95
angle of incidence with mosaic spread
MOVEX
02 Apr 96
move a particle
NEXTRG
29 Jun 95
find region across boundary
N_SOURCE
16 Dec 97
get source neutron
GET_SPACE
get phase space and brightness of source
OPERATE
04 De
c 97
find what happens to particle within region
EXIT_REG
find what happens to particle leaving region
OPTICS
03 Feb 94
photon at region boundary
RFLN
03 Feb 94
do reflection
SNELL
08 Sep 97
apply Snell's law
TESTIN
04 Apr 96
find if within region
WHICHR
03 Feb 94
find what region particle is in
WOBBLE
03 Feb 94
angle of incidence at wavy surface
Material properties and scattering functions
ATTEN_Al
28 Jan 93
attenuation of Al
ATTEN_Be
25 Jan 95
attenuation of Be
ATTEN_X
14 Aug 95
attenuation of single

crystal filter
KERNEL
03 Apr 95
inelastic scattering kernel
PLQSPHR
06 Mar 95
probability from spherical scatterer
POWDER
07 Dec 95
scatter from polycrystalline powder
REFLAYER
04 Mar 95
reflection probability from multiple laye
rs
Random distributions
LORRAND
04 Dec 97
random unit vector in limited solid angle
ORRAND
17 Nov 97
random unit vector in 4
PLCNVL
17 Nov 97
probability from convolution of top

hats
PLEXP
16 Feb 95
probability from exponential
PLIKCARP
17 Nov 97
pr
obability from Ikeda

Carpenter
PLMXWLN
25 Mar 96
probability from Maxwellian
PLNORM
17 Nov 97
probability from normal
PLNRMTBL
13 Feb 95
probability from table vs. normal distribution
PLORENTZ
17 Nov 97
probability from Lorentz (Cauchy)
PLPOISSN
17 No
v 97
probability from Poisson
PLTIME
05 Sep 95
probability from parameterized emission time
PLTRNGL
04 Jan 85
probability from triangle
RAN
various
random number; platform dependent
RAN0
19 Nov 97
desequentialized random number
RNDCRCL
17 Nov 97
random point in unit circle
SRC_PROB
20 Jun 97
probability density /meV/
s from type 91 table
General utilities
DET_2D
07 Jul 97
encode 2

D detector
DIGITS
05 Dec 87
convert number to string
GAMMLN
92
ln(
(x))
GET_BINS
18 Sep 97
find detector/time bin
s
GET_RHO
16 Oct 95
find scattering

length density
GINI
05 Feb 92
standard deviation by Gini statistic
HUNT
13 Apr 87
find index in array
INTERP
84
interpolate in array
LONGOUT
29 Jul 95
output array of
REAL*4
numbers in single record
8
MCLIB
NORM
78
norma
lize unit vector
READ_1D
18 Sep 97
read 1

D block ASCII data file
READ_2D
18 Sep 97
read 2

D block ASCII data file
REALOUT
02 Jan 96
output array of
REAL*4
numbers in block ASCII
Time

dependent devices
GRAV_FOC
07 Apr 94
gravity focuser
XCHOPPER
15 Ma
r 93
disk or blade chopper
Subroutine
OPERATE
and its other entry point
EXIT_REG
play special roles in the transport
process. As ca
n be inferred from the abstract in Appendix B, this subroutine is called when a
neutron enters a new region to determine what happens. Possible results on exit from
OPERATE
include detection or loss of the neutron, exiting the region with reduced statist
ical weight (partial
absorption), splitting the neutron into a transmitted and a scattered particle with the sum of
statistical weights equal to the original weight, or a transform of the coordinate system of the
problem. Coordinate transforms are applied
for region types made up of sub

regions (
e.g.,
Soller slits and benders), elements that change the beam direction (benders and
monochromators), and time

dependent devices (gravity focuser). In order to assure that the
coordinate system is properly restor
ed when the neutron leaves the region,
EXIT_REG
must be
called; this is flagged by using
5 as the integer defining all exit surfaces from such a region.
6. Neutron Source Functions
Subroutines
N_SOURCE
and
PLTIME
are called to generate source neutrons, using table lookup to
select the velocity and a parameterized model to select the emission time. (If the len
gth of the
lookup table is 1, the velocity is a triangular or delta

function distribution; if the thermal decay
constant is
0, emission time is 0.) The table and parameters are generally read from a file
derived either from experimental measurements or
from a detailed Monte Carlo simulation of a
target/moderator/reflector system. Several source files are available at
ftp://azoth.lansce.lanl.
gov/pub/mclib/tables
. The example given below uses a coupled liquid H
2
moderator with a
40

cm thick Be reflector
on a spallation source, computed with MCNP [8], but the procedure
E
(meV)
0.01
0.1
1
10
100
1000
10000
E
×
I
(
E
) (10
8
n / cm
2
/ ster / MW / s)
0.1
1
10
100
1000
10000
Figure
3.
Neutron source energy distribution, as computed with the MCNP code.
Coupled liquid
hydrogen moderator (equal fractions ortho and para), flux

trap geometry, 40

cm Be inner reflector, and Pb
outer reflector (both reflectors are D
2
O cooled). A po
wer law was fitted to the low

energy regime to
interpolate the histogram.
MCLIB
9
would be the same for an experimental spectrum measurement. The energy spectrum from the
MCNP output is shown in fig. 3. Program
MKNRMTBL
is provided to convert the spectrum. To
improve s
ampling at long wavelengths, the data are multiplied by
2
before the cumulative sum
is formed (the neutron is given an initial statistical weight proportional to 1/
2
to account for
this). The cumulative sum of the weighted table is compared to the integ
ral of a Gaussian (
vs
.
log
E
) fitted to the weighted data to form the table; that is, for equal steps in the argument of the
Gaussian, its integral is evaluated and the cumulative sum is interpolated to find the
corresponding neutron energy, which is then
tabulated as log
10
(
E
/ 1
meV). To sample the
neutron energy, a random deviate is chosen from a Gaussian distribution and that value is used to
interpolate in the table (“transform” method [9]). This provides sharp detail over the peak of the
spectrum and
coarser interpolation in the wings, which may extend to more than 5 standard
deviations.
The model for the time distribution is the sum of an epithermal and a thermal term convoluted
with a
Gaussian, added using a switch function [10,11] of the form exp[
(
0
/
)
P
]. This form
was chosen because the convolution is easy in Monte Carlo; note however that a function to
sample an Ikeda

Carpenter pulse shape distribution [10] has now been added to
the library.
Adequate fits can be made to most source terms using two exponentials; the epithermal (and the
two parameters of the Gaussian) are proportional to wavelength, and the thermal time constant is
independent of wavelength. Including the switch f
unction, the total number of parameters in the
model is then six. For the hybrid Be/Pb reflector case shown here, however, the fit with just two
exponentials is not adequate and the thermal term itself must be a linear combination of two
exponentials, for
a total of eight parameters. Neutrons escaping from the moderator surface in
the MCNP run [8] were tallied
vs.
time for 18 (logarithmic) wavelength bins from 0.09 to 15
Å;
H
2
(50/50,64/16), FluxTrap, 40cm Be, Pb
= 1.0 Å
t
(
s)
0
1000
2000
3000
4000
5000
I(t)
t
(
s)
0
10
20
30
I
(
t
)
s
1
)
10
3
10
2
10
1
10
0
10
1
10
2
2
= 170
s (67%)
1
= 610
s (33%)
0
= 0.86 Å
P
= 2.63
= 8.67
s/Å
= 5.68
s/Å
= 2.08
s/Å
Figure 4.
Source time distribution at 1.0 Å.
Points are the MCNP “data,” and the line is the
model fit.
At this wavelength, near the midpoint of the switch function, the early “epithermal” and the two late
“thermal” exponential terms are all seen.
10
MCLIB
fig. 4 shows the data for one such bin, from 0.9 to 1.2
Å.
Since the exponentia
l decay after
1000
s is constant for all wavelengths, the data from 0.9 to 15
Å
were summed and a
least

squares fit was made over the time range 1200
–
6000
s, giving a thermal time constant
1
= 610
s. Since the early (epithermal) parameters are all proportional to
, the early

time data
can be summed after binning
vs.
(
t
/
). (The variable
t
/
is equivalent to the reduced time
vt
used in ref. [10].) A least

squares fit in the
t
/
range 1
–
30
s/Å ga
ve the values epithermal
exponential
= 8.67
s/Å, delay
= 5.68
s/Å, and width
= 2.08
s/Å. Given these fits to
the early and late times, the additional exponential term was estimated (
= 170
s, ratio
R
=
2/3). Then the ratios of thermal to epi
thermal terms were extracted, and the two parameters of
the switch function determined (
0
= 0.86
Å,
P
= 2.63). The line in fig. 4 is an example of the
resulting fit, evaluated at
= 1
.
0
Å where all three components are significant.
The phase space sam
pled is determined by two apertures. Ordinarily the initial trajectory is
chosen by choosing a random point in each aperture. If the radius (or width) of the second
aperture is negative, however, then the second aperture defines a solid angle around the
Z

direction as seen from the first aperture, rather than a position.
7. Sample Program Output
Sample program source codes
LQDGEOM.FOR
and
MC_RUN.FOR
are available wherever the library
is archived, but because of their lengths are not included in this report. Communication between
the two programs is through a geometry file, such as the following example,
LPSS.GEO
. The
first part of the file is not u
sed by the succeeding program, but contains text to identify the
problem and to give the parameters entered by the user.
Geometry code version: Win95, July 9, 1997; P.A.Seeger
LQD Geometry file: lpss.geo
Proton pulse rate = 60.0 Hz, width = 1000.0 us,
max wavelength 15.1 A
Sample at 7.00 m, Detector at 12.95 m from moderator
Sample diameter = 10.00 mm
Spectrum and lineshape file:
\
mclib
\
tables
\
Pb40Be13.tbl
T0 Chop at 2.60m
, open
0.827ms,
close
14.839ms, rms 40.0us, v 56.55m/s
O'lap Chop at 3.10m, open 7.257ms, close 12.224ms, rms 20.0us, v 94.25m/s
Frame Chop at 6.00m, open 16.068ms, close 22.874ms, rms 20.0us, v 94.25m/s
Collimation: entrance aperture radius 7.53 mm at position 2.500 m
exit
4.44 mm 6.800 m
Moderator half

width and half

height = 0.050 m
spot radius 9.32 mm, penumbra radius 14.48 mm
Time range 34370.0

49665.0 us
Sample Transmission at Lmax (15.11 A) is 0.600
Total thickness of Al =
12.0 mm
Detector half

width 0.320 m, pixel 5.00 mm, rms 3.40 mm
penumbra radius 21.70 mm, min angle 4.78 mr
Q

range: 0.00200

0.03134 /A
Gravitational droop 2.21

4.62 mm
The next three lines give the numbers of surfaces and regions
and the length of the parameter
block. They are required, and the remainder of the file must follow in order:
SURFACES 29
REGIONS 28
PARAMETERS 203
MCLIB
11
Each surface definition is given on one line; zero entries may be left null:
,,,,,1/
PZ (plane perpendicular to z

axis) at z = 0
1,,1,,,,

.0036/
CZ (cylinder on z

axis), radius 60 mm
,,,,,1,

2.495/
PZ at z = 2.495 m
486749,,486749,,

1,15.5,

60.0625/
Cone on z

axis
,,,,,1,

2.5/
PZ at z = 2.5 m
,,,,,1,

2.6/
PZ at z = 2.6 m
1,,1,0.6,
,,

.2025/
Cylinder parallel to z

axis, offset y =
0.3 m, radius 0.335 m
,,,,,1,

2.9/
PZ at z = 2.9 m
1,,1,,,,

.0036/
CZ, radius 60 mm
,,,,,1,

3.1/
PZ at z = 3.1 m
1,,1,0.6,,,

.2025/
Cylinder parallel to z

axis, offset y =
0.3 m, radius 0.335 m
,,,,,
1,

3.102/
PZ at z = 3.102 m
,,,,,1,

6/
PZ at z = 6.0 m
1,,1,0.6,,,

.2025/
Cylinder parallel to z

axis, offset y =
0.3 m, radius 0.335 m
,,,,,1,

6.002/
PZ at z = 6.002 m
,,,,,1,

6.003/
PZ at z = 6.003 m
,,,,,1,

6.8/
PZ at z = 6.8 m
585518,,585518,,

1,6.8,

11.56/
Cone on z

axis
,,,,,1,

6.805/
PZ at z = 6.805 m
,,,,,1,

7/
PZ at z = 7.0 m
,,,,,1,

7.001/
PZ at z = 7.001 m
,,,,,1,

7.013/
PZ at z = 7.013 m
1,,1,,,,

0.25/
CZ, radius 0.5 m
,,,,,1,

13/
PZ at z = 13.00 m
,,,,,1,

12.95/
PZ at z = 12.95 m
1,,1,

.00638869,,,

5.51486

4/
Cylinder parallel to z

axis, offset y =
3.2 mm, radius 27 mm
,,,,,1,

13/
PZ at z = 13.00 m
1,,1,,,,

.1024/
CZ, radius 0.32 m
,,,,,1,

13.1/
PZ at z = 13.10 m
Each region is a line with
an entry for each surface:

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1

4

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
1 4

4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1

4

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 4

1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1

4

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 4

4 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1

4

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 4

4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
1

4

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 4

4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0

4 0 0 1

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 4 0 0 1

4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1

4

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 4

4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0

4 0 0 0 0 0 0 1

1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 1

4 0 0 0 0 0 0 0 0 0
0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

4

1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 4

1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

1 0 0 0 0
0 0 0 0 0
12
MCLIB
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1

4

1

11 11 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

1

1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

1

1
The
NAME
associated with each region; the first
NAME
is
the problem title, and the second is the
title of the source file:
7m, 13m, 15A, 60Hz (1000us), Q=0.01/A
H2 (50/50,64/16), Fluxtrap, 40

cm Be, Pb
Shutter penetration
Bulk shield
Circular entrance apertu
re
Collimator entrance
Drift to T0 chopper
T0 Chopper
T0 chopper shielding
Collimator pipe
Collimator shielding
Frame

overlap chopper
Overlap chopper shielding
Collimator pipe
Collimator shielding
Frame

definition chopper
Frame

definition chopper shield
Upstream Monitor
Collimator pipe
Collimator shielding
Circular exit aperture
Collimator exit
Drift to Sample
Sample
Aluminum
Secondary Flight Path
Beamstop
Detector
Next is the list of
INDEX
pointers into the parameter block:
1 20 0 135 0 136 0 137 144 0 145 146 153 0 154 155 162 16
3 0
173 0 174 0 175 178 179 180 181
Finally, the parameter block, including the source wavelength distribution:
90,

.05,.05,

.05,.05,2.5,.007525,6.8,.00444333,,8.9544E

5,.0818145,16666.67,

1000,20,,,,,
91,102,47.6,11.24,5.83333E16,610,170,.667,8
.67,5.68,2.08,.86,2.63,

1.700001,

1.698674,

1.69677,

1.694059,

1.69023,

1.684864,

1.677405,

1.667116,

1.653036,

1.633922,

1.608176,

1.579874,

1.544875,

1.49911,

1.453015,

1.394889,

1.336073,

1.27017,

1.198053,

1.128077,

1.054727,

.978697,

.90117,

.828609,

.75508,

.680891,

.608064,

.540398,

.472038,

.403717,

.342582,

.2799292,

.2202448,

.1641208,

.1068591,

.0565061,

.00302955,.0482707,.0975991,.1475286,.1969256,.2461485,.2951157,.341733,.387269,
.432616,.475812,.519135, .56161,.603647,.6452
07,.685663,.727038,.767961,.808221,
.852178,.893306,.937217,.978778,1.020489,1.062588,1.100546,1.141591,1.177922,1.213727,
1.252556,1.286119,1.326629,1.370402,1.412875,1.467081,1.520122,1.581781,1.650554,
1.722381,1.800968,1.889334,1.981707,2.078442,2.177055,2.278023,2.383607,2.487587,
2.585048,2.668664,2.739018,2.795598,2.847828,2.88669,2.918281,2.943239,2.960935,
2.973372,2.982054,2.988064,2.992237,2.995075,2.996995,2.998247,2.999082,2.999666,3,,,
20,5.65487E

5,707.355,827.355,14839.31,40,16666.67,,,
20,9.42478E

5,117.0279,7257.36,12223.82,20,16666.67,,,
MCLIB
13
20,9.42478E

5,82.747,16067.71,22874.19,20,16666.67,,,,
40,1,1.

4,15871.01,22933.77,37,.01,.1,,16666.67,,,
30,.00223634,.01, 2, 50, 23,
43.01,2
,1.89712,34370,49665,37,.01,.1,,16666.67,,.32,64,.005,.0034,

3.141593,3.141593,
32,.1963495,.0034,,

.00319434,13/
8.
Running the Problem: Program
MC_RUN
As distributed, the program
MC_RUN
has limits on some array sizes. The total space
for 2D (and
3D) histograms is 2,097,152 entries (12 memory bytes per entry), and for 1D histograms 65536
entries (20 bytes each). There are also limits of 16384 bins for any single axis, and a total of 16
detectors. All limits may be adjusted in a
PARAM
ETER
statement at the beginning of the program.
MC_RUN
asks the user for the geometry file name, the number of histories to be detected, and a
starting value for the random

number generator. The user chooses whether the recorded
monitor spectrum is to re
present a sample

in or a sample

out measurement. If the user gives a
surface number (from the geometry definition file) be monitored, an output
.MON
file will be
written with every neutron history that crossed the surface. If the surface number is given a
s
negative, the source neutron corresponding to each neutron crossing the surface will be recorded;
otherwise an additional question is asked for generating a correlation (
.COR
) file with either the
preceding neutron event or a preceding surface crossing.
The
.MON
file may be used as input in
subsequent runs by specifying the type of the source region as 95. A
.BAD
file may also be
generated which will contain complete histories of any neutrons which are flagged as “bad” by
some program module and which s
ubsequently hit the sample or a detector.
If there is enough memory space available, the user may choose to have 3D histograms, for
instance detector X, Y, and t. For very long runs, intermediate
.DAT
files are dumped
periodically to protect against comp
uter crashes. Options are also given to convert the final
output histograms, which are sums of weights, to integers. There are two ways to do this. If
one wishes the error bars on the output to be similar to a Poisson distribution so they look like
real
data, and if the error from the number of histories is small, then the value for each bin can be
replaced by a sample from a Poisson distribution with mean equal to the tally for that bin. If the
error bars are already appropriate, then a simple Russian
roulette procedure is to add a random
number (on the range 0
1) and truncate the result to an integer.
Here is sample output of
MC_RUN
for the geometry file created in the preceding section. The
running time to detect 300000 neutron histories (4.96×10
6
s
tarts) was 14 minutes on a
PentiumPro/200

Mhz PC (a factor of 12 slower on a 486/66); the estimated experiment time is
9.9 seconds (at 1 MW proton beam power).
MONTE CARLO SIMULATION OF A NEUTRON SCATTERING INSTRUMENT
MC_RUN code version: Win95, July 9,
1997, P.A.Seeger
(At this point the text from the front of the geometry file is reproduced)
Successfully read input file with 29 surfaces, 28 regions,
and 203 element parameters.
Region: Upstream Monitor
detector # 1, bins = 1 x 1 ( 0 radial) and 37 times
output data forms: t (spectrum), t (1D)
Region: Detector
detector # 2
, bins = 64 x 32 ( 0 radial) and 37 times
14
MCLIB
output data forms: t (spectrum), t (1D), RT (2D)
A test particle is sent along the instrument axis as a simple sanity check:
Beam elements along the instrument axis:
Surf 1, Z =
.000 m, enter reg 3, Shutter penetration
Surf 3, Z = 2.498 m, enter reg 5, Circular entrance aperture
Surf 5, Z = 2.500 m, enter reg 7, Drift to T0 chopper
Surf 6, Z = 2.600 m, enter re
g 8, T0 Chopper
Surf 8, Z = 2.900 m, enter reg 10, Collimator pipe
Surf 10, Z = 3.100 m, enter reg 12, Frame

overlap chopper
Surf 12, Z = 3.102 m, enter reg 14, Collimator pi
pe
Surf 13, Z = 6.000 m, enter reg 16, Frame

definition chopper
Surf 15, Z = 6.002 m, enter reg 18, Upstream Monitor
Surf 16, Z = 6.003 m, enter reg 19, Collimator pipe
Surf 17, Z = 6.800 m, enter reg 21, Circular exit aperture
Surf 19, Z = 6.805 m, enter reg 23, Drift to Sample
Surf 20, Z = 7.000 m, enter reg 24, Sample
Surf 21,
Z = 7.001 m, enter reg 25, Aluminum
Surf 22, Z = 7.013 m, enter reg 26, Secondary Flight Path
Surf 25, Z = 12.950 m, enter reg 27, Beamstop
Surf 27, Z = 13.000 m, enter reg 28, Detector
** No region found across surface 29, Z = 13.100
Having found the source, the sample, and the detector, it is safe to proceed:
Initial random number = 075BCD15, 300000 neutr
ons to be detected
Moderator phase space to be sampled = .00060 mm**2

ster
x 6.82 e

fold of energy
Source brightness = 5.833E+16
n/ster/m**2/MW/s
t 7m, 13m, 15A, 60Hz (1000us), Q=0.01/A MC__970709_155107
S
ource# Detected X(m) Y(m) Z(m) t(ms) Reg Surf
200 .000

.006 3.100 2.564 12 10
600 .006 .000 3.100 2.717 12 10
1400

.002 .001 6.000 13.8
47 16 13
2085 125

.246 .114 13.000 46.837 28 24
3685

.001

.002 3.100 4.858 12 10
4126 250 .078

.097 13.000 43.527 28 24
(This output, also displayed in r
eal time on the screen, continues with decreasing frequency)
Summary of neutron losses by region:
14.6739% lost in region 3, Shutter penetration
.0450% lost in region 6, Collimator entrance
12.8217% l
ost in region 8, T0 Chopper
71.2384% lost in region 12, Frame

overlap chopper
.6797% lost in region 16, Frame

definition chopper
.0000% lost in region 18, Upstream Monitor
.0003% lost in region 22, Collimator exit
.0000% lost in region 24, Sample
.0675% lost in region 25, Aluminum
.0017% lost in region 27, Beamstop
Total histories started: 4960718 (3.448E+08 neutrons)
Detector 1: Upstream Monitor
beam centered at (

.01, .14) mm, with rms ( 1.9, 1.9) mm
det
ected in Transmission: 312680 (2.296E+03 n, .00%)
MCLIB
15
detected in Scatter Mode: 0 (0.000E+00 n, .00%)
Detector 2: Detector
beam centered at (

.04,

3.12) mm, with rms ( 7.8, 7.8) mm
detected in Transmis
sion: 308449 (1.147E+06 n, .33%)
detected in Scatter Mode: 300000 (4.504E+05 n, .13%)
TOTAL in Scatter Mode: 300000 (4.504E+05 n, .13%)
Beam power on target: 9.90 Mw

s
Bad

Frame neutrons passed: 0.000E+00
Bad

Fr
ame neutrons detected: 0.000E+00
Final random number: 5FBFF52E
The output data files include header information and blocks of ASCII data for bin boundaries,
data, and standard deviations. The values in the blocks are “,” separated and “/” termin
ated as
in the geometry file shown above. Data include the transmitted spectrum histogram and a
2

dimensional histogram of counts integrated in radial zones on the detector for each time slice
(fig. 5). By converting the radial data to
Q
for each individual wavelength (time slice) before
combining, the time

of

flight wavelength resolution is maintained.
9.
Analysis Routines
In general, the user will be expected to use his own analysis rout
ines, the same as if the data
originated from a real instrument. This will generally require a conversion of the data formats.
Sample programs to convert
.DAT
output files are available from the same source as the
MCLIB
library,
ftp://azoth.lansce.lanl.g
ov/pub/mclib/analysis/
. Specifically, the program
SWITCH
will
read a 1D file and convert it to column format, and the program
SPREAD
will convert 2D data to a
2
(°)
0
1
2
3
4
(Å)
11
12
13
14
15
0.1
0.1
1
1
10
10
10
10
1
1
1
0.1
0.1
0.1
0.1
Figure 5.
Relative count rate in detector
vs.
radial position and time of flight.
Time of flight has
been converted to wavelength and detector radial zone to scattering angle (2
). Contours are logarithmic.
16
MCLIB
spreadsheet format. For 3D files, the sample program
QLfromXY
is given; this is a specific
appl
ication for small

angle scattering, based on the data reduction program
SMR
that is used for
the LQD instrument at Los Alamos. Any of these may be modified to fit a user’s needs. For
visualization of the data, programs
QUIKPLOT
and
QUIK2D
are provided fo
r 1D and 2D files
respectively, using the
PGPLOT
graphics library, free from
ftp://astro.caltech.edu/pub/pgplot/
.
While the format of the
.DAT
file represents a real measurement, a
.MON
file includes details that
could not be measured, such as correlatio
ns between actual wavelength and the wavelength
measured by time of flight. To extract this information from a
.MON
file (direct

access binary),
use program
SUPER_KNOW
with modifications for your specific needs. Comments in the source
code indicate where
the user

specific code is to be inserted.
There are also additional subroutines in this folder that are not included in
MCLIB
itself.
NEWHEAD
21 Jan 88
concatenate strings before and after operation string
NEWOP
21 Jan 88
attach operation character t
o operation string
OPLIST
28 Jan 88
concatenate 2 operation strings with operator
PAINT
20 Jun 97
output a 2D plot with
PGPLOT
drivers
PLOT10
06 Oct 92
convert most (38)
PLOT10
primitive calls to
PGPLOT
CHECK_TEXT
decode
PLOT10
(
DISSPLAY
)
special characters to
PGPLOT
CONTOUR
contour plot of an array
CURVE
plot with symbols and various dashed lines
ERRBARS
plot linear or log error bars
LINAXIS
scale ends of linear axis to rounded values
LOGAXIS
scale ends of lo
garithmic axis to rounded values
PREMUL
12 Feb 91
attach multiplier to front of operation string
QBINS
04 Sep 89
generate Q

bin boundaries on log/lin/log scale
RBAR
15 Jun 94
average radius of a rectangular pixel
SPLINE
24 Jun 96
fit cubic spline to points
SPLINT
24 Nov 84
interpolate using cubic spline
WRTABSTR
02 Feb 92
write an LQD standard file abstract
10. Creating a New Type
The process of adding toroidal mirrors will be used as a case study. There were four steps:
1.
Define a new region type (14) and assign parameters and names in
MC_ELMNT.INC
.
2.
Develop needed algorithms and procedures (
DISTORUS
and
ANGTORUS
).
3.
Insert the “methods” for the new type in subroutine
OPERATE
.
4.
Test, debug, and refine the previous t
hree steps.
1.
There is renewed interest in the use of focusing mirrors in neutron instruments, for
example the proposal of Benno Schoenborn for a Laue protein

crystallography instrument [12],
and a recent paper by John Copley [13]. Although an ellipsoid
would be the optically preferred
shape, a toroid may be more practical to manufacture. Since a toroid can not be defined as a
surface
by eq. (1), it was necessary to define it as a
region
. The nature of the region is that it is
divided into two subregion
s by the toroidal surface. The ten parameters can be seen in
Appendix A at type
14. The torus is defined by its major radius
R
, the radius
r
of the generating
circle, the offset of the center of the generating circle from the beam axis, and the longitudinal
position of the center of revolution. Also, the sign of the major radius shows which side of the
beam axis the center is on, and the sign
of the minor radius is used as a flag for orientation: + for
horizontal and
for vertical. (The choice of parameter three as a small offset instead of
distance from the center of the torus was crucial to maintaining precisio
n.) The fifth parameter
is surface roughness, and then three parameters are provided to define a rotation of the
MCLIB
17
instrument beam axis after reflection. Finally, there are two pointers to define the materials
inside and outside the torus. The “inside” ma
terial is in front of the mirror surface, and is
usually void, indicated by a pointer of zero. The “outside” is the mirror surface/substrate, and
the pointer is the offset of its definition in the
PARAM
block relative to the torus parameters.
E.g.
, if th
e material type immediately follows this block, then the pointer value is 12 (or 1 +
NUMBER_14
, the number of cells in
PARAM
used by type 14).
2.
Finding the intersections of a line with a torus (a quartic surface) is not a trivial problem.
Previous codes
[13,14] have solved the problem very generally with sophisticated methods to
retain precision. My approach was to begin with three limiting assumptions and one
complication:
only the near half of the toroid is relevant
only the outer (concave) half

shell of the toroid is relevant
only real solutions for the intersections are relevant
the particle trajectory is a gravitational parabola instead of a straight line.
At this stage I had the assistance of Mike Fitzsimmons and Tom Klugel (MLNSC) to s
tudy
methods of setting up and solving the fourth

order equation. The first method we used to set up
the equation was copied from [14], using analytic derivatives of the trajectory in the coordinate
system of the torus. This ran into serious difficulties
when gravity was included (see §4 below),
and the eventual solution was first to bracket the nearest and furthest intersections (if any!) of the
parabola with bounding cylinders and planes of the relevant portion of the torus, and to solve for
the distanc
e from the surface parametrically
vs.
time of flight at five equally spaced points.
Since the spacing of the points is uniform, multiplication by a predetermined constant matrix
solves for the coefficients of the fourth

order equation of nearest distance
as a function of time.
One or two real roots (if any) are then found by the “safe” Newton

Raphson method discussed in
Numerical Recipes
[15], which uses bisection to keep the roots bracketed whenever the iteration
step would move outside the current brack
et. For further details, see the source code for
DISTORUS
in file
MCLIB.SRC
at
ftp://azoth.lansce.lanl.gov/pub/mclib/
.
Finding the angle of incidence, subroutine
ANGTORUS
, is far simpler. It is essentially a copy of
WOBBLE
with derivatives of the torus equation instead of eq. (1). As in
DISTORUS
, precision is
enhanced by accounting for the fact that the major radius of the torus is large compared to any
other distance.
3.
Inclusion of type

specific code in subroutine
OPER
ATE
requires adherence to the overall
structure of the routine, and will probably be the most difficult (or dangerous) aspect for a user to
add a function to MCLIB. Here is the scenario:
When a particle enters a region,
OPERATE
is called with the follo
wing references:
PART
= record containing description of particle (input/output)
EXDIST
= distance to exit surface particle is aimed at (m) (input/output)
PARAMS
= array with description of what is in this region (input)
GEOM
= structure with all surface a
nd region geometry definitions (input)
IREG
= region number of device, or subregion within device (input/output)
JSURF
= surface number, if particle is initially on surface (input)
KSURF
= surface number that particle is pointed toward (input/output)
NAME
= name of region, may used as file name (
c.f.
type 34) (input)
TRANSMIT
= flag to compute transmission of sample types 30

39 (input)
FLAG
= flag set to .FALSE. if (
e.g.
) chopper in wrong frame (output)
PART_2
= description of particle created by operation (output)
18
MCLIB
DET_WT
= statistical weight of detected particle (output)
IX, IY
= position bin numbers of detected particle (output)
ISEED
= random

number generator seed (input/output)
OPERATE
determines the reg
ion type from
PARAMS
(1) and transfers to the appropriate block
of code in an
IF

THEN

ELSE IF
structure (this will become a
CASE
structure when F90 is
more widespread).
The code may use any variables marked “input” and may change any variables marked
“ou
tput” in the above list. Local variables will have to be declared, and any variables to be
retained across entries must have distinctive names and be placed in a
SAVE
statement.
Parameter names from
MC_ELMNT.INC
and physical constants from
CONSTANT.INC
(see
Appendix A) are available for use.
Valid operations include: moving the particle to the exit from the region with possible
reeducation of statistical weight to account for absorption (set
EXDIST
to 0); selection of a
subregion without moving; refl
ection (set
KSURF
negative); rotate the coordinate system;
statistically split the particle or create a new one (in
PART_2
); scatter either elastically or
inelastically; multiple scattering, keeping control through several motions of the particle till it
r
eaches the region exit; detect the particle and find encoded position in detector. Examples
of all these possibilities can be found in
OPERATE.FOR
.
Algorithms may be deterministic or probabilistic (or both). Various random deviate
functions are availa
ble (and others may be added to the library). For single random
numbers, use
RAN(ISEED)
, or
RAN0(ISEED)
to eliminate correlations.
For compatibility with older Fortran compilers without recursive subroutines, the code
should not make any calls to
OPER
ATE
for processing subregions. Use no numbered
statements. Style should be similar to the existing
OPERATE.FOR
code.
In the case of the toroidal mirror, the actions are to determine the distance to the torus and which
side the particle is on using
DISTORUS
, and to move either to that surface or the region exit
(whichever is closer) through the appropriate material, compute reflection probability using
ANGTORUS
and the ratio of scattering length densities across the surface, and repeat all of these
a
ctions till the region is exited. Note that
only
material types 0 (void), 1 (amorphous
non

magnetic), or 4 (supermirror layer) are allowed, due to the restriction that all code for the
type must be inline. In fact, the code for type 1 and type 4 regions
has been copied into this
block.
4.
The debugging of this module is an example of cooperation with outside users. The
code had been tested by reproducing the results in [13], but in a different coordinate system.
When John Copley tried to use the code,
there were subtle errors. It took one more test to show
that my code was at fault, rather than John’s interpretation of it, two more rounds of modification
and testing to uncover a serious philosophical flaw in the algorithm, and two more rounds to
“perfec
t” the code. This version passes all tests, and is twenty times faster [16] than the
program used previously.
11. A Focusing Mirror in a Small

Angle Scattering Instrument
The original paper on curved

mirror neutron optics [17] suggested small

angle sca
ttering as an
application. The authors considered a bent mirror forming a cylinder with axis transverse to the
beam, and a narrow slit aperture with the solid angle of the mirror matched to the solid angle of
the neutron source. The same geometrical matc
hing can be applied to illuminate a mirror that
MCLIB
19
focuses in two dimensions, as illustrated in Fig. 6. (Note that the transverse dimensions in Fig.
6 are exaggerated by a factor of 10 relative to the longitudinal dimensions.) The MCLIB code
allows us to co
mpare mirror shapes in this geometry: the ideal ellipsoid, the tangent toroid,
cylindrical segments (with longitudinal cylinder axes), or toroidal segments (bent cylinders with
varying curvatures). Compared to a two

aperture collimation system with the sa
me resolution, a
factor of 10 increase in intensity is possible because the full illuminated surface of the moderator
may be viewed. The tradeoff between intensity and resolution depends on the single aperture.
The mirror considered as an example was 3.0 m long and 75 mm wide, with a major axis of 9 m
and t
ransverse (minor) axes of 157 mm. The glancing is 1º at the center (beam axis bent by 2º).
The radius

of

curvature parameters of the tangent toroid [13] are
R
= 257.77 m and
r
= 78.54
mm. The source aperture had a diameter of 2.0 mm, located at one focu
s of the ellipsoid; at this
Figure 7.
Images of a 2

mm diameter aperture for ellipsoidal (left side) and toroidal
(right side) focusing mirrors.
The r e f l e c t i on pl a ne i s ve r t i c a l ( up wa r d) t o c or r e c t f or
gr a vi t y. Ea c h pa t t e r n i s a hi s t o gr a m f or 250000 de t e c t
e d ne ut r ons wi t h wa ve l e ngt h
di s t r i but i on f r o m 2
20 Å, f r om a c oupl e d l i qui d

hydr oge n mode r a t or. Ea c h f i gur e i s 20 mm
s qua r e on t he de t e c t or ( pr i nt e d a t f our t i me s s i z e ), a nd t he i nt e ns i t y s c a l e i s l oga r i t hmi c f r om 1
t o 93250 c o unt s pe r pi xe l.
Fi gure 6.
Smal l

Angl e Scat t eri ng i ns t rument us i ng a f ocus i ng mi rror f or i ncreas ed
i nt ens i t y.
20
MCLIB
eccentricity, the focus is less than 1 mm from the vertex. Also, the transverse diameter of the
ellipsoid at the focus is only 2.7 mm, so the source is “large” in terms of optical properties, and
Monte Carlo is a useful tool a
estimating resolution at the detector placed at the other focus.
The detector pixel size was made small (0.25 mm) and the encoding uncertainty was zero so as
not to affect the computed resolution when compared to the source aperture size, which
contribut
es 0.50 mm to the standard deviation in each detector coordinate. Since the effect of
gravity is always included in the algorithms of MCLIB (for any particles with rest mass), it is
necessary to place the source aperture (and the detector center) below th
e instrument axis
if the
reflection plane is horizontal
by a distance which is proportional to wavelength, such that the
nominal trajectory at the mirror center is horizontal. However,
if the reflection is upwards
(mirror horizontal and concave up), then
no correction is necessary because the excess
downward velocity acquired before striking the mirror is reflected and then canceled by
downward acceleration after reflection (this cancellation is exact only for neutrons reflected at
the center of the mirror
, but is independent of wavelength). For a white source (
i.e.,
spallation)
the vertical reflection plane is much to be preferred.
A comparison of the beam spots for the ellipsoid and toroid is shown in Fig. 7, for vertical
reflection of a cold

moderator
spectrum from 2
20 Å. (A supermirror reflecting layer was
assumed, and reflectivity is good above 3 Å.) The standard deviations for the ellipsoid are 0.64
mm in both the horizontal and vertical directions, and those for the toroid are 1.05 mm horizontal
and 2.29 mm vertical. A minimum Q

value less than 0.001 Å
1
could be achieved with this
geometry in an instrument of total length 13 m (with the ellipsoidal mirror). Further details of
this simulation will be presented elsewhere.
12. Future Directions
It is hoped that the MCLIB code will be generally useful to designers of neutron instruments.
For this goal to be realized, two sets of improvements are being pursued. First, the general user
interface has been moved to a web page (
http://bayberry.lanl.
gov/lansce/Welcome.html
) so that
geometry files may be constructed more easily than is now possible. The MCLIB library will be
modified as necessary to support this project. Second, we are soliciting input of ideas,
algorithms,
and/or
(public) code modul
es from interested or potential users, so that the code will
do what
you
want. All code included will be acknowledged in the subroutine listings and
abstracts. Please communicate with the author by e

mail to
PASeeger@aol.com
if you have any
questions or
suggestions.
It is the intention of the author that MCLIB remain in the public domain (see copyright notice in
Appendix B). The most recent versions of this document and the library and auxiliary codes as
described above will continue to be available to
any interested users from the author, or from the
Manuel Lujan Jr. Neutron Scattering Center at
ftp://azoth.lansce.lanl.gov/pub/mclib
.
Acknowledgments
This effort has been supported historically by the U. S. Department of Energy under contract
W

7405

ENG

36, and the author is grateful to Rex Hjelm, Luke Daemen, and the Manuel Lujan
Jr. Neutron Scattering Center for their continued interest and financial support.
MCLIB
21
References
[1]
P. A. Seeger, “The MCLIB Library: Monte Carlo simulation of neutron scatter
ing
instruments,” 13th Meeting of the International Collaboration on Advanced Neutron
Sources, October 11
–
14, 1995, Paul Scherrer Institut, Villigen, Switzerland, PSI
Proceedings 95

02, pp. 194
–
212.
[2]
P. A. Seeger, “The MCLIB Library: new features,” Work
shop on Methods of Neutron
Scattering Instrument Design,” Lawrence Berkeley Laboratory, Berkeley, CA, Sept 23
–
27,
1996.
[3]
M. W. Johnson and C. Stephanou, “MCLIB: A Library of Monte Carlo Subroutines for
Neutron Scattering Problems,” Rutherford Laboratory
report RL

78

090, September, 1978;
M. W. Johnson, “MCGUIDE: A Thermal Neutron Guide Simulation Program,” Rutherford
and Appleton Laboratories report RL

80

065, December, 1980.
[4]
P. A. Seeger and R. Pynn, “Resolution of Pulsed

Source Small

Angle Neutron
Scattering,”
Nucl. Inst. Methods
A245
(1985) 115
124.
[5]
P. A. Seeger, “Monte Carlo Calculations for the Optimisation of the Beam Optics of the
LOQ Spectrometer,” report to Rutherford Appleton Laboratory, 28 April, 1994.
[6]
P. A. Seeger, G. S. Smith, M.
Fitzsimmons, G. A. Olah, L. L. Daemen, and R. P. Hjelm,
Jr., various contributions, Instrumentation for a Long

Pulse Spallation Source Workshop,
Lawrence Berkeley Laboratory, Berkeley, CA, April 18
21, 1995.
[7]
T. G. Thelliez, L. L. Daemen, P. A. Seeger,
and R. P. Hjelm, Jr., “A user

friendly geometry
interface for the Monte Carlo neutron optics code MCLIB,” 13th Meeting of the
International Collaboration on Advanced Neutron Sources, October 11
–
14, 1995, Paul
Scherrer Institut, Villigen, Switzerland, PSI P
roceedings 95

02, pp. 307
–
311.
[8]
E. J. Pitcher, G. J. Russell, P. A. Seeger, and P. D. Ferguson, “Performance of long

pulse
source reference target

moderator

reflector configurations,” 13th Meeting of the
International Collaboration on Advanced Neutron S
ources, October 11
–
14, 1995, Paul
Scherrer Institut, Villigen, Switzerland, PSI Proceedings 95

02, pp. 323
–
329.
[9]
W. H. Press, S. A. Teukolsky, W. T. Vettering, and B. P. Flannery,
Numerical Recipes,
2
nd
Edition (Cambridge University Press, 1992), sectio
n 7.3.
[10]
S. Ikeda and J. M. Carpenter, “Wide

energy

range, high

resolution measurements of
neutron pulse shapes of polyethylene moderators,”
Nucl. Inst. Methods
A239
(1985)
536
–
544.
[11]
R. A. Robinson and J. M. Carpenter, “
On the use of switch functions in describing pulsed
neutron moderators,”
Nucl. Inst. Methods
A307
(1991) 359
–
365.
[12]
B. P. Schoenborn and E. Pitcher, “Neutron diffractometers for structural biology at
spallation neutron sources,” in
Neutrons in Biology,
B. P. Schoenborn and R. B. Knott,
eds., Plenum Press, New York, pp. 433
444 (1996).
[13]
J. R. D. Copley, “Simulation of neutron focusing with curved mirrors,”
Rev. Sci. Instrum.
67
(1996) 188
–
194.
[14]
G. Hummer, subroutine
TORUS
, Feb. 28, 1995.
[15]
W. H. Press
et al., op. cit.
section 9.4.
[16]
C. Lartigue, private communication.
[17]
H. Maier

Leibnitz and T. Springer, “The use of neutron optical devices on beam

hole
experiments,”
J. Nucl. Energy Parts A/B
17
(1963) 217
–
225.
MCLIB
A

1
Appendix A
Files MC_GEOM.INC, MCGEOM90.INC, MC_ELMNT.INC, and CONSTANT.INC
MC_GEOM.INC
C 08 Feb 1994: converted from COMMON to STRUCTUREs. Added WT, M,
and POL to
C particle record. Require Y direction to be vertical (for
C gravity). [PAS]
C 17 Feb 1994: added BETA to surface record; type 5 boundary; made MAXP
C a parameter [PAS]
C 15 Feb 1995: include structure /MC_ELE
MENT/
formerly
in MC_ELMNT.INC
[PAS]
C 04 Dec 1995: larger MAXR, MAXP; replace POL with Q and (PX,PY,PZ) [PAS]
C 27 Nov 1997: larger MAXS, MAXR, MAXP (total size 65530 bytes) [LLD,PAS]
C
C To change dimensions of all arrays, change maximum numbers of
C
surfaces (MAXS), regions (MAXR), and/or element parameters (MAXP) in
C following PARAMETER statement.
INTEGER*4 MAXS, MAXR, MAXP
PARAMETER (MAXS=152, MAXR=150, MAXP=1620)
C
C Definitions of surface parameters and array of boundaries of s
urfaces:
C SURFACE is a record containing 10 coefficients of a general
C quadratic surface, of the form
C A*(X**2) + B*X + C*(Y**2) + D*Y + E*(Z**2) + F*Z + G +
C P*(X*Y) + Q*(Y*Z) + R*(Z*X) = 0
C BETA = surface roughness param
eter between 0 (smooth) and 1 (cosine);
C negative or >1 is completely random
STRUCTURE /SURFACE/
REAL*4 A, B, C, D, E, F, G, P, Q, R, BETA
END STRUCTURE
C
C REGION is record containing NSURF values of the I*2 variable I
GEOM,
C defining a region by its bounding surfaces:
C + if interior of region is on + side of surface
C

if interior of region is on

side of surface
C 0 if surface is not a boundary of the region
C 1 for ordinary surface described by roughness BETA and possibility
C of refraction or critical reflection
C 2 for totally reflecting surface
C 3 for diffuse scattering surface
C 4 for absorbing surface
C 5 special
action required in previous region before crossing surface
C 6 split particle by 2 after crossing surface, treated as type 1
C +10 for a surface which bounds an embedded region within the region; not
C tested to determine if within regio
n, but tested for exit from region
STRUCTURE /REGION/
INTEGER*2 IGEOM(MAXS)
END STRUCTURE
C
C NSURF = number of surfaces defined
C NREG = number of regions defined
STRUCTURE /MC_GEOM/
INTEGER*4 NSURF,NREG
RECORD /SURFACE/ SURFACE(MAXS)
RECORD /REGION/ REGION(MAXR)
END STRUCTURE
A

2
MCLIB
C
C Definitions of element

parameter block
C Each INDEX points to the beginning of a structure within the contiguous
C block PARAM. The first entry in PARA
M for each element identifies the
C type, followed by a varying number of parameters.
C NAME = 40

character descriptive name of a region or element
C NEXTINDX = pointer to next available location in PARAM
C INDEX = pointer into parameter block; if 0, region is not an "element"
C PARAM = block of element parameters
STRUCTURE /MC_ELEMENT/
CHARACTER*40 NAME(MAXR)
INTEGER*4 NEXTINDX, INDEX(MAXR)
REAL*4 PARAM(MAXP)
END STRUCTURE
C
C Definitions of position and velocity of particle:
C (X,Y,Z) = position of particle (m); note that +Y is UP (for gravity)
C (VX,VY,VZ) = velocity of particle (m/us)
C TOF = time of flight of particle (us)
C M = atomic number o
f particle (e.g., 1 for neutron, 0 for photon)
C Q = charge number of particle (e.g., 0 for neutron, +1 for proton)
C WT = statistical weight of particle
C PX,PY,PZ = average polarization vector, from

1.0 to +1.0
STRUCTURE /PARTICLE/
REAL*4 X, Y, Z, VX, VY, VZ, TOF
REAL*4 M, Q, WT, PX, PY, PZ
END STRUCTURE
C
MCGEOM90.INC
INTEGER (4):: MAXS, MAXR, MAXP
PARAMETER (MAXS=152, MAXR=150, MAXP=1620)
C
TYPE SURFACE
SEQUENCE
REAL (4):: A, B, C, D, E, F, G, P, Q, R, BETA
END TYPE SURFACE
C
TYPE REGION
SEQUENCE
INTEGER (2):: IGEOM(MAXS)
END TYPE REGION
C
TYPE MC_GEOM
SEQUENCE
INTEGER (4):: NSURF, NREG
TYPE (SURFACE) SURFACE(MAXS)
TYPE (REGION) REGION(MAXR)
END TYPE MC_GEOM
C
TYPE MC_ELEMENT
CHARACTER NAME(MAXR)*40
INTEGER (4):: NEXTINDX, INDEX(MAXR)
REAL (4):: PARAM(MAXP)
END TYPE MC_E
LEMENT
C
TYPE PARTICLE
MCLIB
A

3
SEQUENCE
REAL (4):: X, Y, Z, VX, VY, VZ, TOF
REAL (4):: M, Q, WT, PX, PY, PZ
END TYPE PARTICLE
C
MC_ELMNT.INC
C Definitions of beam elements which may occur in regions, and their parameters
C
C P. A. Seeger, April 20, 1994
C 04 Jan 1995: define offsets of parameters within blocks, rather than
C structure with UNIONs [PAS]
C 10 Jan 1995: modified source type 90; added type 91 for source spectrum
C and l
ineshape description [PAS]
Comments 0
Log in to post a comment