THE MCLIB LIBRARY: MONTE CARLO SIMULATION OF NEUTRON SCATTERING INSTRUMENTS

haddockhellskitchenUrban and Civil

Nov 15, 2013 (3 years and 11 months ago)

201 views



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, 40-cm 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]