Keywords

In this section, we will explain the basic MCNP® concepts necessary to prepare the input deck for the MCNP code to solve radiation transport problems. Here we assume that the readers have the MCNP software and a text editor available on their computers.

2.1 Geometry

2.1.1 Simplest Possible Input File

The MCNP code [1,2,3] runs by reading an ASCII input file that describes the problem geometry, materials, sources, tallies, physics, and options. Example 2.1 is the MCNP input file—the simplest possible input file.

Example 2.1 The Simplest Possible Input File

Trivial Example 1 0 -11 imp:n=1 2 0 +11 imp:n=0 11 SPH 0 0 0 5 SDEF

The MCNP input deck consists of card images. Each line is called a card for historical reasons.

The first card of an MCNP input deck is the title card. Here, the title is “Trivial Example.” The title line is followed by three sections of cards: cells, surfaces, and data. These sections are separated by a blank card. In Example 2.1,

1 0 -11 imp:n=1

specifies cell 1 with material 0 (vacuum) inside surface 11 with neutron importance 1 (see section 1.3.4.2 of the MCNP manual).

2 0 +11 imp:n=0

Cell 2 has material 0 and is outside surface 11 with importance 0, which means any particle that passes into cell 2 is no longer followed (“killed”). Thus cell 2 is the outside world. All MCNP geometries extend out infinitely, so cell 1 is everything inside surface 11 and cell 2 is everything outside.

After the blank line delimiter, the surface card(s) are defined.

11 SPH 0 0 0 5

describes surface 11 as a sphere at coordinates (0,0,0) with radius 5 cm (MCNP manual 3.2.2).

After the blank line delimiter, the data cards are defined. MCNP has hundreds of data cards that are defaulted but, at minimum, the source definition, SDEF, must be provided. The defaults are 14 MeV neutrons isotropic at the origin (0,0,0).

The default units for length, energy, time, and mass are cm, MeV, shakes (1 × 10−8 s), and grams, respectively (MCNP manual 7.1.6).

Comment cards begin with a c followed by a space and may be anywhere after the title card. Each card may also have a comment at the end started by “$.” The “$” and what follows are ignored. Capitalization is ignored except when upper or lower case letters appear in the name of a file. Anything entered after column 128 is ignored. A card beginning with five or more blank spaces is a continuation of the previous card. A card that ends with “&” is continued on the next line regardless of whether the five or more blank spaces are present on that line. Entries may appear in any column otherwise.

2.1.2 Running MCNP: The Simplest Case

The input file is run from a Windows (or from a shell in LINUX/UNIX systems) command prompt as:

MCNP6 i=EX11_1

where the name of the text input file in this case is EX11_1 (beware of hidden “txt” extensions in Windows systems). This problem will run until terminated by an interrupt. The possible interrupts are

<ctrl-c><s>

See status of calculation—how many histories are run, etc.

<ctrl-c><q>

Quit cleanly—finishing last history and writing output files

<ctrl-c><m>

Interrupt—to plot geometry or tallies

<ctrl-c><k>

Kill run—do not use unless problem is hung

The calculation will generate an output file, OUTP, and a continue-run file, RUNTPE. If the problem was terminated by <ctrl-c><q>, it can be continued by

MCNP6 c runtpe=RUNTPE

The output file from the continued run would be OUTQ (see section 1.4.1 in the MCNP manual for file naming).

The problem and its output files may also be renamed. For example, the initial run could be

MCNP6 i=EX11_1 n=XYZ.

Then the output file would be XYZ.o and the runtpe would be XYZ.r; the continuation would be

MCNP6 c runtpe=XYZ.r outp=ABC.o

In this example, the runtpe XYZ.r would contain the information from the first run and add, not overwrite, the information from the continuation. The continue-run output would be ABC.o.

A number of messages are printed by MCNP to the command window.

comment. Physics models disabled. comment. using random number generator 1, initial seed = 19073486328125 comment. total nubar used if fissionable isotopes are present.

Comments generally state what defaults are being used and usually can be ignored.

warning. there are no tallies in this problem. warning. no cross-section tables are called for in this problem.

Warnings indicate something that looks wrong but may be intentional.

Fatal errors are user errors found while processing the input file. After the first fatal error, subsequent ones may be spurious because they are caused by an earlier fatal error.

Bad trouble errors are issued while histories are being run and indicate imminent catastrophe. Rather than crash and lose all output, etc., the bad trouble error immediately stops the calculation and puts out what is available so far.

The OUTP output file of this simple problem does not contain much because the problem is so simple. After the copyright notice is

1mcnp version 6 ld=06/09/17 05/15/19 13:31:20 **************************** probid = 05/15/19 13:31:20 i=EX11_1 n=xyz.

The version and load date (ld=06/09/17) uniquely identify the version of MCNP being run followed by the time of the run. The probid (Problem Identification) is the unique identifier of the run, which is the time the run was started. The probid will be attached to all subsequent files for identification with this initial run. Next is the echo of the MCNP command. Here the input file name was EX11_1.

Following is a reprint of the input file and various print tables, which can be turned on and off with a PRINT card in the input file. The input file of Example 2.1 had no PRINT card, so only the minimum printout is provided. The more general print contains information about geometry, surfaces, sources, tallies, physics, data sources, variance reduction, and more.

The problem summary is the balance (creation and loss) of physical particles (weight), energy, and statistical samples (tracks). Absent from this input is what usually follows—namely summaries of all physics and statistical information in every problem cell and tally. The tallies are the answers requested by the user.

2.1.3 Simple Input File

The MCNP input file of Example 2.2 shows how to set materials and additional capabilities of the MCNP code to model a simple nuclear fuel rod. This input file describes a 300 cm long UO2 fuel rod of radius 0.4840 cm encased in zirconium cladding with outer radius 0.5590 cm and is illustrated in Fig. 2.1.

Fig. 2.1
figure 1

X – Y view of Example 2.2 geometry as plotted by the MCNP plotter. Labels 11–13 label the cell numbers and labels 71–72 the surface numbers, as specified in the EX11_2 input file

Example 2.2 A Simple Input File

Fuel Pin c **************************** cells **************************** 11 313 -10.44 -71 u=0 imp:n=1 $ Fuel 12 201 -6.55 -72 71 u=0 imp:n=1 $ Cladding 13 0 72 u=0 imp:n=0 $ Outside World c **************************** surfaces ************************* 71 RCC 0 0 -150 0 0 300 .4840 $ fuel cylinder 72 RCC 0 0 -150 0 0 300 .5590 $ clad cylinder c *************************** data ****************************** SDEF Print NPS 1e6 m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 nlib=80c m313 92235.80c -0.02759 92238.80c -0.85391 8016.80c -0.11850

There are three comment cards in the input file of Example 2.2. The first line of Example 2.2 input file is the title card; the title is “Fuel Pin.”

Cell cards other than comment cards have the following format: cell number, cell material, cell density, surface relationships, and data. Cell, surface, and material numbers are all arbitrary. Densities are positive for units of atom/(barn-cm) and negative for g/cm3.

11 313 -10.44 -71 u=0 imp:n=1 $ Fuel 12 201 -6.55 -72 71 u=0 imp:n=1 $ Cladding 13 0 72 u=0 imp:n=0 $ Outside World

Cell 11 has material 313, density 10.44 g/cm3, inside surface 71, and neutron importance imp:n=1. (The u is explained below.) Cell 12 has material 201, density 6.55 g/cm3, inside surface 72, and outside surface 71, with importance 1. Cell 13 has material 0, which is a void, in which case no density is specified. Cell 13 is outside surface 72 with importance 0. Note that the space outside surface 72 extends forever; all of space out to infinity must be specified. But because the importance is imp:n=0, any neutron that crosses surface 72 into cell 13 escapes the problem and is killed. Here, u=0 is a placeholder; it will be used to build fuel pin assemblies (see Sect. 2.1.7).

Surface cards start with a surface number, followed by a mnemonic surface type, and then the numerical entries for that surface type.

71 RCC 0 0 -150 0 0 300 .4840 $ fuel cylinder 72 RCC 0 0 -150 0 0 300 .5590 $ clad cylinder

Surface 71 is a right circular cylinder, RCC. The base is at x, y, z = 0, 0, −150 cm. The height vector is 0, 0, 300, corresponding to axis direction cosines u, v, w = 0, 0, 1, namely the z-axis with height 300 cm. The radius of surface 71 is 0.4840 cm and the radius of surface 72 is 0.5590 cm.

Thus, the problem geometry is two concentric cylinders with radii 0.4840 cm and 0.5590 cm. These form three cells: the fuel inside surface 71, the clad between surfaces 71 and 72, and the rest of the world outside surface 72.

After the cell and surface cards are the data cards:

SDEF Print NPS 1e6 m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 nlib=80c m313 92235.80c -0.02759 92238.80c -0.85391 8016.80c -0.11850

Data cards may appear in any order. SDEF is the default source definition. Print provides a full rather than partial print in the OUTP output file. NPS runs the problem for one million (1e6) source histories unless the problem is otherwise interrupted. The number, NPS, of source histories, often referred to as “source particles” by MCNP, is the number of times the problem source distribution will be sampled. The number of histories is statistical and has no physical meaning. Note that by default, SDEF has WGT = 1. The particle weight (WGT) represents the physical number of particles modeled per source particle and has no statistical meaning. Thus NPS 1e6 WGT = 1.0 samples the behavior of 1.0 physical source particle a million times; NPS 1000 WGT = 6.02E23 samples the behavior of 6.02E23 physical source particles a thousand times.

Material 201 is specified next and is continued on the following line. Material 313 follows. Materials and cross sections are described in Sect. 2.2.

2.1.4 Running and Plotting MCNP Geometries

Input file EX11_2 of Example 2.2 can be run as

MCNP6 i=EX11_2 n=problem02. IPXRZ

This problem will generate output file problem02.o and continue-run file problem02.r, which are more readily associated with input file EX11_2 than the default OUTP and RUNTPE names. It is recommended that the name option, name= or n=, is used to give more meaningful file names.

IPXRZ is the execution line option: process input I, plot geometry P, get the cross-section data X, run the problem R, and then plot the tallies Z. Usually it is impractical to use all five of these options together. IXR is the default. IP can be used to simply plot the geometry. The X–Y view MCNP plot of Example 2.2 geometry is shown in Fig. 2.1.

To get the plot in Fig. 2.1, run MCNP with

MCNP6 i=EX11_2 n=plot01. IP

When the plot screen pops up, click the following buttons:

XY

Will get X–Y view

L2 Off

Will plot cell labels

10

Upper right corner—will zoom in by a factor of 10

10

Upper right corner—will zoom in by another factor of 10

Alternatively, when the plot screen pops up, click Plot> bottom center, and then enter the following commands after the plot> prompt that appears in the same window where the MCNP code was run:

PZ=0 EX=1 LA=2 2 CEL

Note that the cell labels are twice as large because of the LA 2 2 CEL command; entering plot commands rather than clicking on the plot enables more capabilities. While in the plot> mode, the commands HELP, ?, and OPTIONS provide a help package that lists all possible geometry plot commands. To return to the interactive mode, enter INTERACT. To exit, click END in the interactive mode; type END in the plot> mode.

Plot> commands, such as EX=1, may also be entered at the lower left of the screen where it says “Click here or picture or menu.” If this box is accidentally clicked, then everything is frozen until a command is entered there.

The other commands in the lower left are:

CURSOR

Select a part of the picture in which to zoom in

RESTORE

Return to previous picture

CellLine

Outline cells or show various tally and weight window meshes

PostScript

Make a PostScript format file of the picture

ROTATE

Rotate picture 90 degrees counterclockwise; same as “Theta 90” Plot> command

COLOR

Turn color on or off; click a cell quantity on the right margin, such as DEN, and colors will be by density rather than material

SCALES

Geometric scale in cm: 0 = none, 1 = on boundary; 2 = everywhere (grid)

LEVEL

Repeated structures universe level to view

XY, YZ, ZX

Plot coordinate axis view

L1 sur

On/off surface labels

L2 cel

On/off cell labels. Click on right to label by cel, mat, den, etc.

MBODY

On/off to plot macrobody facet numbers

LEGEND

Numerical values of shaded colors

The region in the upper left is not interactive, but if the screen goes blank this region can be clicked to restore the previous picture.

The controls at the top shift the screen up, down, right, or left. “Origin” moves the center of the plot screen to whichever part of the geometry is clicked. The scale from 0.1 to 10, top right, is a logarithmic scale zoom factor. Click it twice to zoom by that factor. Click it once and then click somewhere in the picture to zoom in or out of the clicked point in the picture.

Plotting is described in more detail in Sect. 2.5.1 and also in the MCNP6.2 code manual under the heading “5.2 THE GEOMETRY PLOTTER, PLOT.”

The output file, problem01.o, like Example 2.1, begins with copyright, code version, date/time, probid, execution line, and input file echo. But because it is no longer a void problem but one with neutron transport physics, and because of the PRINT card, there is much more. The source, materials, volume, masses, surface areas, and surface descriptions are all listed. The cross sections for each material are described in Print Table 100 which appears in every output file if cross sections are used. It is imperative to review Print Table 100 to ensure that the cross sections being used are the ones desired—in some circumstances, the MCNP code will select others!

Print Table 110 lists the first 50 source particles, which helps determine if the source is set up as desired. The problem summary shows the gains and losses of the 1,000,000 statistical tracks and the physics of the average source particle. Additional tables show fission multiplicity, activity in each cell, and activity in each material.

2.1.5 Surfaces and Complicated Cells: Intersections and Unions

Example 2.3 and Example 2.4 describe a 12 cm high radius 5 cylinder inside a 1 cm thick box illustrated in Figs. 2.2 and 2.3. A million statistical histories are run with an isotropic 1 MeV source at x, y, z = 0, 6, 0 cm.

Fig. 2.2
figure 2

Plot of macrobody geometry of Example 2.3

Fig. 2.3
figure 3

Plot of quadratic surface geometry of Example 2.4

Example 2.3 Macrobody Description of Cylinder in a 1 cm Thick Box

Cylinder inside 1-cm thick box 1 0 -1 imp:n = 1 $ cylinder 2 0 1 -3 imp:n = 1 $ inner box 3 0 3 -4 imp:n = 1 $ outer box 4 0 4 imp:n = 0 $ Outside world 1 rcc 0 0 0 0 12 0 5 $ Right circular cylinder, 12 high, radius 5 3 rpp -6 6 0 14 -6 6 $ Right parallelepiped -6<x<6 0<y<14 -6<z<6 4 rpp -7 7 -1 15 -7 7 $ Right parallelepiped -7<x<7 -1<y<15 -7<z<7 sdef pos 0 6 0 erg = 1 nps 1000000

Example 2.4 Quadratic Surface Description of Cylinder in a 1 cm Thick Box

Cylinder inside 1-cm thick box 1 0 -1 2 -3 imp:n = 1 $ cylinder 2 0 (1:-2:3) 5 -6 2 -7 8 -9 imp:n = 1 $ inner box 3 0 (-5:6:-2:7:-8:9) 10 -11 12 -13 14 -15 imp:n = 1 $ outer box 4 0 -10:11:-12:13:-14:15 imp:n = 0 $ Outside world 1 cy 5.0 2 py 0 3 py 12 c inner box 5 px -6 6 px 6 7 py 14 8 pz -6 9 pz 6 c outer box 10 px -7 11 px 7 12 py -1 13 py 15 14 pz -7 15 pz 7 sdef pos 0 6 0 erg = 1 nps 1000000

MCNP surfaces may be described by either macrobodies or quadratic surface equations. In a cell description, blanks between the surface numbers are the intersection operator, “AND,” whereas the union operator, “:”, between the surface numbers is “OR.” Cells using the union operator are “complicated cells.” Whereas macrobodies are internally converted to quadratic surfaces, cells using macrobody surfaces are also complicated cells.

The quadratic surface equations are described in Table 3-4 of the MCNP manual [2, 3]. They consist of planes P, PX, PY, PZ (arbitrary, normal to X, normal to Y, and normal to Z axis); spheres SO, S, SX, SY, SZ (centered at origin, anywhere, or along X, Y, or Z axis); cylinders C/X, C/Y, C/Z, CX, CY, CZ (parallel or on X, Y, Z axis); cones K/X, K/Y, K/Z, KX, KY, KZ (parallel or on X, Y, or Z axis); special quadratic SQ (ellipsoid, hyperboloid, paraboloid parallel to X, Y, or Z axis); general quadratic (cylinder, cone, ellipsoid, hyperboloid, paraboloid with arbitrary axis); and tori TX, TY, TZ (elliptical or circular torus on X, Y, or Z axis).

Surfaces may also be defined by coordinate points. Quadratic surfaces rotationally symmetric about the major axes are specified by the X, Y, and Z surface cards in (x,R) pairs, where x is the x, y, or z coordinate and R is the radius. One pair makes a PX, PY, or PZ plane; two pairs make a CX, CY, CZ cylinder or KX, KY, KZ cone; three pairs make an ellipsoid, paraboloid, or hyperboloid. Skew planes may be specified on the general plane P surface card by providing three triplets of x, y, z coordinates.

It is easier and more common to describe MCNP surfaces by macrobodies described in Section 3.2.2.4 of the MCNP manual. These are BOX (planes specified by the three vectors at a corner); rectangular parallelepiped RPP (box of planes normal to X, Y, Z axes); sphere SPH (specifying X, Y, Z center and radius); right circular cylinder RCC (specifying X, Y, Z base and h1, h2, h3 height vector, and radius); and right hexagonal prism RHP, right elliptical cylinder REC, truncated right cone TRC, ellipsoid ELL, wedge WED, and arbitrary polyhedron ARB. The macrobodies are internally converted by the MCNP code into quadratic surface facets. Figure 2.2 shows the facets of macrobody surfaces 1, 3, and 4. The MCNP geometry plotter shows these by setting the MBODY option off.

The MCNP cell description uses senses to indicate whether a cell is inside or outside a surface. In the quadratic surface description of Example 2.4:

1 0 -1 2 -3 imp:n = 1 $ cylinder

Cell 1 with material zero (vacuum) is inside y-axis cylinder surface 1 cy 5.0, outside (above) y-axis plane surface 2 py 0, and below y-axis plane surface 3 py 12.

2 0 (1:-2:3) 5 -6 2 -7 8 -9 imp:n = 1 $ inner box

Cell 2 uses the union operator “:” to describe everything outside 1 or below 2 or above 3 (1:−2:3) combined with what is above 5 and below 6 and above 2 and below 7 and above 8 and below 9.

3 0 (-5:6:-2:7:-8:9) 10 -11 12 -13 14 -15 imp:n = 1 $ outer box

Similarly, cell 3 is everything below 5 or above 6 or below 2 or above 7 or below 8 or above 9, combined with what is above 10 and below 11 and above 12 and below 13 and above 14 and below 15.

4 0 -10:11:-12:13:-14:15 imp:n = 0 $ Outside world

Cell 4 is the outside world extending to infinity and is everything below 10 or above 11 or below 12 or above 13 or below 14 or above 15.

The sense of a surface to a cell can be determined by picking any point in the cell, x, y, z and then using that point to evaluate the sense to the surface. For example, the z-axis CZ cylinder has the equation x2 + y2 – R2 = 0; for radius R = 10 and x, y, z = 1,2,3, then 12 + 22 – 100 = −95 < 0, so the cell containing this point has a negative sense relative to the surface; for R = 10 and x, y, z = 8, 9, 0, then 82 + 92 – 100 = 45 > 0, so the cell containing this point has a positive sense relative to the surface. (Surface equations are given in the MCNP manual 3.2.2.1.)

2.1.6 Duplicate Cells, Complements, and Translations: LIKE n BUT and TRCL

Example 2.5 and Example 2.6 take the geometry of Example 2.3 and Example 2.4 and duplicate it in a 30° rotation translated 20 cm away in the X-direction. The source is now a 1 MeV isotropic point emitting at x, y, z = 0, 6, 0 and x, y, z =20, 6, 0. Cells 11, 12, and 13 are the translated and rotated cells. See Figs. 2.4 and 2.5.

Fig. 2.4
figure 4

Cylinder inside a 1 cm thick box duplicated into a second rotated box using LIKE BUT and cell rotation/translation (Example 2.5)

Fig. 2.5
figure 5

Cylinder inside a 1 cm thick box duplicated into a second rotated box specifying 2nd box and using surface rotation/translation (Example 2.6)

Example 2.5 Cylinder Inside a 1 cm Thick Box Duplicated into a Second Rotated Box Using LIKE BUT and Cell Rotation/Translation

Cylinders inside 1-cm thick box 1 0 -1 imp:n = 1 $ cylinder 2 0 1 -3 imp:n = 1 $ inner box 3 0 3 -4 imp:n = 1 $ outer box 4 0 4 #11 #12 #13 imp:n = 0 $ Outside world 11 like 1 but *TRCL=(20 0 0 30 90 120 90 0 90 -60 90 30) $ 2nd cylinder 12 like 2 but TRCL=100 $ 2nd inner box 13 like 3 but TRCL=100 $ 2nd outer box 1 rcc 0 0 0 0 12 0 5 $ Right circular cylinder, 12 high, radius 5 3 rpp -6 6 0 14 -6 6 $ Right parallelepiped -6<x<6 0<y<14 -6<z<6 4 rpp -7 7 -1 15 -7 7 $ Right parallelepiped -7<x<7 -1<y<15 -7<z<7 sdef pos=d1 erg = 1 si1 L 0 6 0 20 6 0 sp1 1 1 nps 1000000 *TR100 20 0 0 30 90 120 90 0 90 -60 90 30

Example 2.6 Cylinder Inside a 1 cm Thick Box Duplicated into a Second Rotated Box Specifying the 2nd Box and Using Surface Rotation/Translation

Cylinders inside 1-cm thick box

1 0 -1 imp:n = 1 $ cylinder

2 0 1 -3 imp:n = 1 $ inner box

3 0 3 -4 imp:n = 1 $ outer box

21 0 -11 imp:n = 1 $ 2nd cylinder

22 0 11 -13 imp:n = 1 $ 2nd inner box

23 0 13 -14 imp:n = 1 $ 2nd outer box

4 0 4 14 imp:n = 0 $ Outside world

1 rcc 0 0 0 0 12 0 5 $ Right circular cylinder, 12 high, radius 5

3 rpp -6 6 0 14 -6 6 $ Right parallelepiped -6<x<6 0<y<14 -6<z<6

4 rpp -7 7 -1 15 -7 7 $ Right parallelepiped -7<x<7 -1<y<15 -7<z<7

11 100 rcc 0 0 0 0 12 0 5 $ Right circular cylinder, 12 high, radius 5

13 100 rpp -6 6 0 14 -6 6 $ Right parallelepiped -6<x<6 0<y<14 -6<z<6

14 100 rpp -7 7 -1 15 -7 7 $ Right parallelepiped -7<x<7 -1<y<15 -7<z<7

sdef pos=d1 erg = 1

si1 L 0 6 0 20 6 0

sp1 1 1

nps 1000000

*TR100 20 0 0 30 90 120 90 0 90 -60 90 30

In Example 2.5, cells 11, 12, and 13 are like cells 1, 2, and 3 but translated and rotated. The rotation translation is

O1 O2 O3 XX’ YX’ ZX’ XY’ YY’ ZY’ XZ’ YZ’ ZZ’ m

where X, Y, Z are the Cartesian coordinate axis vectors and X′,Y′,Z′ are the rotated axis vectors. If m = −1, then the coordinate systems are reversed. Thus

*TR100 20 0 0 30 90 120 90 0 90 -60 90 30

is a 30° rotation in the XZ plane translated 20 cm in the X-direction. TR and TRCL are specified in cosine; *TR and *TRCL are specified in degrees. The transformation can be either on the cell card, TRCL=(…) as in cell 11, or on a TR card in the data section of the input and specified on the cell card by TRCL=100, as in cells 12 and 13.

4 0 4 #11 #12 #13 imp:n = 0 $ Outside world

The complement operator, #, is now required in cell 4 to indicate that LIKE BUT cells 11, 12, and 13 are excluded from the outside world in addition to everything outside surface 4.

When a cell is translated, it must create new surfaces. In Fig. 2.4, surfaces 11001, 12003, and 13004 are added for surface 1 of cell 11, surface 3 of cell 12, and surface 4 of cell 13. Like macrobody surfaces 1, 2, and 3, these surfaces have facets. Note that the Y-axis planes in cells 11, 12, and 13 are the same as those in cells 1, 2, and 3. Consequently the Y-axis facets, 1.2, 1.3, 3.3, 3.4, 4.3, and 4.4 are used in both cells 1, 2, 3 and 11, 12, 13. New facets are created only for the X-axis and Z-axis facets. The complement, #, operator is much more convenient than using the unions of the new and old named facets in the description of cell 4.

Example 2.6 duplicates the cylinder in a 1 cm thick box by simply adding cells 21, 22, 23 and surfaces 11, 13, and 14. There is no LIKE BUT or cell rotation/ translation, but now there is surface rotation/translation. Surfaces 11, 12, and 14 refer to transformation TR100 by putting the transformation number after the surface number.

2.1.7 Filled Cells: Universes

Example 2.7 illustrates universe/fill geometries. The cylinder and region outside it from Example 2.3/Example 2.4 are made into their own universe by adding u = 1 to cells 1, 2, and 3. Figure 2.6 shows the full geometry. It is composed of the “real world” universe, u = 0, geometry illustrated in Fig. 2.7. Note that there is no longer cell 4, the outside world, because universe 1 cell 3 extends out to infinity—namely everything outside surface 3, as illustrated in Fig. 2.8. Instead, cell 14 is added, which is bounded by everything inside surface 4 and filled with universe 1. Cell 11 is also filled with universe 1. Cell 15 is now the outside world.

Fig. 2.6
figure 6

MCNP plot of Example 2.7, cylinder inside a 1 cm thick box duplicated into a second rotated box by universe fill

Fig. 2.7
figure 7

MCNP plot of Example 2.7 level 0 (universe u = 0)

Fig. 2.8
figure 8

MCNP plot of Example 2.7 level 1 (universe u = 1)

Example 2.7 Cylinder Inside a 1 cm Thick Box Duplicated into a Second Rotated Box by Universe Fill

Cylinders inside 1-cm thick box 1 0 -1 u=1 imp:n=1 $ cylinder 2 0 1 -3 u=1 imp:n=1 $ inner box 3 0 3 u=1 imp:n=1 $ outside inner box 14 0 -4 fill=1 imp:n=1 $ outer box 11 like 14 but TRCL=100 $ 2nd cylinder in box 15 0 #14 #11 imp:n=0 $ outside world 1 rcc 0 0 0 0 12 0 5 $ Right circular cylinder, 12 high, radius 5 3 rpp -6 6 0 14 -6 6 $ Right parallelepiped -6<x<6 0<y<14 -6<z<6 4 rpp -7 7 -1 15 -7 7 $ Right parallelepiped -7<x<7 -1<y<15 -7<z<7 sdef pos=0 6 0 erg = 1 cel=d2 si2 L (1<14) (1<11) sp2 1 1 nps 1000000 PRINT *TR100 20 0 0 30 90 120 90 0 90 -60 90 30

Cell 11 is the outer box translated 20 cm in the X-axis and tilted 30°:

14 0 -4 fill 1 imp:n=1 $ outer box 11 like 14 but TRCL=100 $ 2nd cylinder in box

The outside world can be specified either of two ways:

15 0 #14 #11 imp:n=0 $ outside world

or

15 0 4 #11 imp:n=0 $ outside world

Cells 11, 14, and 15 in Fig. 2.6 comprise universe 0, “the real world,” and u = 0 on the cell cards is optional. Note how the surfaces of LIKE n BUT cell 11 have been renamed 11004 for facets 1,2,5, and 6. The output file describes what happens to Y-axis top and bottom facets 3 and 4:

surface 1.3 and surface 3.4 are the same. 3.4 will be deleted. surface 4.3 and surface 11004.3 are the same. 11004.3 will be deleted. surface 4.4 and surface 11004.4 are the same. 11004.4 will be deleted.

The filled cells, 1, 2, and 3, retain their names, surface numbers, and coordinates. Consequently, the source must be changed to a distribution of cells:

sdef pos=0 6 0 erg = 1 cel=d2 si2 L (1<14) (1<11) sp2 1 1

The source location, pos=0 6 0, is the same for cell 1 level 1 whether it is in the universe below cell 11 or cell 14 (level 0). The first 50 histories, Print Table 100 in the output file, show the source starting in cell 1 at x, y, z = 0, 6, 0 all the time. The source is in either cell 11 or 14 about half the time and is at x, y, z = 20, 6, 0 in cell 11 and at x, y, z at 0, 6, 0 in cell 14.

nps x y z cell surf u v w 1 2.000E+01 6.000E+00 0.000E+00 11 8.004E-01 5.963E-01 6.178E-02 0.000E+00 6.000E+00 0.000E+00 1 0 6.623E-01 5.963E-01 4.537E-01 2 0.000E+00 6.000E+00 0.000E+00 14 -2.775E-01 -1.028E-01 9.552E-01 0.000E+00 6.000E+00 0.000E+00 1 0 -2.775E-01 -1.028E-01 9.552E-01

Print Table 128 in the output file shows the activity in each cell and its path:

1neutron activity in each repeated structure / lattice element print table 128 source entering collisions path 500488 500488 0 1 < 14 0 442581 0 2 < 14 0 500488 0 3 < 14 499512 499512 0 1 < 11 0 441725 0 2 < 11 0 499512 0 3 < 11 0 0 0 15 1000000 2884306 0

Print Table 128 is sometimes turned off even when the PRINT card is in the input deck. Print Table 128 can require excessive amounts of computer time and storage, locating the path to thousands of elements at lower universe levels such as in the case of a nuclear reactor model or human phantom model consisting of millions of pixels. When the MCNP code turns off Print Table 128, the following warning is issued:

warning. universe map (print table 128) disabled.

2.1.8 Lattice Geometries

Example 2.8 illustrates square and hexagonal lattices and is shown in Fig. 2.9. The geometry consists of an RPP box (cell 14) filled with universe 11 and another one like it but translated 20 cm in the X-direction and filled with universe 12, as shown in Fig. 2.10. Both are contained in a huge sphere, surface 20.

Fig. 2.9
figure 9

X – Z view of Example 2.8 with both square and hex lattices

Fig. 2.10
figure 10

Example 2.8 level u = 0 universe. Note how surface facets 11004.1 and 11004.2 are created because cell 11 is like cell 14 but surfaces 4.1 and 4.2 remain with cell 14

Example 2.8 Square and Hex Lattice File

Square and hex lattices 1 313 -10.44 -1 u=1 imp:n=1 $ UO2 2 201 -6.55 -2 1 u=1 imp:n=1 $ Zr clad 3 609 -1.00 2 u=1 imp:n=1 $ H2O 21 0 -31 lat=1 fill=1 u=11 imp:n=1 $ Square lattice 22 0 -32 lat=2 fill=1 u=12 imp:n=1 $ Hexagonal lattice 14 0 -4 fill 11 imp:n=1 $ Square lattice box 11 like 14 but TRCL=(20 0 0) fill 12 $ Hexagonal lattice box 15 0 #14 #11 -20 imp:n=1 $ outside boxes 16 0 20 imp:n=0 $ outside world 1 rcc 0 0 0 0 12 0 .5 $ Right circular cylinder 2 rcc 0 0 0 0 12 0 .6 $ Right circular cylinder 4 rpp -7 7 0 12 -7 7 $ Right parallelepiped 14x12x14 31 rpp -1 1 0 12 -1 1 $ Right parallelepiped 2x12x2 32 rhp 0 0 0 0 12 0 1 0 0 $ Right hex 12x1 20 sph 0 0 0 50 $ Sphere containing lattice boxes m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 nlib=80c m313 92235.80c -0.02759 92238.80c -0.85391 8016.80c -0.11850 M609 1001.80c 2 8016.80c 1 MT609 LWTR sdef pos=10 6 0 nps 1000000 PRINT

At universe level 1, universe 11 fills cell 14 and universe 12 fills cell, as illustrated in Fig. 2.11. Cell 21 of universe 11 is a square lattice, LAT = 1, and cell 22 of universe 12 is a hexagonal lattice, LAT = 2.

Fig. 2.11
figure 11

Example 2.8 level 1 universe showing LAT = 1 rectangular lattice filling cell 21 and LAT = 2 hexagonal lattice filling cell 22

Figure 2.12 illustrates universe level 2, where both the rectangular and hexagonal lattices are filled with universe 1. Universe 1 is a UO2 fuel rod (cell 1, material M313) in Zr cladding (cell 2, material M201) in an infinite extent of water (cell 3, material 609) outside surface 2.

Fig. 2.12
figure 12

Example 2.8 level 2 universe showing cells 1, 2, and 3 repeated in the lattices

The fuel rods can be differentiated according to which lattice element they fill. The indices are plotted in Fig. 2.13 using the LA=0 1 IJK command. With the interactive plotter, click IJK on the right margin and then label L2 in the lower left menu box.

Fig. 2.13
figure 13

Example 2.8 level 2 universe showing cells 1, 2, and 3 repeated in the lattices and labeled by lattice index

The source of Example 2.8 is a 14 MeV isotropic neutron (default) midway between the two lattice boxes at POS = 10 6 0. From the problem summary in the output file, 0.13 fission neutrons and 0.07 (n,xn) neutrons are produced per source particle.

Print Table 128 in the output file shows the activity in each cell and its path:

1neutron activity in each repeated structure / lattice element print table 128 source entering collisions path 0 708051 216431 1 < 21 < 14 0 1587248 67168 2 < 21 < 14 0 2736613 1330851 3 < 21 < 14 0 921642 251038 1 < 22 < 11 0 1868093 79206 2 < 22 < 11 0 2972581 1227971 3 < 22 < 11 1000000 1801413 0 15 0 0 0 16 1000000 12595641 3172665

The lattices in Example 2.8 are simple lattices because all elements are filled with the same universe.

2.1.9 Fully Specified Lattice Geometries

A lattice in which the elements are filled with different universes is a “fully specified lattice.” In a reactor fuel assembly, these would be for water holes, fuel rods, or fuels of different compositions. In a 3He detector assembly, these would be 3He tubes of different atmospheres. In a human phantom model, the elements would be filled with different human organs. Example 2.9 is an MCNP input file with a fully specified lattice geometry.

Example 2.9 Fully Specified Square and Hex Lattices with Multiple Universes

Square and hex lattices 1 313 -10.44 -1 u=1 imp:n=1 $ UO2 2 201 -6.55 -2 1 u=1 imp:n=1 $ Zr clad 3 609 -1.00 2 u=1 imp:n=1 $ H2O 21 111 -4.00 -31 lat=1 u=11 imp:n=1 $ Square lattice fill = -3:3 0:0 -3:3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 1 1 1 1 1 1 22 111 -4.00 -32 lat=2 u=12 imp:n=1 $ Hexagonal fill = -5:5 -4:4 0:0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 12 1 1 1 1 1 1 1 1 1 1 14 0 -4 fill 11 imp:n=1 $ Square lattice box 11 like 14 but TRCL=(20 0 0) fill 12 $ Hexagonal lattice box 15 0 #14 #11 -20 imp:n=1 $ outside boxes 16 0 20 imp:n=0 $ outside world 1 rcc 0 0 0 0 12 0 .5 $ Right circular cylinder, 12h,.5r 2 rcc 0 0 0 0 12 0 .6 $ Right circular cylinder, 12h,.6r 4 rpp -7 7 0 12 -7 7 $ Right parallelepiped 14x12x14 31 rpp -1 1 0 12 -1 1 $ Right parallelepiped 2x12x2 32 rhp 0 0 0 0 12 0 1 0 0 $ Right hex 12x1 20 sph 0 0 0 50 $ Sphere containing lattice boxes m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 nlib=80c m313 92235.80c -0.02759 92238.80c -0.85391 8016.80c -0.11850 M609 1001.80c 2 8016.80c 1 MT609 LWTR M111 6000 1 1001 1 7014 1 sdef pos=10 6 0 nps 1000000 PRINT

The simple lattice definition of Example 2.8 was

21 0 -31 lat=1 fill=1 u=11 imp:n=1 $ Square lattice 22 0 -32 lat=2 fill=1 u=12 imp:n=1 $ Hexagonal lattice

To change one element the fully specified lattice specification is

21 111 -4.00 -31 lat=1 u=11 imp:n=1 $ Square lattice fill = -3:3 0:0 -3:3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 1 1 1 1 1 1 22 111 -4.00 -32 lat=2 u=12 imp:n=1 $ Hexagonal lattice fill = -5:5 -4:4 0:0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 12 1 1 1 1 1 1 1 1 1 1

A new material, M111, is added to cells 21 and 22 so that the different elements in Fig. 2.14 contain different universes and are plotted in a different color.

M111 6000 1 1001 1 7014 1

Fig. 2.14
figure 14

Fully specified rectangular and hex lattices of Example 2.9, with one element filled differently in each

Instead of FILL=1 for both lattices, the FILL is replaced with the range of lattice elements needed to fill the filling cells 21 and 22:

fill = -3:3 0:0 -3:3 fill = -5:5 -4:4 0:0

In the square lattice, element [−3,0,3] is filled with universe 11. Because universe 11 is the same universe as its lattice cell, cell 21 u = 11, the lattice is “filled with itself” with material M111 instead of the fuel pins of universe 1. Note that this is the upper left in the MCNP plot but the lower left in the lattice map of the FILL entry.

In the hexagonal lattice, element [−5,4,0] is filled with universe 12. Because universe 12 is the same universe as its lattice cell, cell 22 u = 12, the lattice is “filled with itself” with material M111 instead of the fuel pins of universe 1. Note that unlike the rectangular lattice, the universe 12 element is the lower left in both the MCNP plot and in the lattice map of the FILL entry. In the hexagonal lattice there are two constraints: the first and second indices must increase across adjacent surfaces and the third index must increase in one or the other direction along the length of the hexagonal axis. A rectangular lattice can completely fit into a rectangular box, but a hexagonal lattice cannot completely fit in a rectangular box or any other MCNP microbody surface shape. Consequently, the hexagonal lattice of Fig. 2.14 is filled with 99 elements −5:5 −4:4 0:0, but only 45 are full hexes, 22 are cut in half, and the remaining 32 are cut off altogether. A further description of the MCNP LAT and FILL capability is in the MCNP Manual [2] Section 3.3.1.5.2 and 3.3.1.5.3.

The lattice elements are most easily recognizable from geometry plots that have the label LA=0 1 IJK or by clicking on ijk in the right margin and then L2 in the lower left menu. These indices are shown in Figs. 2.15 and 2.16.

Fig. 2.15
figure 15

Square lattice with ijk lattice indices. The extra (−1,0,3) in the upper left is an MCNP6.2 plotting bug that repeats the label of the upper right. The I indices increase left to right, which is the X-direction for RPP surface 31 of lattice cell 21. The J indices are from 0:0 and they increase in the Y-direction (out of the plane of the plot) for RPP surface 21 of lattice cell 21. The K indices increase from bottom to top, which is the X-direction for RPP surface 31 of lattice cell 21

Fig. 2.16
figure 16

Hexagonal lattice with ijk lattice indices. The extra (−2,4,0) in the lower left is an MCNP6.2 plotting bug that repeats the label of the lower right. The I indices increase at 60° from lower left to upper right. The J indices increase at 120° from upper left to lower right. The K indices increase from bottom to top of the hex axis—from 0:0 in this case—because the hex is only one element high

The fully specified hexagonal lattice is more commonly specified in the MCNP input file with the elements of the hexagonal lattice cut off by the filling cell (cell 11 filled with universe 12, which is the hex lattice cell) being filled with universe 0:

22 111 -4.00 -32 lat=2 u=12 imp:n=1 $ Hexagonal fill = -5:5 -4:4 0:0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 12 1 1 1 1 1 1 0 0 0 0

The above hexagonal lattice is also used later in Example 2.11, Example 2.12, and Example 2.19.

2.2 Materials and Cross Sections

2.2.1 Specifying Materials

An MCNP input file consists of three sections: cells, sources, and data. The data section of Example 2.8 is shown in Example 2.10:

Example 2.10 Data Cards from Example 2.6

m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 nlib=80c m313 92235.80c -0.02759 92238.80c -0.85391 8016.80c -0.11850 M609 1001.80c 2 8016.80c 1 MT609 LWTR sdef pos=10 6 0 nps 1000000 PRINT

Materials are specified on the materials cards, but additional input in the data section is often required. Material cards consist of ZAID-fraction pairs. Negative fractions are mass fractions; positive fractions are atom fractions. The fractions are normalized to unity for all materials. The ZAID numbers consist of XXXAAA.nnx, where ZZZ is the atomic number, AAA is the atomic mass, nn is the specific evaluation, and x is the data type.

Material numbers are arbitrary. For material m201, ZAID = 40090 indicates Z = 40 zirconium isotope A = 90. Thus, material 201 consists of five zirconium isotopes—90Zr, 91Zr, 92Zr, 94Zr, and 96Zr—with mass fractions 0.505239, 0.110180, 0.168413, 0.170672, and 0.027496. The optional nlib=80c sets nnx for each nuclide to 80c, which, in this case, indicates ENDF/B-VII.1 continuous-energy neutron cross sections.

Material 313 is specified as uranium dioxide, UO2, 3.17 at% or 3.13 wt% enriched in 235U.

2.2.2 Neutron Cross Sections

Cross sections are available only for common isotopes. If the cross section for an isotope specified on the material card does not exist, then the MCNP code will terminate with a fatal error. For example, if there is no tritium neutron transport cross section available,

M609 1003.80c 2 8016.80c 1

then the code terminates with the message

fatal error. cross-section tables missing for zaid = 1003.80c

Even worse, in problems with protons or other charged particles or ions, missing cross-section data will cause the MCNP code to use theoretical models for those data instead. The code will also use model physics for energies higher than the highest tabulated energy in any cross-section library it is using. These models are good for neutrons with energies >150 MeV and particles other than neutrons; the models are generally poor for neutrons with energies <150 MeV. It is important to check the cross sections used to ensure that they are in the range of particle energy simulated in each problem.

Further, different cross-section evaluations are available for different computers and different MCNP versions.

Wrong cross sections give wrong answers—the most frequent user error. There is a document listing the most available data [4], but these are not necessarily the ones available to the MCNP code on a particular computer and code installation. Consequently, the MCNP code frequently selects unexpected cross sections; therefore, it is imperative to always check in the MCNP OUTP output files to see which cross sections are actually being used by the code. The cross sections used are described in Print Table 100:

40090.80c 96741 Zr90 ENDF71x (jlconlin) Ref. see jlconlin (ref 09/10/2012) mat4025 12/18/12 Energy range: 1.00000E-11 to 2.00000E+01 MeV. particle-production data for protons being expunged from 40090.80c particle-production data for deuterons being expunged from 40090.80c particle-production data for tritons being expunged from 40090.80c particle-production data for helions being expunged from 40090.80c particle-production data for alphas being expunged from 40090.80c probability tables used from 5.3500E-02 to 1.0000E-01 mev.

The cross section used in this case is 90Zr from an experimental ENDF/B-VII library, ENDF71x, and in the energy range up to 20 MeV. Whereas the problem runs only neutrons, the production data for protons, deuterons, tritons, helions (3He), and alphas are deleted to save storage; these reactions are treated as capture. The probability tables from 53.5 to 100 keV provide better physics in the unresolved cross-section energy range. Often the temperature of the cross-section data is also provided.

Figure 2.17 is an MCNP plot of material 313. The MCNP code is run in the cross-section plotting mode:

MCNP6 I=Ex11_6a name=junk. IXZ

and the MCNP plot command is

xs m313 cop xs 92238.80c cop xs 92235.80c

Fig. 2.17
figure 17

Neutron cross section of material M313, UO2 (black), 238U (blue), and 235U (red) in Example 2.10

MCNP cross-section plotting is further described in Sect. 2.5. The cross-section table starts at 1e−11 MeV (1e−5 eV); the data are horizontally extrapolated at lower energies, though all neutrons are essentially absorbed before their energies become that low. From 1e−11 to 1e−5 MeV the neutron cross-section has a 1/v tail, where v is the neutron velocity. Above that is the resonance cross-section region up to 2 keV for 235U and 20 keV for 238U. Above that is the unresolved resonance cross-section region—where the resonances are so close together that they cannot be resolved and are represented by either a smooth curve or probability tables. Above 500 keV are inelastic and reaction cross sections such as secondary particle production. The uranium cross sections end at 20 MeV, and physics models are used above, or the cross section is merely extrapolated. The material m313 cross section goes up to 150 MeV because of the oxygen cross section.

2.2.3 Low-Energy Neutron Problems: Thermal Free Gas Treatment

Neutrons at very low energies and all collisions with hydrogen and very low-Z nuclides are affected by the thermal motion of the nuclei with which they interact. MCNP treats the kinematics of neutron scatter with moving nuclei with the thermal free gas treatment.

The thermal free gas treatment assumes room temperature. For problems at other temperatures, the temperature, kT, in MeV units must be specified.

kT (MeV)

Unit of T

8.617 × 10−11 T

Kelvin

8.617 × 10−11 (T + 273.15)

degrees C

4.787 × 10−11 T

degrees R

4.787 × 10−11 (T + 459.67)

degrees F

The temperature is usually specified on cell cards (using TMPn cell parameter). For cells 1 and 2 in Example 2.10, assume the temperatures are 300 °C and 600 °C. Then

1 313 -10.44 -1 u=1 imp:n=1 $ UO2 2 201 -6.55 -2 1 u=1 imp:n=1 $ Zr clad

becomes

1 313 -10.44 -1 u=1 tmp1=4.939e-8 imp:n=1 $ UO2 2 201 -6.55 -2 1 u=1 tmp1=7.524e-8 imp:n=1 $ Zr clad

Room temperature (20.5 °C = 2.5301e−8 MeV) is assumed for all cells without a TMP1 entry. Temperatures may be time-dependent using the THTME card and TMP1, TMP2, TMP3, etc.

Alternatively, temperatures may be specified in the data section of input for the first two cells among the cell cards and using the default (room temperature) for the other 7 cells in the problem.

TMP1 3.216e-8 4.939e-8 7J

Temperatures may also be time-dependent:

THTME 1e8 6e9 TMP1 3.216e-8 4.939e-8 7J TMP2 4.939e-8 7R

The above example has 100 °C and 300 °C for the first two cells and room temperature for the remaining seven cells from time 0 to 1 s (1e8 shakes). The temperature is 300 °C for all cells at 60 s (6e9 shakes) and later. The temperature is linearly interpolated from 1 s < time < 60 s. Note the use of J (jump) notation to take the default values and the R (repeat) notation to repeat the temperatures at 60 s.

Cell temperatures impact elastic scattering cross sections and collisions kinematics. The temperature on the TMP card does not affect inelastic or reactions cross sections, thermal scattering kernels, or resonances in scattering. Cross-section resonance broadening is done only by choosing cross-section data files at the desired temperature.

2.2.4 Low-Energy Neutron Problem Data: S(α,β) Thermal Treatment

Molecular binding of nuclides into molecules also affects the scattering of neutrons at low energies (few eV). The MCNP code treats these scatters with the S(α,β) thermal treatment [2]. By default, the S(α,β) treatment is off because the code does not know the compositions of molecules. The MCNP input must explicitly turn it on. If low-energy neutrons are important (almost always the case if fission is present), failure to specify S(α,β) will result in wrong answers. There is almost no computer time or storage cost for including S(α,β), so it is recommended to specify S(α,β) for all neutron problems where the limited amount of S(α,β) data is available.

M609 1001 2 8016 1 MT609 LWTR

As illustrated in Example 2.10, material M609 is light water, so its S(α,β) molecule is specified on the above MT609 card.

S(α,β) data are available for different temperatures of a limited number of molecules of hydrogen, deuterium, beryllium, graphite, oxygen, polyethylene, benzene, and more. The most reliable way to tell what is available is to examine the XSDIR data-directory file that comes with each version of the MCNP code and then to verify what the MCNP code is actually used by examining the output file. Further documentation is also available [2, 4].

As another example, consider material 111 from Example 2.7:

M111 6000.80c 1 1001.80c 1 7014.80c 1 MT111 POLY.10T BENZ.10T GRPH.10T

By adding the MT111 card to M111 in Example 2.7, the nitrogen, 7014.80c, will have no S(α,β) because none of the MT111 molecules contains nitrogen data. The hydrogen, 1001.80c, will use POLY.10T for S(α,β) because POLY appears before BENZ on the MT111 card and is the first MT111 molecule containing hydrogen data. Carbon, 6000.80c, will use GRPH.10T because POLY and BENZ apply only to the hydrogen in POLY and BENZ, not the carbon; GRPH is the only MT111 molecule containing carbon data. The carbon, 6000.80c, and graphite S(α,β), GRPH.10T, are plotted in Example 2.11.

Figure 2.18 is an MCNP plot of material M111 in Example 2.7. MCNP is run in the cross-section plotting mode:

MCNP6 I=Ex12_1b name=junk. IXZ

Fig. 2.18
figure 18

Carbon 6000.80c continuous-energy total (MT = −1) neutron cross section (black) co-plotted with GRPH.10T S(α,β) cross-section total (MT = 1) (blue), elastic (MT = 2) (red), and inelastic (MT = 4) (green)

The MCNP plot command is

xlim 1e-11 1e-5 xs 6000.80c mt 1 cop xs grph.10t mt 1 cop mt 2 cop mt 4

MCNP cross-section plotting is further described in Sect. 2.5. The S(α,β) cross sections extend from about 1e−11 MeV to the S(α,β) cutoff, generally about 1e−5 MeV. Above the cutoff, only the continuous-energy data are used. Below the cutoff, only the S(α,β) data are used. S(α,β) data should be used whenever possible. When neutrons do not scatter below the cutoff, the S(α,β) data are ignored. When neutrons do scatter below the cutoff, then the molecular thermal effects are important.

2.2.5 Photon Cross Sections

The MODE card is required whenever anything other than a neutron is transported. In Example 2.11, the M0 causes all neutrons to use the .70c ENDF/B-VII data and all photons to use the .04p photon data. If NLIB or PLIB were missing, then MCNP would use the default data libraries, which is whatever libraries MCNP finds first on the particular computer for the particular version of MCNP. Generally the default photon cross sections are good enough, but the output file should always be checked to see if the ones selected are appropriate.

Example 2.11 Photon Problem

Square and hex lattices 1 313 -10.44 -1 u=1 imp:n=1 $ UO2 2 201 -6.55 -2 1 u=1 imp:n=1 $ Zr clad 3 609 -1.00 2 u=1 imp:n=1 $ H2O 21 111 -4.00 -31 lat=1 u=11 imp:n=1 $ Square lattice fill = -3:3 0:0 -3:3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 1 1 1 1 1 1 22 111 -4.00 -32 lat=2 u=12 imp:n=1 $ Hexagonal fill = -5:5 -4:4 0:0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 12 1 1 1 1 1 1 0 0 0 0 14 0 -4 fill 11 imp:n=1 $ Square lattice box 11 like 14 but TRCL=(20 0 0) fill 12 $ Hexagonal lattice box 15 0 #14 #11 -20 imp:n=1 $ outside boxes 16 0 20 imp:n=0 $ outside world 1 rcc 0 0 0 0 12 0 .5 $ Right circular cylinder,12h,.5r 2 rcc 0 0 0 0 12 0 .6 $ Right circular cylinder,12h,.6r 4 rpp -7 7 0 12 -7 7 $ Right parallelepiped 14x12x14 31 rpp -1 1 0 12 -1 1 $ Right parallelepiped 2x12x2 32 rhp 0 0 0 0 12 0 1 0 0 $ Right hex 12x1 20 sph 0 0 0 50 $ Sphere containing lattice boxes MODE N P m0 nlib=.70c plib=.04p m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 nlib=80c m313 92235 -0.02759 92238 -0.85391 8016 -0.11850 M609 1001 2 8016 1 MT609 LWTR M111 6000 1 1001 1 7014 1 MT111 POLY BENZ GRPH sdef pos=10 6 0 nps 1000000 PRINT

Most photon reactions are photoatomic and involve interactions with the atom’s electrons. Consequently, photoatomic data are elemental rather than isotopic. The MCNP code will convert

m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 m313 92235 -0.02759 92238 -0.85391 8016 -0.11850 M609 1001 2 8016 1 M111 6000 1 1001 1 7014 1

to

m201 40000 -0.505239 40000 -0.110180 40000 -0.168413 40000 -0.170672 40000 -0.027496 m313 92000 -0.02759 92000 -0.85391 8000 -0.11850 M609 1000 2 8000 1 M111 6000 1 1000 1 7000 1

Note that the Zr in material M201 will still be sampled for five different mass fractions of the same Zr element 40,000; uranium in material M313 will be sampled for two different mass fractions for 92,000.

Figure 2.19 is an MCNP plot of hydrogen, oxygen, zirconium, and uranium total photon cross sections. MCNP is run in the cross-section plotting mode:

MCNP6 I=Ex12_2a name=junk. IXZ

and the MCNP plot command is (see section 5.3.3 of the MCNP manual):

ylim .01 1e7 par p mt -5 xs 1000.04p cop xs 8000.04p & cop xs 40000.04p cop xs 92000.04p

Fig. 2.19
figure 19

Total photon cross section of H, O, Zr, and U

Photons are considered absorbed and no longer transported below 1 keV. Photon cross sections go up to 100 MeV. The five photoatomic reactions in MCNP are (1) incoherent (Klein–Nishina Compton scatter plus form factors); (2) Coherent (Thomson scatter plus form factors); (3) photoelectric (with k- and l-shell fluorescence); (4) pair production, in which the positron immediately decays into two photons; and (5) bremsstrahlung, generated by electrons. The first four reactions are treated rigorously, but bremsstrahlung is approximate unless electrons are also transported. Electron transport is a thousand times slower than photon transport and is approximate. Thus, the default photon treatment is the thick-target bremsstrahlung (TTB) approximation. TTB produces all of the photons that would be produced by full electron transport at the point of interaction and in the direction of the photoatomic collision without transporting the electrons. In a thick geometric region, this is a good approximation: after many electron interactions, the electrons remain in the same region, move insignificantly, and have essentially random directions. In a thin geometric region where the electrons lose only a few percent of their energy while slowing down, TTB is a poor approximation and a full electron transport is required.

Figure 2.20 is an MCNP plot of the basic uranium photon cross sections. MCNP is run in the cross-section plotting mode:

MCNP6 I=Ex12_2a name=junk. IXZ

and the MCNP plot command is

ylim .01 1e7 par p xs 92000.04p mt -5 cop mt -1 cop mt -2 cop mt -3 cop mt -4

Fig. 2.20
figure 20

Photon cross sections for uranium: total (MT = −5), incoherent (MT = −1), coherent (MT = −2), photoelectric (MT = −3), and pair production (MT = −4)

Bremsstrahlung is not included in the photon libraries because it is generated by electron models.

Photonuclear reactions are those in which photons interact with the nucleus. Because they are generally low probability, the MCNP code ignores them unless specifically requested with a PHYS:P card. Photonuclear reactions produce neutrons, photons, and other particles, and a limited number of isotopic photonuclear cross sections are available with MCNP. The ZAID number ends in a u, as in 6012.70u. Photonuclear reactions should be modeled only if secondary photonuclear particles are of particular interest.

2.2.6 Electron-Stopping Powers for Coupled Photon and Electron Problems

Example 2.11 becomes a coupled photon-electron problem by replacing the MODE card and M0 data cards:

MODE P E M0 PLIB=04p ELIB=03e ESTEP=10 COND GAS

In the above example, photons use the .04p data, electrons use the .03e data, and the material is treated electrically as a gaseous conducting material. ESTEP will be described soon.

Coupled electron-photon transport is possible below 10 keV with the single-event electron transport method by specifying data from the EPRDATA14 photon library (.14p) and setting the 15th entry on the PHYS:E card. The single-event electron-photon transport models every interaction of the electron and photon random walk and is very slow and approximate. Generally, the condensed history method is used, and photons and electrons below 1 keV are simply absorbed.

Electron transport in MCNP for energies E > 1 keV utilizes the condensed history method.

The condensed history method models all interactions of an electron over an energy substep. An electron passing through matter will interact with each atom along its trajectory via ionization, angular scattering, and production of secondaries. Modeling this interaction exactly by the single-event model is time prohibitive. The condensed history algorithm attempts to average the effect of all of these interactions into aggregate quantities, which can be approximated over a series of major steps and several substeps. The major steps are controlled by EFAC, the 14th entry on the PHYS:E card. Each step is at an energy EFAC = 0.917 times the energy of the previous (higher) step, which provides eight energy steps for a factor of two energy loss. Setting EFAC = 0.99 makes each step at an energy 0.99 of the previous step, which provides a much finer and more accurate energy grid but increases computational time by a factor of about ten. The substeps are controlled by ESTEP on the materials card and other factors such as the particular material and energy. The default is ESTEP = 3, or 3 substeps per major or full step. The step size can be different for each material.

Here, ESTEP = 10 specifies 10 substeps for each full step. M0 applies this to all materials for which ESTEP is not otherwise specified. The other parameters specify that photons use the .04p data, electrons use the .03e data, and the material is treated electrically as a gaseous conducting material.

The highest energy, EMAX, is the 1st entry on the PHYS:E card (100 MeV by default). An electron starting at some energy E will use the electron data (stopping powers, slowing-down parameters, etc.) at its energy range and travel the distance needed to lose a substep’s worth of energy. When it loses a full step of energy, or 1.0–0.917 = 8.3% of its energy, it will then use the data of that step. For example, by default, an electron at 1 MeV will use the 1 MeV data until its energy falls below 917 keV. It will then use these data until its energy falls below 814 keV, etc. In each of these ranges, it will have three substeps that determine how far the electron travels and its deflection direction. At each step, bremsstrahlung photons, photon-induced electrons, electron induced X-rays, and knock-on electrons are produced and these production mechanisms may be biased with the PHYS:E card.

Electron transport is approximately 1000 times slower than photon transport. Changing ESTEP from 3 to 20 slows down electron transport by another factor of 10. Changing EFAC to the maximum resolution of EFAC = 0.99 on the PHYS:E card slows down electron transport by another factor of 10. But the default values of ESTEP = 3 and EFAC = 0.917 may result in too coarse an energy grid structure and wrong answers. If electron transport is unimportant, then the thick-target bremsstrahlung approximation should be sufficient. If the calculation is sensitive to the electron transport, then we recommend checking if higher resolution (increasing EFAC and ESTEP) affects results. Turning off knock-on electrons (RNOK=0, 8th entry) will often speed calculations by a factor of seven without adverse effects. The effect of knock-on electrons is problem dependent. It is recommended to run first with knock-on electrons turned on (default) to include knock-on effects; then turn off knock-ons to see if tallies are unaffected and the calculation runs faster.

Positron and electron physics are identical except when they stop (fall below energy cutoff). Electrons deposit energy, whereas positrons generate annihilation photons. At high energies, positrons behave like electrons. At low energies (<1 MeV), stopping powers, bremsstrahlung, knock-ons, and annihilations are increasingly poor. Positron sources are allowed:

SDEF par = -e SDEF par = f

Positrons may be tallied separately:

FTn ELC 3

The electron-stopping powers and other parameters are printed for each energy step and each material in the problem. We recommend printing this information in the output file only when changing materials to avoid excessively long output files; otherwise, eliminate this printout by using

Print -85 -86

Figure 2.21 is an MCNP plot of basic three electron materials. The MCNP code is run in the cross-section plotting mode:

MCNP6 I=Ex12_2c name=junk. IXZ

Fig. 2.21
figure 21

Total electron-stopping (MT = 3) power for material M609 (water), M201 (Zr), and M313 (UO2)

And the MCNP plot command is

ylim 1 200 par e xs m609 cop xs m201 cop xs m313

The total, collisional, and radiative stopping powers are plotted in Fig. 2.22 from the same input file with the following command:

par e xs m313 ylim .5 30 mt 3 cop mt 1 cop mt 2

Fig. 2.22
figure 22

Electron-stopping powers for material M313, UO2: total (MT = 3), collisional (MT = 1) and radiative (MT = 2)

The stopping powers are for materials but not nuclides. Thus, in the cross section, plot command XS = M313 works, but XS = 92000.03e does not. Electron-stopping powers also depend on the material density. Thus, if a material is at more than one density, the stopping powers for the first density will be correct but wrong for all others. Each material of a different density requires a new material number even though the nuclides are the same.

2.2.7 Data and Models for Ions and Charged Particles

The MCNP code models 36 light ions and subatomic particles and also transports 2200+ heavy ions. See “MCNP6 Particles” Table 2-2 in the MCNP code manual. Though there are a few data libraries for protons, deuterons, tritons, alphas, and helions, most nuclides and all the other particle types (deuteron o, triton r, alpha a, helion s) and heavy ions use physics models and continuous slowing-down theory (like electrons). A number of physics models have been developed internationally and these can be selected along with key parameters on the LCA, LCB, LEA, LEB, and LCC input cards. The principal MCNP models are CEM, LAQGSM (default) and Bertini, Isabel, and INCL, which have their own specific data files (such as BERTIN, FALPHA, and CINDERGL) provided in the MCNP code package. See the MCNP manual [2] section 3.3.3.7.1.

The following example shows how the MODE, M0, and M111 cards might be changed for a charged particle problem:

MODE N P E H D T S A | / Z m0 nlib=.70c plib=.04c elib=.03e alib=.00a slib=.00s MT111 POLY BENZ GRPH M111 6000 1 1001 1 7014 1 MX111:P 6012.24U 0 7014.70U $ Photonuclear MX111:H 6012.70H 1001.71H 7014.24H $ Proton MX111:D 6012.00o 1001.00o 7014.00o $ Deuteron MX111:T 6012.00R 1001.00R 7014.00R $ Triton MX111:S 6012.00S MODEL 7014.00S $ Helion MX111:A 6012.00A MODEL 8016.00A $ Alpha

The MODE card will enable transport of neutrons (N), photons (P), electrons (E), protons (H for hadron), deuterons (d), tritons (T), helions (S for 3He), and alphas (A). Good practice is to add negative muons (|), positive pions (/), and neutral pions (Z) to the MODE card when energies in the problem reach the level of rest mass for these particles (~100–150 MeV). Add positive kaons (K) for energies >~500 MeV. Electron neutrinos (U) and muon neutrinos (V) often are added.

The M0 card is the default for all materials cards. NLIB, PLIB, ELIB, ALIB, and SLIB are for neutron, photoatomic, electron, alpha, and helion data.

The MX cards will enable nuclide substitution. No substitutions are allowed for photoatomic (P) and electron (E) data because these data are elemental rather than isotopic. Carbon on the M111 card is specified as 6000 for all modern carbon cross-section libraries, but carbon is available only as 6012 for photonuclear, proton, and other light nuclei libraries, thus requiring the MX111 cards. Also, isotopes rather than elements should be specified for all model physics. MX111:P is for photonuclear data and specifies that no cross section is used for photonuclear hydrogen because hydrogen has no photonuclear reactions. The library types for deuterons, tritons, helions, and alphas are O, R, S, and A. MX111:S and MX111:A are for helions and alphas and specify model physics rather than data libraries, which are used for hydrogen. For neutrons, photoatomic, photonuclear, and electrons, the data libraries are best in the energy range for which they are provided. For protons, the data libraries—if they are available—are sometimes better for secondary particle production but not transport, where models are sometimes better. For the remaining charged particles and ions, the models are generally better; data libraries, such as TENDL (https://tendl.web.psi.ch/tendl_2019/tendl2019.html) must be used with caution.

2.2.8 Additional Data Diagnostics and Recommendations

As illustrated in the previous sections, MCNP data can be plotted. See Sect. 2.5 for details.

The MCNP code has a “first collision” feature for studying the physics in isolation. It is activated by the 8th entry on the LCA card (MCNP manual section 3.3.3.7.2):

LCA 7J -2.

Source particles do not transport and instead initiate an immediate collision. Products of the reaction are then transported out as if in a void. The option works for library and model physics. With this option, F1 or other tallies can be used to tally angle dependence, energy dependence, multiplicities, and secondary particle production of any particle type that collides with any nucleus.

The MCNP code also has a GENXS option. With a beam of particles or nuclei impinging upon a thin target and the GENXS option, the MCNP code can easily calculate double-differential and angle/energy integrated spectra, multiplicities, yields of emitted particles, fission cross sections, and fission yields. The GENXS option is like using the event-generator of the MCNP6 code as a stand-alone code—without any real transport through the matter. A small and simple second auxiliary input file is required and is described along with the GENXS option in the MCNP manual (section 3.3.3.9).

The following recommendations avoid common user errors:

  • Mn material card charged particle steps:

    • ESTEP (electrons) and hstep (all other charged particles) have a default of 3 substeps, which may be too small; try something larger like 10 or 20 to see if results—particularly secondary energy distributions—change, although this will run slower;

    • EFAC on the PHYS card for each charged particle type controls the slowing-down energy step size and the default, EFAC = 0.917, may be too small; try something larger like 0.99 to see if results change, although this will run slower;

  • Upper energy bounds: The PHYS:x for all particle types has Emax as the first entry, which must be set above maximum energy of the problem to prevent wrong answers. It is recommended to set it just above to avoid inefficient effects of loading unneeded data;

  • MODE n h d t s a | / z k --- add d t s a as needed; | / z for E > 100 MeV; k for E > 500 MeV.

  • Lower energy cutoff. On the CUT:x card, set the 2nd entry, Emin, to 0.001 to go below the high default E lower energy cutoff. The default cutoff is the A number of the projectile. Thus, alpha particles transported below 4 MeV will be absorbed unless Emin is lowered.

  • Consider analog capture, 3rd and 4th entry = 0, on CUT:x, to avoid 0 weight window complications.

2.3 Sources

2.3.1 SDEF Fixed Sources

The input in Example 2.12 is similar to Example 2.6 with additional materials specifications. It is illustrated in Fig. 2.14 through Fig. 2.16.

Example 2.12 Fully Specified FILL Square and Hex Lattices. Square and Hex Lattices

1 313 -10.44 -1 u=1 imp:n=1 $ UO2 2 201 -6.55 -2 1 u=1 imp:n=1 $ Zr clad 3 609 -1.00 2 u=1 imp:n=1 $ H2O 21 111 -4.00 -31 lat=1 u=11 imp:n=1 $ Square lattice fill = -3:3 0:0 -3:3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 1 1 1 1 1 1 22 111 -4.00 -32 lat=2 u=12 imp:n=1 $ Hexagonal lattice fill = -5:5 -4:4 0:0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 12 1 1 1 1 1 1 0 0 0 0 14 0 -4 fill 11 imp:n=1 $ Square lattice box 11 like 14 but TRCL=(20 0 0) fill 12 $ Hexagonal lattice box 15 0 #14 #11 -20 imp:n=1 $ outside boxes 16 0 20 imp:n=0 $ outside world 1 rcc 0 0 0 0 12 0 .5 $ Right circular cylinder 2 rcc 0 0 0 0 12 0 .6 $ Right circular cylinder 4 rpp -7 7 0 12 -7 7 $ Right parallelepiped 14x12x14 31 rpp -1 1 0 12 -1 1 $ Right parallelepiped 2x12x2 32 rhp 0 0 0 0 12 0 1 0 0 $ Right hex 12x1 20 sph 0 0 0 50 $ Sphere containing lattice boxes MODE N P m0 nlib=.70c plib=.04c elib=.03e m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 nlib=80c m313 92235 -0.02759 92238 -0.85391 8016 -0.11850 M609 1001 2 8016 1 MT609 LWTR.10t M111 6000 1 1001 1 7014 1 MT111 POLY.10t BENZ.10t GRPH.10t sdef pos=10 6 0 nps 1000000 PRINT

The MCNP standard source definition, SDEF, may have no entries, fixed entries, distributions, and dependent distributions.

The SDEF in Example 2.12 places the source point at position POS=10 6 0 instead of the default POS=0 0 0 cm

sdef pos=10 6 0

Example 2.12 SDEF is the same as

SDEF X=10 Y=6 Z=0 ERG=14 TME=0 PAR=N CEL=15 EFF=.01 WGT=1.0

Dimensions (POS, X, Y, Z, RAD, EXT) are in cm, energies (ERG) in MeV, and times (TME) in shakes, 10−8 s.

Particle types are listed in Table 2-2 of the MCNP manual. In addition, heavy ions may be source particles, as in PAR = 6012, which requires heavy ions to be specified on the MODE card as MODE #. Antiparticles are specified as negative: PAR = −E is a positron source. Particle types may also be PAR = SF spontaneous fissions, PAR = −SF spontaneous fission neutrons, or PAR = SP spontaneous photons.

The default cell is wherever the source is located. If the cell is specified, CEL=15, then any part of the source not in that cell is rejected. Consequently, CEL specifies a rejection criterion rather than a cell location. If more than 1% of the source locations are rejected, then the calculation is stopped because the default source rejection efficiency is EFF=0.01. In some calculations, it is desirable to set the source rejection efficiency lower, such as EFF=0.0001.

The particle weight, WGT, is the physical number of particles modeled in the source. The default is 1 source particle so that all results may be easily normalized to the true source strength. All results are normalized to the source weight, WGT.

Additional SDEF scalar variables are SUR, NRM, ARA, TR, and CCC (see table 3-62 of the MCNP manual).

Particles starting exactly on a surface must have that surface name specified with SUR. If the particles are directed inward or outward relative to the source normal, then NRM is specified as NRM=1 in the direction of the positive source normal and NRM=−1 in the direction of the negative source normal.

If a source is mono-directional, then the surface area of the source, ARA, must be specified if there are point detector tallies, F5, or DXTRAN variance reduction.

The location of source particles may also be translated and rotated, where TR=n refers to the TRn card translation and rotation.

Source particles may also be rejected if they are outside a region defined by a cookie cutter cell (CCC), which is not necessarily a cell in the problem.

2.3.2 SDEF Source Distributions

Additional SDEF keywords are vectors VEC and AXS, which are associated with distributions RAD and EXT.

In the following example, only the region inside the hexagonal box is sampled after starting particles uniformly in volume in an overlaying sphere. Replace the SDEF in Example 2.12 with

SDEF POS=10 6 0 RAD=D11 CCC=11 EFF=.0001 SI11 0 20 SP11 -21 2

The combination of POS and RAD samples sources particle location in a sphere. Radius is now a distribution, RAD=D11. The source information card, SI11, distributes the radius between 0 and 20 cm. The source probability, SP11, invokes a power law (distribution type −21), where the radius is raised to the second power in order to sample uniformly in volume of the sphere. The second power is needed because the differential volume of a sphere is

$$ dV=4\pi {r}^2 dr $$

Note that CCC rejects any point sampled in the sphere volume that is not in the cookie cutter cell, CCC=11, which is the hexagonal lattice box. The efficiency of sampling is 0.0696. Whereas some source particles do not find a source location within the box after 100 tries, MCNP stops with the error message

expire parameter is the sampling efficiency in source cell 11 is too low

The problem immediately terminates with a “bad trouble” error; therefore, the efficiency must be set lower. In this case, EFF=0.0001 works.

To sample a cylindrical volume, the combination POS, RAD, EXT is used (replacing SDEF in Example 2.12):

SDEF POS=0 6 0 RAD=D12 EXT=D13 AXS=0 1 0 SI12 0 10 SP12 -21 1 SI13 H -6 6 SP13 0 1

The radius is sampled (SI12) between 0 and 10 cm as a power law raised to the 1st power (SP12 −21 1) because the differential volume of a cylinder with respect to r is

$$ dV=2\pi rh\kern0.5em dr $$

The extent is sampled in the Y-direction (AXS=0 1 0) as a histogram (SI13 H −6 6) from −6 < Y < 6 cm relative to the position (POS=0 6 0). The distribution is linear with probability (SP13 0 1) 0 for Y < −6, 1 between −6 < Y < 6. The H designator is one of the several possibilities:

  • H—histogram bin upper boundaries (default)

  • L—discrete values follow

  • A—points where probability density distribution is defined

  • S—distribution numbers follow

The SP options are as follows:

  • D—bin probabilities for an H or L distribution (default)

  • C—cumulative bin probabilities for an H or L distribution

  • V—probability is proportional to cell volume

  • W—intensities for a mix of particle sources

Particles may be started only in cell 1 of the square lattice (replacing SDEF in Example 2.12) by making CEL a distribution; however, the full path through the universe levels must be specified when the cell is in a universe. Particles are started in all 48 locations of cell 1:

SDEF POS=0 6 0 RAD=D12 EXT=D13 AXS=0 1 0 CEL=D14 EFF=1e-5 SI12 0 10 SP12 -21 1 SI13 H -6 6 SP13 D 0 1 SI14 L (1<21<14) SP14 1

Note that EFF had to be set because the efficiency is 0.0025, less than the default 0.01.

Particles may also be started in all 66 locations of cell 1 of the hexagonal lattice by specifying which filling universes of cell 22 to use. Because the universe (u = 1) is specified, the probabilities of the individual cell 1s must be specified (SP14) and can therefore be different if desired:

SDEF POS=0 6 0 RAD=D12 EXT=D13 AXS=0 1 0 CEL=D14 SI12 0 .5 SP12 -21 1 SI13 -6 6 SP13 0 1 SI14 L (1<22[u=1]<11) SP14 1 65R

In both this and the previous case, the radius may simply be the cell 1 radius (SI12 0 .5) and no EFF is needed because the source efficiency is 100% in each cell 1.

The most common specification of cells in a lattice specifies the path to each individual repetition of the cell. In the following, particles are started in all 66 locations of cell 1 of the hexagonal lattice but with different probabilities depending upon whether the cell is cut off or not (Fig. 2.14):

SDEF POS=0 6 0 RAD=D12 EXT=D13 AXS=0 1 0 CEL=D14 SI12 0 .5 SP12 -21 1 SI13 H -6 6 SP13 0 1 SI14 L (1<22[-5 3 0]<11) (1<22[-4 1:4 0]<11) (1<22[ -3 -1:4 0]<11) (1<22[-2 -3:4 0]<11) (1<22[-1:1 -4:4 0]<11) (1<22[ 2 -4:3 0]<11) (1<22[ 3 -4:1 0]<11) (1<22[ 4 -4:-1 0]<11) (1<22[ 5 -4:-3 0]<11) SP14 .5 .5 1 1 .5 .5 1 3R .5 .5 1 5R .5 .5 1 6R .5 .5 1 6R .5 .5 1 6R .5 .5 1 5R .5 .5 1 3R .5 .5 1 1 .5 .5 .5

The square brackets indicate the lattice index numbers. The “:” notation on the SI entries expands the 9 entries to the 66 lattice elements. The SP entries are .5 for elements that are cut off by the filling cell in the universe level above.

2.3.3 SDEF Dependent Distributions: DS

The SDEF definition and source distributions from Example 2.19 in Sect. 2.4.3 and shown in Example 2.13 illustrate dependent distributions. The MCNP code allows one level of dependency, but this example illustrates how to model three levels of dependency: energy is dependent upon time which is dependent upon location dependent upon the particle type.

Example 2.13 Dependent Source Description Taken from Example 2.19

SDEF PAR=D42 CEL=FPAR=D51 POS=FPAR=D53 RAD=FPAR=D54 EXT=FPAR=D55 AXS=FPAR=D56 TME=FPAR=D57 ERG=FPAR=D58 SI42 L SF N P P P P $ PAR SP42 .3 .3 .1 .1 .1 .1 $ PAR DS51 S 46 47 48 48 48 48 $ CEL DS53 L 0 6 0 0 6 0 10 3 0 10 3 0 10 7 0 10 7 0 $ POS DS54 S 43 43 0 0 0 0 $ RAD DS55 S 44 44 0 0 0 0 $ EXT DS56 S 45 45 0 0 0 0 $ AXS DS57 L 0 1e8 2e8 3e8 4e8 5e8 $ TME DS58 S 0 92 93 94 95 96 $ ERG SC46 Spontaneous fission cells in XYZ lattice SI46 L (1<21[-2:3 0 3]<14) (1<21[-3:3 0:0 -3:2]<14) SP46 1 5R 1 41R SC47 Neutron source cells in HEX lattice SI47 L (1<22[-5 3 0]<11) (1<22[-4 1:4 0]<11) (1<22[ -3 -1:4 0]<11) (1<22[-2 -3:4 0]<11) (1<22[-1:1 -4:4 0]<11) (1<22[ 2 -4:3 0]<11) (1<22[ 3 -4:1 0]<11) (1<22[ 4 -4:-1 0]<11) (1<22[ 5 -4:-3 0]<11) SP47 .5 .5 1 1 .5 .5 1 3R .5 .5 1 5R .5 .5 1 6R .5 .5 1 6R .5 .5 1 6R .5 .5 1 5R .5 .5 1 3R .5 .5 1 1 .5 .5 .5 SC48 Point photon source cell in water between lattices SI48 L 15 SP48 1 c particle-dependent radius, extent, axis SI43 0 .5 $ Radii SP43 -21 1 $ Radii SI44 -6 6 $ Extent SP44 -21 0 $ Extent SI45 L 0 1 0 $ Axis SP45 1 $ Axis c energies dependent upon time, position, and particle type SP92 -3 1.175 1.0401 $ 252Cf Frohner Watt parameters SI93 50 59 $ Photon, XYZ water, TME 2e8, energies SP93 0 1 SI94 60 69 $ Photon, XYZ water, TME 3e8, energies SP94 0 1 SI95 L 70 $ Photon, XYZ water, TME 4e8, energies SP95 1 SI96 L 80 $ Photon, XYZ water, TME 5e8, energies SP96 1

The SDEF definition has particle type (PAR) as an independent distribution and then dependent distributions for cell (CEL), position (POS), radius (RAD), extent (EXT), axis (AXS), time (TME), and energy (ERG).

The particle types are spontaneous fission (SF), neutrons (N), and then photons (P). The photons are repeated four times (SI42) to enable more levels of dependencies: there is an energy distribution for each time, times are sampled from two locations, and locations are sampled from three particle types.

The DS cards point to the values or distributions for each of the six particle types.

DS51 provides the distributions of cells for each particle type. Spontaneous fission occurs in the XYZ lattice cells, distribution 46. Note that each spontaneous fission will produce between zero and ten spontaneous fission neutrons. The neutron source is in the HEX lattice cells, distribution 47. The photon source is in cell 15, distribution 48. Note that each of these distributions has a title description on the source comment (SC) card.

DS53 provides the source position for each particle type. Note that the position in the XYZ lattice and the HEX lattice is the same, namely the [0 0 0] lattice element is a POS=0 6 0 in the coordinate system of the filling square lattice boxes, which are 20 cm apart. The photons start between the lattice boxes at POS=10 3 0 or POS=10 7 0.

DS54, DS55, and DS56 point to the distributions for radius, extent, and axis for the fuel pin cylinders. Whereas the photon source is at points, there is no radius, extent, or axis description, so the photon distributions are set to distribution zero or no distribution.

DS57 sets the time of spontaneous fission (0 s), neutron source (1 s), and four photon times, 2, 3 4 and 5 s at POS=10 3 0 and 4 and 5 s at POS=10 7 0.

DS58 specifies the energy distributions. None is specified for spontaneous fission because the appropriate Watt spectrum for each spontaneous fission is automatically sampled by the MCNP code. Neutrons are started in the HEX lattice with a 252Cf Frohner Watt spectrum [5], SP92. Note that no SI92 is used for special Watt distribution (−3) on SP92. The photon energies are

  • 50 < E < 59 MeV for TME=2 s and POS=10 3 0;

  • 60 < E < 69 MeV for TME=3 s and POS=10 3 0;

  • E = 70 MeV for TME=4 s and POS=10 7 0; and

  • E = 80 MeV for TME=5 s and POS=10 7 0

Less common options for the dependent distribution DS card are H for dependent continuous histogram and T and Q for dependent values of independent values. For example,

ERG = D1 DIR = FERG = D3 SI1 .001 .01 .1 .4 1.3 2.6 5 14 DS3 Q .3 4 1.8 5 14 6

The DS3 card with the Q option specifies

From 0–0.3 MeV

Get DIR from distribution SI4

From 0.3–1.8 MeV

Get DIR from distribution SI5

From 1.8–14 MeV

Get DIR from distribution SI6

2.3.4 Criticality Sources

MCNP criticality calculations are used in safeguards to ensure criticality safety. Nuclear criticality is the determination of keff, which is the number of neutrons in a system produced by an initial number of neutrons. If keff <1, the system is subcritical, and more neutrons leak out or are captured than those produced by fission. If keff = 1, then the system is critical, and the number of neutrons produced by fission precisely balances the number lost to leakage and other physical processes [6]. See Example 2.14.

Example 2.14 Four 20 wt% Enriched Uranium Cans (Radius 8 cm, Height 40 cm) in Air (Radius 25 cm, Height 42 cm) Inside an Aluminum Casing (Radius 40 cm, Height 44 cm) with Criticality (KCODE) Source

Four uranium cans in air and aluminum 11 101 -18.0 -51 imp:n=1 $ uranium can 12 like 11 but trcl=(20 0 0) imp:n=1 $ uranium can 13 like 11 but trcl=( 0 20 0) imp:n=1 $ uranium can 14 like 11 but trcl=(20 20 0) imp:n=1 $ uranium can 21 102 -1.0 -52 #11 #12 #13 #14 imp:n=1 $ air tank 22 103 -2.7 -53 52 imp:n=1 $ Al liner 23 0 53 imp:n=0 $ outside 51 RCC 0 0 -20 0 0 40 8 $ fuel cylinder 52 RCC 10 10 -21 0 0 42 25 $ air cylinder 53 RCC 10 10 -22 0 0 44 40 $ aluminum liner m101 92238 -.8 92235 -.2 nlib=.80c m102 7014.80c .8 8016.80c .2 m103 13027.80c 1 print kcode 1e4 1.0 30 130 sdef axs=0 0 1 rad=d31 ext=d32 pos=d33 si31 0 8 sp31 -21 1 si32 -20 20 sp32 0 1 si33 L 0 0 0 20 0 0 0 20 0 20 20 0 sp33 1 3r FC1 Neutrons crossing surfaces F1:n (51.1 51.2 51.3) (52.1 52.2 52.3) (53.1 53.2 53.3) C1 0 1 T FQ1 F C FC4 Fissions and neutron production in uranium cells F4:n 11 12 13 14 T FM4 (-1 101 -6) (-1 101 -6 -7) FQ4 F M SD4 1 4r

Figures 2.23 and 2.24 show the axial and cut-away views of the criticality problem of Example 2.14.

Fig. 2.23
figure 23

X – Y view of Example 2.14, four uranium cans in air

Fig. 2.24
figure 24

X – Z view of Example 2.14, four uranium cans in air

Example 2.14 is a contrived model of four uranium cans in air surrounded by aluminum. The source is specified by a KCODE card

KCODE 1000 1.0 30 130

The KCODE default values are shown above.

The first KCODE entry is the number of histories to run in a generation, or “batch” or “cycle.” This number should be as large as possible in order to sample all geometry, particularly fissionable regions. Because each generation of fission neutrons is based on the previous generation, the statistical sampling of keff is correlated and thus violates the central limit theorem of statistics. Consequently, the estimated standard deviation of keff is underpredicted, which can yield incorrect results. Further, nuclear systems must be “delayed critical,” which means that the part of keff from delayed neutrons—less than 1% of the fission neutrons produced—requires that the relative error allowed for keff is << 0.01. Further, variance reduction cannot be used to better converge results because keff is based on the balance of the whole system and variance reduction works by focusing on important parts of phase space at the expense of others.

The second KCODE entry is the guess of keff for the first cycle only. It does not need to be predicted closely and can be a factor of two off. If it is too low, then MCNP will issue the message

new source has overrun old source. will start over with larger k.

If the guess is too high, then MCNP will waste a lot of storage space. So the default of 1.0 is fine.

The third KCODE entry is the number of settle cycles to run, that is, the number of cycles to be skipped before beginning tally accumulation. Whereas the initial fission distribution must be specified by KSRC, SDEF, or SRCTP, a number of cycles are needed to better distribute the fission source before beginning the calculation of keff and any tallies. The default of 30 should be the minimum. If keff or tallies of interest do not converge, then more may be needed.

The fourth KCODE entry is the total number of cycles to run altogether. The total number of histories in the problem is approximately the product of the batch size and the total number of cycles. With these default values, the first 30,000 histories are discarded—after hopefully converging the fission distribution—and then keff and the tallies are calculated from the next 100,000 histories, for a total of about 130,000 histories. Note that the batch size varies slightly from cycle to cycle to correctly adjust for particle weighting; consequently, the final number of histories will not be exactly the product of the requested batch size and requested number of cycles.

The initial fission distribution for the first cycle only is determined by KSRC or SDEF cards or a SRCTP file. Example 2.14 uses

sdef axs=0 0 1 rad=d31 ext=d32 pos=d33 si31 0 8 sp31 -21 1 si32 -20 20 sp32 0 1 si33 L 0 0 0 20 0 0 0 20 0 20 20 0 sp33 1 3r

In this case, the initial fission distribution is uniform in volume in each of the four cans of uranium. Alternatively, the initial fission distribution could be specified by points by replacing the above with one fission event in the center of each uranium can

KSRC 0 0 0 20 0 0 0 20 0 20 20 0

which is the same as

sdef pos=d33 si33 L 0 0 0 20 0 0 0 20 0 20 20 0 sp33 1 3r

Alternatively, and best, is to use the SRCTP file from one run as the initial source for the next if all the fission cells have the same geometry. The SRCTP has the source locations and energies from a previous calculation which is better than an approximate guess. In this case, the name of the source file is specified on the MCNP execution line and neither KSRC nor SDEF is present.

Example 2.14 also has two tallies specified after the source description. These are described in more detail in Sect. 2.4, Output and Tallies, and are the number of neutrons crossing surfaces 51, 52, and 53 and the number of fissions and number of fission neutrons in cells 11, 12, 13, and 14.

Because keff must be known very precisely, with an error << 1%, and because the relative errors of tallies and standard deviations of keff are underestimated by violating the central limit theorem of statistics, the MCNP code calculates thousands of values of keff in each run. The MCNP code calculates keff in seven different ways using three uncorrelated estimators (collision, absorption, track length) and four additional combinations (collision/absorption, absorption/track length, track length/collision, collision/absorption/track length). These estimators and combinations are for every cycle individually, as simple averages and as covariance-weighted averages. Fortunately, the “the final estimated combined collision/absorption/track-length keff” is printed in a box in the output file. How well this final keff is converged is indicated by how well the keff values agree and deviate from the final.

We highly recommend tallies for KCODE calculations to provide confidence in the overall calculation of keff. Tally 4 in Example 2.14 uses the FM card (see Sect. 2.4.3.6) to multiply the flux in each fissionable can by ρσf and ρνσf to get the number of fissions and number of fission neutrons produced. The tally 4 output is

cell number of fissions number of fission neutrons 11 7.77321E-02 0.0018 1.99838E-01 0.0018 12 7.68256E-02 0.0018 1.97500E-01 0.0018 13 7.74996E-02 0.0018 1.99243E-01 0.0018 14 7.79120E-02 0.0018 2.00343E-01 0.0018 total 3.09969E-01 0.0007 7.96924E-01 0.0007

The total number of fission neutrons agrees well with “the final estimated combined collision/absorption/track-length keff = 0.79691 with an estimated standard deviation of 0.00056”; however, the tallies between the different fissionable cells should be the same due to symmetry. Instead, they are 1.44% apart with 0.0018 relative errors, which implies that the fission distribution is not converged. We would recommend running ten times as many histories per cycle, KCODE 1e5, to seek better convergence of the fission distribution. But since keff << 1.0, further accuracy may not be needed because the system is subcritical below any reasonable margin of error.

Section 2.5 describes how to plot tally output. Figure 2.34 in that section is an MCNP tally plot showing the convergence of this Example 2.14 criticality calculation. The combined C-A-T—collision, absorption, and track length—tallies are plotted together as a function of cycle number and the plot commands are provided.

2.3.5 Surface Source Write and Read (SSW, SSR)

The surface source enables particles that cross various surfaces in a problem to be recorded (SSW—surface source write run) and then rerun in a subsequent problem (SSR—surface source read run), omitting all surfaces leading up to it. Thus, the subsequent SSR run can be used for multiple parameter studies without having to repeat the tracking through the initial part of the problem.

A frequent use in safeguards is to model a bremsstrahlung source with full and time-consuming electron transport in the first SSW run, writing the resulting photons to a surface source file. In the subsequent SSR run, only the photons from the surface source file are followed without the expensive initial electron run. See Example 2.15.

Example 2.15 SSW Input File to Write a Surface Source. The Geometry is the Same as That for the Criticality Example 2.14: Four 20 wt% Enriched Uranium Cans (Radius 8 cm, Height 40 cm) in Air (Radius 25 cm, Height 42 cm) Inside an Aluminum Casing (Radius 40 cm, Height 44 cm) with a Fixed SDEF Source

Four uranium cans in air and aluminum 11 101 -18.0 -51 imp:n=1 $ uranium can 12 like 11 but trcl=(20 0 0) imp:n=1 $ uranium can 13 like 11 but trcl=( 0 20 0) imp:n=1 $ uranium can 14 like 11 but trcl=(20 20 0) imp:n=1 $ uranium can 21 102 -1.0 -52 #11 #12 #13 #14 imp:n=1 $ air tank 22 103 -2.7 -53 52 imp:n=1 $ Al liner 23 0 53 imp:n=0 $ outside 51 RCC 0 0 -20 0 0 40 8 $ fuel cylinder 52 RCC 10 10 -21 0 0 42 25 $ air cylinder 53 RCC 10 10 -22 0 0 44 40 $ aluminum liner m101 92238 -.8 92235 -.2 nlib=.80c m102 7014.80c .8 8016.80c .2 m103 13027.80c 1 print nps 100000 sdef axs=0 0 1 rad=d31 ext=d32 pos=d33 erg=d34 si31 0 8 sp31 -21 1 si32 -20 20 sp32 0 1 si33 L 0 0 0 20 0 0 0 20 0 20 20 0 sp33 1 3r sp34 -3 .965 2.29 FC1:n Neutrons crossing surfaces F1:n (51.1 51.2 51.3) (52.1 52.2 52.3) (53.1 53.2 53.3) C1 0 1 T FQ1 F C FC4 Fissions and neutron production in uranium cells F4:n 11 12 13 14 T FM4 (-1 101 -6) (-1 101 -6 -7) FQ4 F M SD4 1 4r SSW 52.1 52.2 52.3

Example 2.15 is the same geometry as the criticality Example 2.14 illustrated in Figs. 2.23 and 2.24. The geometry of Example 1.3-3 is used here to avoid introducing an entirely new geometry. But the criticality source (KCODE) is replaced with surface source write (SSW) and surface source read (SSR) making Example 2.15 an entirely different problem for illustrative purposes. The criticality source

kcode 1e4 1.0 30 130 sdef axs=0 0 1 rad=d31 ext=d32 pos=d33

is replaced with

nps 100000 sdef axs=0 0 1 rad=d31 ext=d32 pos=d33 erg=d34 sp34 -3 .965 2.29

where the energy distribution, SP34, is a Watt spectrum (special function −3) with the Watt parameters of the initial default KCODE cycle fission distribution (a = 0.965 MeV, b = 2.29 MeV−1).

In addition, the surface source write card is added

SSW 52.1 52.2 52.3

which writes all particles crossing surface 52 between the air and the aluminum liner to a surface source write file, WSSA. Note that the surface must be specified as facets and that surfaces created with LIKE BUT, such as 12051.1, may not be used. Note that the name option will write a file with .w appended

MCNP6 I=inputfile name=outSSW.

will rename the WSSA file as outSSW.w

The resulting tally for surfaces 51, 52, and 53 is

1tally 1 nps = 100000 Neutrons crossing surfaces surface a is (51.1 51.2 51.3) surface b is (52.1 52.2 52.3) surface c is (53.1 53.2 53.3) angle : 0.0000E+00 1.0000E+00 total surface a 9.72549E-01 0.0079 1.99828E+00 0.0071 2.97083E+00 0.0073 b 1.29713E+00 0.0071 3.72160E+00 0.0062 5.01873E+00 0.0064 c 0.00000E+00 0.0000 2.41904E+00 0.0060 2.41904E+00 0.0060

The SSR run input file that utilizes the resulting surface source file written by the SSW run is shown in Example 2.16.

Example 2.16 SSR Run Input File Reading the Surface Source File Produced by Example 2.15

Air tank and aluminum casing 21 102 -1.0 -52 imp:n=0 $ air tank 22 103 -2.7 -53 52 imp:n=1 $ Al liner 23 0 53 imp:n=0 $ outside 52 RCC 10 10 -21 0 0 42 25 $ air cylinder 53 RCC 10 10 -22 0 0 44 40 $ aluminum liner m102 7014.80c .8 8016.80c .2 m103 13027.80c 1 print nps 200000 SSR OLD = 52.1 52.2 52.3 FC1:n Neutrons crossing surfaces F1:n (52.1 52.2 52.3) (53.1 53.2 53.3) C1 0 1 T FQ1 F C

Example 2.16 is the input file that utilizes the surface source produced by Example 2.14. It is run with an execution line that names the surface source file to be used.

MCNP6 I=exSSR name=outSSR. RSSA=outSSW.w

Note that the interior of the air tank is deleted, along with its surfaces and materials. The importance of the air (cell 21) is set to zero so that any neutron backscattering out of the aluminum liner is killed. This important albedo consideration is because all particles scattering from the liner back into the air were accounted for already in the surface source and must not be counted twice. The surface source read card

SSR OLD = 52.1 52.2 52.3

specifies the surfaces to be read. They may also have new names in the subsequent SSR run as long as they have the same geometric location.

The full SSR card is

SSR old=S1 S2 new=S3 S4 cel=C1 C2 col=n WGT=w TR=m PSC=x

where old and new surface names are given; for example,

SSR OLD = 52.1 52.2 NEW = 62.1 62.2 TR=123

reads the tracks crossing surfaces 52.1 and 52.2, which are started on identical surfaces 62.1 and 62.2 located elsewhere according to surface transformation TR123 (not shown). Surface 52.3 is not read.

For both the SSW and SSR, cel=C1 C2 is used to record all fission sites in a KCODE calculation. Additional options are available on the SSR card (MCNP manual section 3.3.4.8). The entry col=n filters n = −1/0/1 uncollided, all, collided particles. WGTcan renormalizes particle weights. TR maps the surface elsewhere, or it can be a distribution to repeat it many places in the subsequent run. PSC is used only for point detector (F5) tallies and DXTRAN variance reduction.

The surface source that is written in the initial run is described in the output file of the run that reads it. For Example 2.16, the SSR output file describes the surface source it is reading

1summary information from this use of surface-source file outSSW.w print table 170 number of original histories 100000 number of records read 677171 number of recorded histories 99271 number of histories read 99271 source multiplication factor 2.00000 total source weight 3.72160E+05 total rejected weight 0.00000E+00 fraction weight rejected 0.00000E+00

In Example 2.15, 100,000 histories were run, causing many fissions to produce 677,171 tracks crossing surface 52. The number of recorded histories is 99,271, indicating that 729 histories and the resulting fission progeny did not cross surface 52. The source multiplication factor of 2.0 is because NPS is 100,000 in the SSW Example 2.15 and 200,000 in the SSR Example 2.16. If NPS is different in the SSW writing and SSR reading runs, the tracks crossing the surface will be split or rouletted accordingly.

Particles or tracks written to the surface source file may be rejected if they are for a particle not listed on the subsequent run MODE specification, if they are filtered by col=n, if they are not on one of the listed OLD surfaces, or if they are outside the time or energy cutoffs of the subsequent problem.

The output of Example 2.16 summarizes the surface source read as

total number of particles accepted 677171 total weight accepted 3.72160E+05 fraction of total weight accepted 0.10000E+01

That is, 677171 tracks crossed surface 52, with a total weight of 3.72160E5. The summary table of Example 2.16 output states

source 1354342 3.7216E+00 8.1280E-01

There are 1354342 source tracks because 200000 source histories were requested (SSR run) from the initial 100000 SSW run: 2 × 677171 = 1354342. The average source weight is 3.7216E+00, and the average source energy is 8.1280E–01 MeV based on a source weight of 1.0 and 100000 histories in the initial SSW run. All other summary information and tallies in the output file are also based on the initial source weight of 1.0 and 100000 histories. The results of the current tally in the SSR run are

1tally 1 nps = 99271 Neutrons crossing surfaces number of histories used for normalizing tallies = 100000.00 surface a is (52.1 52.2 52.3) surface b is (53.1 53.2 53.3) angle : 0.0000E+00 1.0000E+00 total surface a 1.29752E+00 0.0066 0.00000E+00 0.0000 1.29752E+00 0.0066 b 0.00000E+00 0.0000 2.41864E+00 0.0062 2.41864E+00 0.0062

The SSR run results are statistically the same as SSW run results. For every 1 starting neutron, 1.297 neutrons cross surface 52 from the aluminum into the air, where they are killed. Then 3.7126 neutrons enter from the air on surface 52 in the SSW run and from the surface source in the SSR run. Then 2.419 cross surface 53 to escape the geometry.

In this example, the SSR run (Example 2.16) had the same result as the SSW run model of the entire geometry (Example 2.15), but the Figure of Merit (FOM) was 4.5 times higher—that is, the same result was achieved 4.5 times faster. The FOM is described in detail in Sect. 4.1 on Variance Reduction and is the sole measure of relative problem efficiency.

Following is an example of the SSW/SSR CEL option and how a KCODE criticality calculation can be on a simple cell that is then translated to multiple locations. A single uranium can has its KCODE criticality source neutrons written to a surface source file with the SSW CEL option in Example 2.17. Then this file is read in the four-can geometry with the SSR CEL option and translated to all four can locations in Example 2.18. The result ignores the small backscattered induced fission, but is a good approximation of running the KCODE calculation on the four-can geometry in Example 2.14.

Example 2.17 New Geometry Consisting of a Single Uranium Can from Example 2.14. The SSW CEL Option Records All Fission Neutrons from All Active Cycles of the Criticality Calculation in the One Uranium Can, Cell 24. Because the Can Is Surrounded by a Zero-Importance Cell, No Neutrons from Other Cans Enter. The Resulting Neutrons from Fission Sources Are Written to the SSW Surface Source Write File

Single uranium can in air and aluminum 24 101 -18.0 -50 imp:n=1 $ uranium can 23 0 50 imp:n=0 $ outside 50 RCC 0 0 -20 0 0 40 8 $ fuel cylinder m101 92238 -.8 92235 -.2 nlib=.80c print kcode 1e4 1.0 30 130 sdef axs=0 0 1 rad=d31 ext=d32 si31 0 8 sp31 -21 1 si32 -20 20 sp32 0 1 SSW CEL 24

Example 2.18 Four-Can geometry of Example 2.14 with the Fission Source from One-Can Geometry Example 2.17, Read with SSR CEL and Translated to All Four Can Positions. The Approximation Introduced by the NONU Card Is Discussed Below

Four uranium cans in air and aluminum 11 101 -18.0 -51 imp:n=1 $ uranium can 12 like 11 but trcl=(20 0 0) imp:n=1 $ uranium can 13 like 11 but trcl=( 0 20 0) imp:n=1 $ uranium can 14 like 11 but trcl=(20 20 0) imp:n=1 $ uranium can 21 102 -1.0 -52 #11 #12 #13 #14 imp:n=1 $ air tank 22 103 -2.7 -53 52 imp:n=1 $ Al liner 23 0 53 imp:n=0 $ outside 51 RCC 0 0 -20 0 0 40 8 $ fuel cylinder 52 RCC 10 10 -21 0 0 42 25 $ air cylinder 53 RCC 10 10 -22 0 0 44 40 $ aluminum liner m101 92238 -.8 92235 -.2 nlib=.80c m102 7014.80c .8 8016.80c .2 m103 13027.80c 1 print nps 100000 NONU SSR CEL 24 TR=D40 SI40 L 0 42 43 44 SP40 1 3R TR42 20 0 0 TR43 0 20 0 TR44 20 20 0 FC1:n Neutrons crossing surfaces F1:n (51.1 51.2 51.3) (52.1 52.2 52.3) (53.1 53.2 53.3) C1 0 1 T FQ1 F C FC4 Fissions and neutron production in uranium cells F4:n 11 12 13 14 T FM4 (-1 101 -6) (-1 101 -6 -7) FQ4 F M SD4 1 4r

The surface tally output of Example 2.18 follows:

1tally 1 nps = 1296681 Neutrons crossing surfaces number of histories used for normalizing tallies = 1000000.00 surface a is (51.1 51.2 51.3) surface b is (52.1 52.2 52.3) surface c is (53.1 53.2 53.3) angle : 0.0000E+00 1.0000E+00 total surface a 2.22166E-01 0.0022 4.40912E-01 0.0013 6.63078E-01 0.0015 b 2.88259E-01 0.0019 8.16748E-01 0.0008 1.10501E+00 0.0011 c 0.00000E+00 0.0000 5.27283E-01 0.0006 5.27283E-01 0.0006

The surface tally output of Example 2.14 follows:

1tally 1 nps = 99770 Neutrons crossing surfaces number of histories used for normalizing tallies = 1000591.00 surface a is (51.1 51.2 51.3) surface b is (52.1 52.2 52.3) surface c is (53.1 53.2 53.3) angle : 0.0000E+00 1.0000E+00 total surface a 2.10443E-01 0.0076 4.21114E-01 0.0052 6.31556E-01 0.0057 b 2.95602E-01 0.0068 8.18688E-01 0.0039 1.11429E+00 0.0045 c 0.00000E+00 0.0000 5.21861E-01 0.0036 5.21861E-01 0.0036

Example 2.14 used a KCODE fission source. Example 2.18 used the SSR surface source instead with the CEL option and ran eight times faster and demonstrated its efficacy by achieving nearly the same result.

The results are comparable. The differences are because

  • induced neutrons in each can from neutrons scattering in from other cans or reflected are ignored. Consequently the surface 51 tally, surface “a,” is 5% lower;

  • neither KCODE criticality calculation was fully converged, and the relative errors of the SSR calculation, Example 2.18, assume the surface source file from Example 2.17 was perfectly converged;

  • the single can KCODE calculation Example 2.17 included neutrons from cycles that were not fully settled; and

  • the four-can KCODE calculation Example 2.14 had different numbers of fissions in each can even though, by symmetry, they should have been the same.

Other points of interest are as follows:

  • The NONU option was required in Example 2.18 because all fission neutrons in the single can were already accounted for in the surface source. However, this approach is approximate because induced fissions from in-scattering from other cans or back scatter from the aluminum liner are neglected. To avoid this incorrect albedo neglecting induced fissions, the SSR should use the sources on the surface rather than the CEL option on the SSR card and then allow fissions in the cells by omitting the NONU card.

  • The SSR cell name is CEL = 24 even though there is no cell 24 in the SSR run of Example 2.18 because of SSW CEL 24 in Example 2.17.

  • The source can be moved and duplicated anywhere with TR on the SSR card.

  • Repeating the single-can source from Example 2.17 in the four can locations of Example 2.18 ensures that the fissions are the same in each can, which they should be from symmetry.

  • Use of the CEL option on SSW/SSR is usually approximate. But in cases where there are many fissionable sources that are only loosely coupled to each other the SSW/SSR with the CEL option can be a very effective approximation. In safeguards science applications include barrels of nuclear waste, arrays of tanks of enriched or depleted uranium, and in some cases, reactor fuel pins.

The advantages of using SSW and SSR in safeguards runs are as follows:

  • Calculations may be faster and more efficient when separated into several parts.

  • Calculations that require analog sampling, such as multiplicity counters and pulse-height detectors, cannot use variance reduction but can be separated into parts that can take advantage of variance reduction.

  • Parameter studies, particularly small perturbations with the MCNP perturbation (PERT) capability, can be applied to the subsequent SSR problem without having to run the entire SSW run initial problem.

  • Calculations with a similar source, such as a bremsstrahlung electron accelerator, may be separated from subsequent SSR runs, modeling a variety of different targets.

Some disadvantages of using SSW and SSR in safeguard runs are

  • the surface source write/read files may be very large;

  • care must be taken for albedo effects—re-entering particles should not be duplicated or ignored;

  • relative errors are not propagated; SSR assumes the SSW file has zero relative error and, thus, the SSR relative errors are underestimated;

  • there may be streaming effects—the only directions sampled in the SSR run are those written in the SSW run, which may miss directions of interest; and

  • point detectors and DXTRAN variance reduction require a user-input approximate surface direction (PSC option on the SSR card); this approximation is needed; particle tracks read from an SSR file have no recollection of the source or collision preceding the surface crossing.

2.3.6 Checking sources

The MCNP output file, OUTP, is the first way to check the source. All too often, what the user thinks was requested in the input file is not what the code thinks it was told! Print Tables 10, 110, 128, and 170 specifically describe the source. The problem summary and weight balance tables also provide important source information. The descriptions in Sect. 2.4.1 should be reviewed every time the source is modified.

Where possible, a tally should be constructed that confirms the source. The SCD and SCX special tally treatments are often useful. Tallies are described in Sect. 2.4.

Finally, it is possible to plot source distributions with the TMESH mesh tally and special options of the FMESH mesh tally, which are described in Sect. 2.6.

2.4 Output and Tallies

2.4.1 Output Files

MCNP output files, OUTP, are generated by every MCNP run. The basics of MCNP generated files are given in Sect. 2.1.1.

The output file is organized into Print Tables. If the PRINT card is in the input file or the option “PRINT” appears on the MCNP execution line, then the maximum print of the OUTP is made. Otherwise, only the minimum set tables are provided (default). The PRINT card has the form

PRINT 110

The above input card will print the minimum set of tables plus table 110.

PRINT -85 -86 -128 -70 -98

The above input card (MCNP code manual section 3.3.7.2.1) will print the possible tables except print tables 85, 86, 128, and 70. Tables 85 and 86 are the electron and charged particle range tables, which are very lengthy and unchanged from calculation to calculation unless materials are changed. Print Table 128 is the “lattice universe map,” which not only requires a lot of storage but slows the calculation by recording every source, entering, and collision in every cell and every repeated structures and lattice sub-cell. See the example in Sect. 2.1.7, “Filled Cells: Universes.” Print Table 70 gives the surface coefficients of quadratic and macrobody surfaces and is of interest in only some circumstances. Print Table 98, giving basic constants, takes up little space but is the same for every run.

All of the OUTP output file should be reviewed the first time a problem is run, and then the parts that change should be reviewed every time the input file is changed or the files that the MCNP code is reading—such as the RSSA file from the surface source—are changed. Some of the more important Print Tables are described as follows.

Source Check: Optional Print Table 10 lists all SDEF parameters, showing the defaults when they have not been set in the input file. It also describes each source distribution and its probability, along with the KCODE criticality source parameters and what is read from the surface source read file, RSSA. See the example in Sect. 2.3.5, “Surface Source Read and Write (SSR, SSW).”

Variance-Reduction Summaries: Print Table 20 describes the MCNP interpretation of the input file weight windows when present. Print Table 62 describes forced collision and exponential transform variance-reduction input. The Problem Summary, particle activity by cell (Print Table 126), and particle weight balance by cell (Print Table 130) also contain variance-reduction activity summaries. Print Tables 180 and 190 summarize the weight window generator for variance reduction. Print Table 190 summarizes weight window generation from multigroup fluxes from alternative codes, if used. Generated cell-based weight windows, if any, are provided at the very end of the output file.

Geometry/Mass Check: Print Tables 40, 50, and 60 describe the cell volumes, masses, nuclide fractions, and importances. Masses should be checked against the physical masses being modeled. Track length tallies (F4, F6, F7) are normalized by volumes and masses, and these normalizations should be verified. In repeated structures or lattices, only the master cell volume and mass are calculated and a multiple of them must be used for the tally normalizations.

Materials Check: Print Table 100 is essential reading. About 10% of all user errors are due to not having the desired materials. Print Table 100 lists what the MCNP code is using, which can change from platform to platform and from time to time depending on which computing platform is hosting the MCNP calculation. Also, the intuitive specification of nuclides—such as 6012 for carbon neutron cross sections—causes the MCNP code to use decades-old data rather than the correct 6000 elemental specification for modern 12C. A listing of available data is in LA-UR-20709 [4] for the MCNP6.2 release, but the MCNP code uses whatever is in the user’s XSDIR data directory, and Print Table 100 should always be checked to know which data are actually used.

Source Check: Optional Print Table 110 lists the basic parameters of the first 50 histories: position (x,y,z), cell, surface, direction (u,v,w direction cosines), energy, weight, and time. If the source consists of multiple particle types, then the particle type is provided. If repeated structures and universes are used, then the source position, cell, and direction are provided at each level of the universe, along with lattice index, if appropriate. The energy, weight, and time columns are not shown in the following, which is from Example 2.13:

print table 110 nps par x y z cell lattice[i j k] surface u v w 1 p 1.000E+01 3.000E+00 0.000E+00 15 0 2.033E-01 7.180E-01 6.657E-01 2 sf -2.104E+00 2.156E+00 -1.874E+00 14 7.962E-01 -4.029E-01 -4.513E-01 sf -1.039E-01 2.156E+00 1.258E-01 21[-1 0 -1] 7.962E-01 -4.029E-01 -4.513E-01 sf -1.039E-01 2.156E+00 1.258E-01 1 0 7.962E-01 -4.029E-01 -4.513E-01 3 p 1.000E+01 3.000E+00 0.000E+00 15 0 -5.363E-01 6.280E-01 -5.638E-01 4 n 1.334E+01 1.033E+01 -1.825E+00 11 8.493E-01 4.573E-01 2.637E-01 n 3.446E-01 1.033E+01 -9.270E-02 22[-4 1 0] 8.493E-01 4.573E-01 2.637E-01 n 3.446E-01 1.033E+01 -9.270E-02 1 0 8.493E-01 4.573E-01 2.637E-0

Problem Summary: This essential summary of particle and energy creation and loss balance should be reviewed in every calculation. It is provided for each particle type and provides a wealth of important information. The following problem summary is from Example 2.13:

neutron creation tracks weight energy neutron loss tracks weight energy (per source particle) (per source particle) source 898096 1.0000E+00 1.8331E+00 escape 1039960 1.0468E+00 1.2701E+00 nucl. interaction 0 0. 0. energy cutoff 0 0. 0. particle decay 0 0. 0. time cutoff 0 0. 0. weight window 0 0. 0. weight window 0 0. 0. cell importance 0 0. 0. cell importance 0 0. 0. weight cutoff 0 2.4484E-03 9.2821E-05 weight cutoff 13869 2.4638E-03 8.1124E-05 e or t importance 0 0. 0. e or t importance 0 0. 0. dxtran 0 0. 0. dxtran 0 0. 0. forced collisions 0 0. 0. forced collisions 0 0. 0. exp. transform 0 0. 0. exp. transform 0 0. 0. upscattering 0 0. 1.7092E-08 downscattering 0 0. 8.9672E-01 photonuclear 0 0. 0. capture 0 6.5082E-02 1.0656E-02 (n,xn) 765 7.4328E-04 4.8758E-04 loss to (n,xn) 382 3.7116E-04 3.1146E-03 prompt fission 259293 1.8546E-01 3.7613E-01 loss to fission 105888 7.5412E-02 2.9833E-02 delayed fission 1945 1.4743E-03 7.7146E-04 nucl. interaction 0 0. 0. prompt photofis 0 0. 0. particle decay 0 0. 0. tabular boundary 0 0. 0. tabular boundary 0 0. 0. tabular sampling 0 0. 0. elastic scatter 0 0. 0. total 1160099 1.1901E+00 2.2105E+00 total 1160099 1.1901E+00 2.2105E+00 number of neutrons banked 472911 average time of (shakes) cutoffs neutron tracks per source particle 1.2917E+00 escape 3.4851E+07 tco 1.0000E+33 neutron collisions per source particle 1.0777E+01 capture 3.2604E+07 eco 0.0000E+00 total neutron collisions 9678987 capture or escape 3.4719E+07 wc1 -5.0000E-01 net multiplication 1.1150E+00 0.0005 any termination 3.2597E+07 wc2 -2.5000E-01 photon creation tracks weight energy photon loss tracks weight energy (per source particle) (per source particle) source 400222 4.4563E-01 6.7254E+01 escape 1447030 1.7452E+00 5.0610E+01 nucl. interaction 0 0. 0. energy cutoff 0 0. 3.3191E-03 particle decay 0 0. 0. time cutoff 0 0. 0. weight window 0 0. 0. weight window 0 0. 0. cell importance 0 0. 0. cell importance 0 0. 0. weight cutoff 0 0. 0. weight cutoff 0 0. 0. e or t importance 0 0. 0. e or t importance 0 0. 0. dxtran 0 0. 0. dxtran 0 0. 0. forced collisions 0 0. 0. forced collisions 0 0. 0. exp. transform 0 0. 0. exp. transform 0 0. 0. from neutrons 484434 9.9354E-01 2.0178E+00 compton scatter 0 0. 3.9870E+00 bremsstrahlung 8435674 9.6553E+00 1.5243E+01 capture 12018341 1.4424E+01 3.3933E+00 p-annihilation 624882 7.1252E-01 8.1705E-01 pair production 312441 3.5626E-01 2.7897E+01 photonuclear 0 0. 0. photonuclear abs 0 0. 0. electron x-rays 0 0. 0. loss to photofis 0 0. 0. compton fluores 0 0. 0. muon capt fluores 0 0. 0. 1st fluorescence 3241292 3.9805E+00 5.3427E-01 2nd fluorescence 591308 7.3783E-01 2.4329E-02 cerenkov 0 0. 0. (gamma,xgamma) 0 0. 0. tabular sampling 0 0. 0. prompt photofis 0 0. 0. total 13777812 1.6525E+01 8.5891E+01 total 13777812 1.6525E+01 8.5891E+01 number of photons banked 10136298 average time of (shakes) cutoffs photon tracks per source particle 1.5341E+01 escape 3.1190E+08 tco 1.0000E+33 photon collisions per source particle 1.6852E+01 capture 2.3323E+08 eco 1.0000E-03 total photon collisions 15134406 capture or escape 2.4173E+08 wc1 -5.0000E-01 any termination 2.4404E+08 wc2 -2.5000E-01 computer time so far in this run 1.29 minutes maximum number ever in bank 94 computer time in mcrun 1.24 minutes bank overflows to backup file 0 source particles per minute 8.0335E+05 random numbers generated 337454313 most random numbers used was 9187 in history 689906

Physics Summary: Optional Print Table 126 provides the activity in each cell for each particle type. The eight columns of information can also be plotted using the PAC, N, and PAR commands on the right-hand side of the interactive geometry plotter. The eight columns of information are:

  • PAC1: tracks entering (including source particles and secondary particles)

  • PAC2: population (does not include re-entrant particles)

  • PAC3: collisions

  • PAC4: collisions * weight

  • PAC5: number weighted average energy

  • PAC6: flux weighted average energy

  • PAC7: average track weight

  • PAC8: average track mean free path (MFP)

These quantities can be plotted after a run using the RUNTPE continuation file:

mcnp6 i=inp r=runtpefile z

or during the course of a run using the interrupt:

<ctrl-c><m>

Once the MCPLOT prompt appears, type “PLOT” to get the geometry plot. The command “LABEL 0 1 PAC3” in the bottom left block will put the number of collisions in each cell of the plot. Alternatively, click PAC, N, N, N, and L2. Click PAR to get the correct particle type if there is more than one.

Repeated Structures Physics Summary: Optional Print Table 128 provides the activity in each repeated structure/lattice element. A sample from Example 2.13 follows:

1neutron activity in each repeated structure / lattice element print table 128 source entering collisions path 12371 23463 7105 1 < 21[-3 0 -3] < 14 0 36450 2045 2 < 21[-3 0 -3] < 14 0 53043 39364 3 < 21[-3 0 -3] < 14 12338 29259 9835 1 < 21[-2 0 -3] < 14 0 49460 2789 2 < 21[-2 0 -3] < 14 0 75996 65336 3 < 21[-2 0 -3] < 14 12467 29552 10042 1 < 21[-2 0 3] < 14 0 50204 2786 2 < 21[-2 0 3] < 14 0 77396 72336 3 < 21[-2 0 3] < 14 2677 6333 1304 1 < 22[ 1 4 0] < 11 0 7440 360 2 < 22[ 1 4 0] < 11 0 11781 6326 3 < 22[ 1 4 0] < 11 0 1125238 0 15 0 0 0 16 898096 18269879 9678987

Cell Weight Balance Summary: Optional Print Table 130 is the particle weight balance in each cell. The particle weight information from the summary table is provided for each individual cell in the problem.

Nuclide Activity by Cell Summary: Optional Print Table 140 is the activity of each nuclide in each cell, which includes the total collisions, collisions by weight, weight lost to capture, weight gain by fission, weight gain by (n,xn), and secondary particle production information. In coupled neutron-photon problems, it also includes several summary tables of the photoatomic activity of each nuclide in each cell, per source particle.

Source Check: Optional Print Table 170 describes each source distribution, how often it was sampled, and how often it was expected to be sampled. For surface source SSR runs, the RSSA file utilization is described.

The output file concludes with the tallies, which are described in the rest of Sect. 2.4. The tallies are user-specified and so general that they can be used to answer almost any question of what is going on in the calculation. For each tally, there is a tally fluctuation chart at the end of the output file and a “yes” or “no” answer to whether the tally passes each of ten statistical tests for convergence. There is also an optional more detailed statistical analysis of tally convergence (Print Table 160), optional printed plots of the probability density function (Print Table 161), and cumulative normed and unnormed tally density function plots (Print Table 162).

2.4.2 MCNP Estimators and Tally Types

The Monte Carlo method does not give answers. It provides estimates with statistical errors. In the MCNP code, the statistical errors are relative error fractions—one standard deviation of the mean divided by the mean value. The four possible estimators are as follows:

  • Track length estimator: flux = weight * track length/volume. The MCNP code uses the track length estimator for cell-averaged flux, heating, keff criticality, and reaction rates.

  • Surface estimator: the MCNP code uses the surface estimator for the current or flux of particles crossing a surface, energy balance, and pulse counts in cells.

  • Collision estimator: the MCNP code uses the collision estimator for reaction rates, keff criticality, problem summary and energy deposition.

  • Pseudo-deterministic or next-event estimators: the MCNP code uses next-event estimators for flux and reaction rates.

Flux vs. Fluence: Fluence is the time-integrated flux. The MCNP code is fully time-dependent even if time information is not requested in the tallies. If the source units are particles per second, then the tally quantity is flux in units of particles/cm2/s. If the source is one source particle over all time, then the tally quantity is fluence, or particles/cm2; that is, the tally units are determined by the source. If the source is in particles/s, then the tally will be per second.

Tally Types: The MCNP tally types are:

Tally

Type

Estimator

Fn Units

*Fn Units

F1:

Surface Current

W

#

MeV

F2:

Surface Fluence

W/(|μ|*A)

#/cm2

MeV/cm2

F4:

Cell Fluence

/V

#/cm2

MeV/cm2

F5:

Detector Fluence

Wp(μ)eσx/2πR2

#/cm2

MeV/cm2

F6:

Energy Deposition

WλσT(E)H(E)ρa/ρg

MeV/g

jerks/g

F7:

Fission Energy Deposition

Wλσf(E)a/ρg

MeV/g

jerks/g

F8:

Pulse Height

WS in bin EDW/WSs

pulses

MeV

The symbols are as follows:

  • W = particle weight; Ws = source weight

  • E = particle energy (MeV); ED = Energy deposited (MeV)

  • μ = cosine of angle between surface normal and trajectory

  • A = surface area (cm2); V = cell volume (cm3)

  • λ = track length (cm)

  • p(μ) = probability density function: μ = cosine of angle between particle trajectory and detector

  • σ = total mean free path to the detector

  • R = distance to the detector (cm)

  • σT(E) = microscopic total cross section (barns)

  • H(E) = heating number (MeV/collision)

  • ρa = atom density (atoms/barn-cm)

  • ρg = gram density (g/cm3)

  • M = cell mass (g) = ρa * V

  • σf(E) = microscopic fission cross section (barns)

  • Q = fission heating Q-value (MeV)

The historical (Manhattan Project) unit “jerk” is equal to 109 Joules or 6.241 × 1021 MeV, or approximately one quarter of a metric ton of TNT.

  • Track length estimator: F4, F6, F7, keff, reaction rates

  • Surface estimator: F1, F2, F8, energy balance

  • Collision estimator: F6, keff, reaction rates, energy balance

  • Next-event estimator: F5 flux, reaction rates

2.4.3 Basic Tally Format

2.4.3.1 The Eight Tally Dimensions: FDUSMCET

There are eight dimensions of MCNP tallies:

  • F: surface numbers (F1, F2), cell numbers (F4, F6, F7, F8), and point detector parameters (F5);

  • D: surface flagging (SF), cell flagging (CF), or direct/total detector;

  • U: user bins;

  • S: segment bins;

  • M: multiplier bins;

  • C: angle bins;

  • E: energy bins;

  • T: time bins.

2.4.3.2 F: Surface and Cell Tallies

The basic tally format for cell and surface tallies (not F5 detector tallies) is

Fn:<pl> S1 S2 . . . Sk Fn:<pl> S1 (S2 . . . S3) (S4 . . . S5) S6

where n is a positive integer with the following meaning. The last digit of n is the tally type. If it ends in 1, then it is an F1 tally. For example, F111:N. The particle type, <pl>, is N for neutrons, P for photons, E for electrons, H for protons, etc. The Sn entries are surfaces for F1 and F2 and cells for F4, F5, F6, F7, and F8. Parentheses can be used to sum over several surfaces or get the average over several cells. If the surfaces are macrobodies, the facets must be specified, as in Example 2.14, Example 2.15, and Example 2.18.

F1:n (51.1 51.2 51.3) (52.1 52.2 52.3) (53.1 53.2 53.3) F4:n 11 12 13 14 T

The T at the end of the F4 tally averages over all cells, which is equivalent to

F4:n 11 12 13 14 (11 12 13 14)

In some earlier versions of the MCNP code, surfaces—and particularly macrobody facets—that are eliminated for being the same as others cannot tallied upon.

2.4.3.3 F: Tallies in Lattices and Repeated Structures

Tallying in lattices and repeated structures requires specifying not only the cells or surfaces to be tallied but also the paths to them through various levels of universes. Example 2.19, which follows, uses the lattice/repeated structure geometry of Example 2.13. This geometry is illustrated in Fig. 2.9. Surfaces 51, 52, and 53 have been added so that both the surface tally 901 and energy deposition tally 906 can illustrate tally segmenting.

Example 2.19 Repeated Structures/Lattice Geometry from Example 2.12 and Example 2.13 Illustrating Tallies in Repeated Structures and Lattices and Segment Tallies. (Note: The MCNP6 Manual (Appendix C) Recommended Parameters for the 252Cf Neutron Energy Watt Spectrum Are 1.18 (MeV) and 1.03419 (MeV−1)).

Square and hex lattices 1 313 -10.44 -1 u=1 imp:n=1 $ UO2 2 201 -6.55 -2 1 u=1 imp:n=1 $ Zr clad 3 609 -1.00 2 u=1 imp:n=1 $ H2O 21 111 -4.00 -31 lat=1 u=11 imp:n=1 $ Square lattice fill = -3:3 0:0 -3:3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 1 1 1 1 1 1 22 111 -4.00 -32 lat=2 u=12 imp:n=1 $ Hexagonal lattice fill = -5:5 -4:4 0:0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 12 1 1 1 1 1 1 0 0 0 0 14 0 -4 fill 11 imp:n=1 $ Square lattice box 11 like 14 but TRCL=(20 0 0) fill 12 $ Hexagonal lattice box 15 0 #14 #11 -20 imp:n=1 $ outside boxes 16 0 20 imp:n=0 $ outside world 1 rcc 0 0 0 0 12 0 .5 $ Right circular cylinder 2 rcc 0 0 0 0 12 0 .6 $ Right circular cylinder 4 rpp -7 7 0 12 -7 7 $ Right parallelepiped 14x12x14 31 rpp -1 1 0 12 -1 1 $ Right parallelepiped 2x12x2 32 rhp 0 0 0 0 12 0 1 0 0 $ Right hex 12x1 20 sph 0 0 0 50 $ Sphere containing lattice boxes 51 PX 0 52 PX 8 53 PX 12 MODE N P m0 nlib=.70c plib=.04c elib=.03e m201 40090 -0.505239 40091 -0.110180 40092 -0.168413 40094 -0.170672 40096 -0.027496 nlib=80c m313 92235 -0.02759 92238 -0.85391 8016 -0.11850 M609 1001 2 8016 1 MT609 LWTR.10t M111 6000 1 1001 1 7014 1 MT111 POLY.10t BENZ.10t GRPH.10t SDEF PAR=D42 CEL=FPAR=D51 POS=FPAR=D53 RAD=FPAR=D54 EXT=FPAR=D55 AXS=FPAR=D56 TME=FPAR=D57 ERG=FPAR=D58 SI42 L SF N P P P P $ PAR SP42 .3 .3 .1 .1 .1 .1 $ PAR DS51 S 46 47 48 48 48 48 $ CEL DS53 L 0 6 0 0 6 0 10 3 0 10 3 0 10 7 0 10 7 0 $ POS DS54 S 43 43 0 0 0 0 $ RAD DS55 S 44 44 0 0 0 0 $ EXT DS56 S 45 45 0 0 0 0 $ AXS DS57 L 0 1e8 2e8 3e8 4e8 5e8 $ TME DS58 S 0 92 93 94 95 96 $ ERG SC46 Spontaneous fission cells in XYZ lattice SI46 L (1<21[-2:3 0 3]<14) (1<21[-3:3 0:0 -3:2]<14) SP46 1 5R 1 41R SC47 Neutron source cells in HEX lattice SI47 L (1<22[-5 3 0]<11) (1<22[-4 1:4 0]<11) (1<22[ -3 -1:4 0]<11) (1<22[-2 -3:4 0]<11) (1<22[-1:1 -4:4 0]<11) (1<22[ 2 -4:3 0]<11) (1<22[ 3 -4:1 0]<11) (1<22[ 4 -4:-1 0]<11) (1<22[ 5 -4:-3 0]<11) SP47 .5 .5 1 1 .5 .5 1 3R .5 .5 1 5R .5 .5 1 6R .5 .5 1 6R .5 .5 1 6R .5 .5 1 5R .5 .5 1 3R .5 .5 1 1 .5 .5 .5 SC48 Point photon source cell in water between lattices SI48 L 15 SP48 1 c particle-dependent radius, extent, axis SI43 0 .5 $ Radii SP43 -21 1 $ Radii SI44 -6 6 $ Extent SP44 -21 0 $ Extent SI45 L 0 1 0 $ Axis SP45 1 $ Axis c energies dependent upon time, position, and particle type SP92 -3 1.175 1.0401 $ 252Cf Frohner Watt parameters SI93 50 59 $ Photon, XYZ water, TME 2e8, energies SP93 0 1 SI94 60 69 $ Photon, XYZ water, TME 3e8, energies SP94 0 1 SI95 L 70 $ Photon, XYZ water, TME 4e8, energies SP95 1 SI96 L 80 $ Photon, XYZ water, TME 5e8, energies SP96 1 nps 100000 PRINT -162 -85 -86 -140 c ***** Tallies ************** FC901 Segmented surface tally F901:n (4.1 4.2) (11004.1 11004.2) (2.1<21<14) (2.1<22<11) (11004.1 4I 11004.6) (2.1<21[0 0 0]<14) FS901 -52 -53 T $ 4 segment bins FQ901 F S $ Print cells down, segments across FC906 Segmented energy deposition tally +F906 (1<21<14) (1<22<11) (21<14) (22<11) (1<21[-3:3 0:0 -3:3]<14) (1<22[-5:5 -4:4 0:0]<11) FS906 -51 -52 -53 T $ 4 segment bins FQ906 F S $ Print cells down, segments across SD906 1 759R $ MeV units, 102 cells x 5 segments

Tallying in lattices and repeated structures is illustrated in Example 2.19 with the surface tally F901 and cell tally F906.

The “<” separates different universe levels. The cells or surfaces on each side must be in the universe at the next level and not in the same universe.

The F901 tally provides the format for F1 and F2 type surface tallies in lattices or repeated structures. In this example, it counts neutrons in the following bins:

  • Bin 1 (4.1 4.2): either facet 1 or 2 of surface 4;

  • Bin 2 (11004.1 11004.2): either facet 1 or 2 of surface 4 translated to cell 11;

  • Bin 3 (2.1<21<14): any of the 49 appearances of surface 2.1 u=1 filling rectangular lattice cell 21 u=11, which fills cell 14 u=0;

  • Bin 4 (2.1<22<11): any of the 99 appearances of surface 2.1 u=1 filling hexagonal lattice cell 22 u=12, which fills cell 11 u=0;

  • Bin 5 (11004.1 4I 11004.6): surfaces 4.1 and 4.2 of surface 4 translated to cell 11. In addition, bin 5 also counts all neutrons crossing surfaces 1.2, 1.3, 4.5, and 4.6. because these surfaces are the same as 11004.3, 11004.4, 11004.5, and 11004.6. Further, the particles that cross surfaces 1.2, 1.3, 4.5, and 4.6 include particles that cross these surfaces anywhere in the problem and not just in cell 11; and

  • Bin 6 (2.1<21[0 0 0]<14): surface facet 2.1 u=1 filling lattice element 0,0,0 in cell 21 u=11 filling cell 14 u=0.

The +F906 tally has no particle type because +F6 energy deposition tallies include the energy deposition from all particle types in the problem from whichever estimator they use—either track length or collision. The F906 tally provides the format for F4, F6, and F7 type cell tallies in lattices or repeated structures. Again, the “<” separates universe levels; in this example, it tallies the energy deposition in the following bins from:

  • Bin 1 (1<21<14): the sum of all 49 occurrences of cell 1 u=1 filling cell 21 u=11 filling cell 14 u=0;

  • Bin 2 (1<22<11): the sum of all 99 occurrences of cell 1 u=1 filling cell 22 u=12 filling cell 11 u=0;

  • Bin 3 (21<14): the 1 occurrence of cell 21 u=11 filling cell 14 u=0;

  • Bin 4 (22<11): the 1 occurrence of cell 22 u=12 filling cell 14 u=0;

  • Bins 5–53 (1<21[-3:3 0:0 -3:3]<14): each of the 49 occurrences of cell 1 u=1 filling cell 21 u=11 filling cell 14 u=0; and

  • Bins 54–152 (1<22[-5:5 -4:4 0:0]<11): each of the 99 occurrences of cell 1 u=1 filling cell 22 u=12 filling cell 11 u=0.

Note that bins 3 and 4, (21<14) and (22<11), are needed because the (1<21[-3 0 3]<14) bin in which cell 21 is filled with itself and the (1<22[-5 4 0]<11) bin in which cell 22 is filled with itself have zero tallies because cell 1 is not in these lattice positions.

The U=n notation (not shown) may also be used at any level when there are multiple cells at that level, but this capability does not work in some cases, such as the above (see the MCNP code manual section 3.3.5.1.4).

2.4.3.4 C: Angle Bins and FC Tally Comments

Surface tallies can have angle bins and nearly all tallies can have time and/or energy bins. Also, a descriptive tally title is strongly encouraged. From Example 2.14, Example 2.15, and Example 2.18:

FC1:n Neutrons crossing surfaces F1:n (51.1 51.2 51.3) (52.1 52.2 52.3) (53.1 53.2 53.3) C1 0 1 T

The title on the tally comment card, FC1, has the description “Neutrons crossing surfaces” everywhere in the output or in a plot where tally 1 occurs.

The angle of particles that cross the surface is specified by the C1 cosine card. The cosine bins are −1 < μ < 0, 0 < μ < 1, and −1< μ < 1 from the “T” option, which is the total over all other bins. The cosines are relative to the surface outward normal. Thus, the cosine bins are everything crossing inward, everything crossing outward, and everything crossing in either direction. If *C1 is used, then the angles are in degrees, and the specification goes from 180°. For example, one-degree angle bins would be

*C1 179 177I 0 T

The lower limit of the bins, 180°, is not specified. Between 179° and 0°, there are 177 interpolated points.

2.4.3.5 E, T: Energy and Time Bins, Default Bins, and FQ and TF Format Control

Energy bins are specified in MeV with the upper energy bound, as in Example 2.20 for the heating tally F16:

e16 0.2 48I 10

The total over all energy or time—the “T” option—is the default for time and energy bins and can be turned off with the “NT” option. The above is the same as

e16 0.2 48I 10 T

Example 4.1 shows a number of capabilities:

t0 1e2 263ilog 2e5 e0 1e-4 1e-3 1e-2 1e2 fq0 T E

The time bins are in shakes (10−8 s) and go from 1 μs to 2 ms with logarithmic interpolation. The default energy bins are one bin over all energy; the default time bins are one bin over all time. The E0 and T0 change the default for all tally time and energy bins unless they are set for specific tallies, as in e16 above.

By default, energy bins are printed vertically in the output file and time bins horizontally. The FQ0 changes the default for all tallies to print the time bins down and the energy bins across. These are then in blocks of C bins, in blocks of M bins, in blocks of S bins, in blocks of U bins, in blocks of D bins, and in blocks of F bins. The FQ changes the order of printing (see MCNP manual section 3.3.5.6).

From Example 2.14, Example 2.15, and Example 2.18:

FQ1 F C FQ4 F M

Tally card FQ1 changes the default from

F D U S M C E T

to

FQ1 D U S M E T F C

Here, the output file prints C bins across and F bins down within blocks of T bins within blocks of E bins, etc.

Tally card FQ4 changes the default to

FQ1 D U S C E T F M

Here, the output file prints M bins across and F bins down within blocks of T bins within blocks of E bins, etc.

The fq0 T E input of Example 4.1 (in Sect. 4.1.2) prints energy bins across and time bins down.

The eight tally dimensions are also used with the tally fluctuation chart control, TF, which points to which bin of a tally is to be used for the tally fluctuation chart at the end of the output. The TF card bin is used to obtain the statistical analysis of a quantity of real interest in the simulation (for example, the user may not be interested in the last time bin of a problem) but also is crucial for the weight window generation (Sect. 4.1) target tally bin. The default TF entries are

TFn 1 1 L L 1 L L L

These entries point to the 1st F (cell, surface, or detector) bin, 1st D (flagging) bin, last U (user) bin, last S (segment) bin, 1st M (multiplier) bin, last C (cosine) bin, last E (energy) bin, and last T (time) bin.

In Example 3.4:

t4 10e8 98i 1000e8 tf4 7j 7

The above TF4 card uses the default F, D, U, S, M, C, E bins but the 7th T bin instead of the last 101st total bin.

2.4.3.6 M: Multiplier Bins and SD Divisors

The multiplier bins, which are the M dimension of tallies, are often used to convert flux tallies to reaction rates, were described in Example 2.14:

FC4 Fissions and neutron production in uranium cells F4:n 11 12 13 14 T FM4 (-1 101 -6) (-1 101 -6 -7) FQ4 F M SD4 1 4r

Tally 4 in Example 2.14 multiplies the flux in each fissionable can by ρσf and ρνσf to get the number of fissions and number of fission neutrons produced. The basic format of the tally multiplier (FM) card is

FMn (C1 m1 R1) (C2 m2 R2) . . . T

  • n = tally number

  • Ci = multiplicative constant (if Ci= −1, for type 4 tallies only, the tally is multiplied by the atom density of the cell, ρa)

  • mi = material number identified on an Mm card

  • Ri = a combination of ENDF reaction numbers

The reaction number values follow the ENDF (evaluated nuclear data file) standard and have the special MCNP values, as shown in Table 2.1:

Table 2.1 The Principal MCNP Reaction Number (MT) Values (PN is photonuclear; rxn is reaction). The S(α,β) reactions are: 1 = total, 2 = elastic, 4 = inelastic. The neutron cross sections do not include thermal effects. The photon cross sections do not include form factors for coherent and incoherent scattering

Other examples are

F2:N 1 2 $ 36 tally bins FM2 (1.0) (2.0) (3.0) $ Constant multipliers E2 .5 1 2 4 10 T $ Energy bins F4:N (1 2) 3 T$ 6 tally bins FM4 (-1 1 -6 -7) $ Track-length estimate of keff (-1 2 1 -4)$ Neutron Heating (MeV/cm3)

The SD4 1 4r in the above from Example 2.14 changes the units by overriding the normal tally divisor. The SD card replaces the usual divisor for F (cell/surface) and S (segment) bins. The usual divisors for F2, F4, F6, and F7 are Area, Volume, Mass, and Mass. Thus, “SDn 1” would divide a tally bin by unity instead of the area (for an F2 tally) effectively multiplying the tally by the area. The new units would be particles, particles, MeV, and MeV rather than 1/cm2, 1/cm2, MeV/g, and MeV/g. The number of SD entries is the number of cells or surfaces on the F card times the number of segment bins.

2.4.3.7 D: Direct and SF, CF Flagging Bins

The D bins are described along with point detector tallies in Sect. 2.4.6. They are also used for surface, SF, and cell, CF, flagging. For example,

SF4 11.3 12

Tally 4 then has two D bins—the total, or usual tally, and that part of the tally that was made after particle tracks crossed surfaces 11.3 and 12.

For example,

CF1 20 21

Tally 1 then has two D bins—the total, or usual tally, and that part of the tally that was made after particle tracks entered cells 20 or 21.

2.4.3.8 U: User Bins

User bins, FU, enable input into tallies programmed by users who are altering the MCNP source code. They are more commonly used in conjunction with special tally treatments described in Sect. 2.4.4.

2.4.3.9 S: Segment Bins

Tally segmenting is shown in Example 2.19. The tally portion of that example is repeated here:

c ***** Tallies ************** FC901 Segmented surface tally F901:n (4.1 4.2) (11004.1 11004.2) (2.1<21<14) (2.1<22<11) (11004.1 4I 11004.6) (2.1<21[0 0 0]<14) FS901 -52 -53 T $ 4 segment bins FQ901 F S $ Print cells down, segments across FC906 Segmented energy deposition tally +F906 (1<21<14) (1<22<11) (21<14) (22<11) (1<21[-3:3 0:0 -3:3]<14) (1<22[-5:5 -4:4 0:0]<11) FS906 -51 -52 -53 T $ 5 segment bins FQ906 F S $ Print cells down, segments across SD906 1 759R $ MeV units, 152 cells x 5 segments

Note that surfaces 51, 52, and 53

51 PX 0 52 PX 8 53 PX 12

are used only for segmenting; they are not part of the problem geometry because these surfaces do not appear within the cell card section of the input. They only appear on FS tally segmenting cards.

FS901 -52 -53 T $ 4 segment bins

FS901 causes the surfaces of tally 901 to tally in 4 segment bins:

  • Segment bin 1: −52: all surface crossings with x < 8 cm, namely the square lattice;

  • Segment bin 2: −53: all surfaces not in the first bin but crossing with 8 < x < 12, namely the geometry between the square lattice and the hexagonal lattice;

  • Segment bin 3: everything else bin, which would include the hexagon lattice;

  • Segment bin 4: total over all segments, “whole,” which is the same tally that F901 would give without FS901 segment bins.

The F901 tally output is

1tally 901 nps = 100000 + Segmented surface tally tally type 1 number of particles crossing a surface. particle(s): neutrons surface b is (11004.1 11004.2) surface e is (11004.1 11004.2 1.2 1.3 4.5 4.6) surface f is (2.1<21[0 0 0]<14) segment 1: -52 segment 2: 52 -53 segment 3: 52 53 segment 4: whole surface segment: 1 2 3 whole surface (4.1 4.2) 2.55726E-01 0.0080 0.00000E+00 0.0000 0.00000E+00 0.0000 2.55726E-01 0.0080 b 0.00000E+00 0.0000 0.00000E+00 0.0000 1.79995E-01 0.0087 1.79995E-01 0.0087 (2.1<21<14) 3.62434E+00 0.0074 0.00000E+00 0.0000 0.00000E+00 0.0000 3.62434E+00 0.0074 (2.1<22<11) 0.00000E+00 0.0000 0.00000E+00 0.0000 2.09843E+00 0.0082 2.09843E+00 0.0082 e 5.02658E-01 0.0061 0.00000E+00 0.0000 4.58745E-01 0.0052 9.61403E-01 0.0038 f 1.08731E-01 0.0184 0.00000E+00 0.0000 0.00000E+00 0.0000 1.08731E-01 0.0184

The 1st, 3rd, and 5th surface bins of F901 are in the square lattice and thus score in segment bin 1. The 2nd, 4th, and 6th surface bins of F901 are in the hexagonal lattice and thus score in segment bin 3. The surface tallies FS segment bins are applied at the u=0 “real world” universe level coordinate system, not the lower universe level coordinate system of the surfaces being tallied.

FS906 -51 -52 -53 T $ 5 segment bins

FS906 causes the cells of tally 906 to tally energy deposition in five segment bins:

  • Segment bin 1: −51: all cell energy deposition with x < 0 cm, namely the −x portion of the square lattice;

  • Segment bin 2: −52: all cell energy deposition not in the first bin but with 0 < x < 8, namely the +x portion of the square lattice;

  • Segment bin 3: −53: all cell energy deposition not in the first two bins but crossing with 8 < x < 12, namely the geometry between the square lattice and the hexagonal lattice;

  • Segment bin 4: everything else bin, which would include the hexagon lattice;

  • Segment bin 5: total over all segments, “whole,” which is the same tally that F906 would give without FS906 segment bins.

The F906 tally output has 760 bins, so only parts are shown here:

1tally 906 nps = 100000 + Segmented energy deposition tally tally type 6+ energy deposition units mev/gram particle(s): neutrons cell e is (1<21[-3 0 -3]<14) cell f is (1<21[-2 0 -3]<14) cell g is (1<21[-1 0 -3]<14) cell eo is (1<22[-2 4 0]<11) cell ep is (1<22[-1 4 0]<11) cell eq is (1<22[0 4 0]<11) cell er is (1<22[1 4 0]<11) segment: 1 2 3 4 whole cell (1<21<14) 6.01161E+00 0.0113 6.85838E+00 0.0102 0.00000E+00 0.0000 0.00000E+00 0.0000 1.28700E+01 0.0101 (1<22<11) 5.39556E+00 0.0105 3.81918E+00 0.0131 0.00000E+00 0.0000 0.00000E+00 0.0000 9.21474E+00 0.0107 (21<14) 6.70886E+00 0.0106 7.93863E+00 0.0093 0.00000E+00 0.0000 0.00000E+00 0.0000 1.46475E+01 0.0093 (22<11) 6.28464E+00 0.0099 4.47032E+00 0.0119 0.00000E+00 0.0000 0.00000E+00 0.0000 1.07550E+01 0.0100 e 2.54071E-02 0.1040 3.18441E-02 0.0944 0.00000E+00 0.0000 0.00000E+00 0.0000 5.72512E-02 0.0869 f 5.09948E-02 0.0925 5.45180E-02 0.0800 0.00000E+00 0.0000 0.00000E+00 0.0000 1.05513E-01 0.0753 g 6.50667E-02 0.0749 6.09488E-02 0.0761 0.00000E+00 0.0000 0.00000E+00 0.0000 1.26016E-01 0.0678 eo 1.56021E-02 0.1168 1.42059E-02 0.1184 0.00000E+00 0.0000 0.00000E+00 0.0000 2.98080E-02 0.1034 ep 1.08991E-02 0.1180 1.17894E-02 0.1248 0.00000E+00 0.0000 0.00000E+00 0.0000 2.26885E-02 0.1062 eq 9.31922E-03 0.1227 9.11473E-03 0.1529 0.00000E+00 0.0000 0.00000E+00 0.0000 1.84340E-02 0.1178 er 4.57305E-03 0.1316 5.53235E-03 0.1952 0.00000E+00 0.0000 0.00000E+00 0.0000 1.01054E-02 0.1476

All energy deposition is scored in the 1st and 2nd segment bins, namely x < 0 and 0 < x < 8. The hexagonal lattice cells that fill cell 11 with x > 12 do not score in the 4th segment bin with x > 12 because cell segmenting is performed at the universe level of the cell. Both the square lattice that fills cell 14 and the hexagonal lattice that fills cell 11 have their centers at the origin in the u = 1 lowest level universe. Therefore, cell segments apply to the coordinate system of the level of the cell being tallied.

The segments are printed horizontally and the surfaces and cells vertically because of the FQ901 and FQ906 tally format specifications.

SD906 has 760 entries: 152 cells × 5 segments. The value of 1 replaces the normal FD6 divisor of grams, changing the units from MeV/g to MeV.

2.4.4 Special Tally Treatments

Special tally treatments may be applied to various tallies:

Form: FTn id1 p1,1 p1,2 … id2 p2,1 p2,2… … n = tally number Id = Special tally treatments given below pi,j = parameter j for the ith tally treatment.

The special tally treatments are

  • FRV—Fixed arbitrary reference direction for tally 1 cosine binning

  • GEB—Gaussian energy broadening

  • TMC—Time convolution

  • INC—Identify the number of collisions

  • ICD—Identify the cell from which each detector score is made

  • SCX—Identify the sampled index of a specified source distribution

  • SCD—Identify which of the specified source distributions was used

  • ELC—Electron current tally

  • PTT—Put different multigroup particle types in different user bins

  • RES—Residual nuclei

  • TAG—Tally tagging

  • LET—Tally stopping powers instead of energy

  • ROC—Receiver-operator characterization

  • PDS—Point detector sampling

  • FFT—First fission tally

  • COM—Compton image tally

  • CAP—Coincidence capture

  • PHL—Pulse-height light tally

2.4.4.1 FRV: Fixed Arbitrary Reference Direction for Tally F1 or F2 Cosine Binning

FTn FRV u v w

The default cosine bins are relative to the surface normal. The FRV treatment makes the cosine binning relative to the vector (u,v,w).

2.4.4.2 GEB: Gaussian Energy Broadening

The tallying of a gamma-ray pulse height spectrum in the MCNP code, by default, is equivalent to using a detector with infinite pulse height resolution. This is a very useful feature for studying the properties of the incoming gamma-rays. For simulation of experiments it may be necessary to include the actual resolution of the detector, for example, in triggering coincidences (see Sect. 2.4.4.18). Also, a user may simply want the simulated spectrum to resemble an experimentally measured spectrum. In these cases the MCNP tally can be broadened with the built-in tally treatment, GEB, to simulate the finite resolution of a real detector (an example is given in Sect. 3.3.4). The variation of resolution of the detector as a function of energy is accommodated by three parameters. The usage is as follows:

FTn GEB a b c

The parameters specify the full width at half maximum (FWHM) of the observed energy broadening in a physical radiation detector: \( \textrm{FWHM}(E)=a+b\sqrt{E+c{E}^2} \), where E is the energy of the particle. The units of a, b, and c are MeV, MeV1/2, and 1/MeV, respectively. The energy actually scored is sampled from a Gaussian distribution with that FWHM (from MCNP manual Sect. 3.3.5.18).

2.4.4.3 TMC: Time Convolution

FTn TMC a b

The tally scores are made as if the source were a square pulse starting at time a < time < b. All particles should be started at time zero in the source definition (SDEF).

2.4.4.4 INC: Identify the Number of Collisions

FTn INC FUn 0 1 8I 10 1e8 T UNC:N k

The FUn bins are the number of collisions that have occurred in the track. The first bin has uncollided particles; the second has those with exactly one collision; the third entry interpolates so that those particles having exactly 2, 3, 4, 5, 6, 7, 8, 9, 10 collisions are in those bins. The last entry, 1e8, should be a large number to count all collisions. The UNC card controls whether the collision count starts over for each secondary particle (k = 1, default) or continues (k = 0).

2.4.4.5 ICD: Identify the Cell from Which Each F5 Detector Score Is Made

FTn ICD FUn 11 12 13 … 88 89 T

The FUn bins are the names of some or all of the cells in the problem. Thus, the detector tally is subdivided into bins according to which cell had the source or collision that resulted in the detector score.

2.4.4.6 SCX: Identify the Sampled Index of a Specified Source Distribution

FTn SCX k

One user bin is created for each bin of source distribution SIk plus a total bin, thus subdividing the tally into parts corresponding to each bin of the source distribution.

2.4.4.7 SCD: Identify Which of the Specified Source Distributions Was Used

FTn SCD FUn k1 k2 k3 . . .

Each FUn user bin of the tally is the portion of the tally from source distribution ki. If more than one of the source distributions listed on the FU card is used for a given source particle, only the first one used will score.

2.4.4.8 ELC: Electron Current Tally

FTn ELC c

ELC specifies how the charge of a particle is scored to an F1 tally. Without the ELC treatment, the F1 tally gives particle current without regard for the particle charge. The possible values for c are

  • c = 1 to cause negatively charged particles to make negative scores,

  • c = 2 to put charged particles and antiparticles into separate user bins, and

  • c = 3 for the effect of both c = 1 and c = 2.

If c = 2 or 3, three user bins (e.g., positrons, electrons, and total) are created.

2.4.4.9 PTT: Put Different Multigroup Particle Types in Different User Bins

FTn PTT FUn a1 a2 a3 . . .

Each FUn user bin is the atomic weights in units of MeV of particles that masquerade as neutrons in a multigroup data library (MGOPT multigroup option turned on). Thus a1,a2 = 0.511, 0 separates the tally into electrons and photons.

2.4.4.10 RES: Residual Nuclei

Fn:# a1 a2 a3 . . . FTn RES 8016 20040 26000 92238

For heavy ion (PAR = #) surface current (F1), surface flux (F2), cell flux (F4), and cell heating (F6) tallies, n = 1, 2, 4, 6, the tally is subdivided into user bins according to which nuclide heavy ion, 13O, 40Ca, all Fe isotopes, and 238U made the tally in each cell or surface a1,a2,a3.

F8:# a1 a2 a3 FT8 RES z1 z2

For F8 tallies, the weight of all ions with Z numbers z1 < z < z2 produced by any reaction with any particle in the problem are recorded. If z1, z2 are ZAID identifiers, then the production of those nuclides is tallied.

MODE N H F128:# 1101 1102 FT128 RES 5 8 FQ128 U F

This F128 tally type 8 would list the production of all possible isotopes of B, C, N, and O. The production would come from primary reactions with either neutrons (N) or protons (H) and not secondary reactions if the produced nuclide decays (requiring the ACT activation card). The FQ128 tally format card lists the cells (F) horizontally and residual nuclei (U) vertically.

F258:# 1103 9I 1113 FT258 RES 1002 1003 2003 2004 FQ258 F U

This F258 tally type 8 lists the production of deuterons, tritons, helions, and alphas from all problem particle types in cells 1103–1113. The FQ258 format lists the four light ions horizontally and the eleven cells vertically in the output.

2.4.4.11 TAG: Tally Tagging

FTn TAG a FUn b1 b2 b3 b4 . . .

Tally tagging separates a tally into bins by how and where the scoring particle was produced and works for all but type 8 tallies. The TAG options are

  • a=1: collided particles lose their tag; bremsstrahlung and annihilation photons are included in the scattered FU 0 bin;

  • a=2: collided particles lose their tag; bremsstrahlung and annihilation photons are given special tags for segregation;

  • a=3: all collided particles retain their production tag (most useful option); and

  • a=4: all collided particles except Compton photoatomic retain their production tag.

Each FUn bin, bi, has the form CCCCCZZAAA.RRRRR

  • CCCCC = cell number or 00000 to tag all cells

  • ZZAAA = target nuclide identifier

  • RRRRR = reaction identifier (e.g., 00102 for n,γ ) or residual nuclide ZAID

Tally tagging FUn bin, bi, special cases:

  • −0000000001 or −1 source particle tag for all cells

  • −CCCCC00001 source (i.e., uncollided) particle tag for cell CCCCC

  • 0000000000 or 0 scattered particle tag

  • 10000000000 or 1e10 everything else tag

Tally tagging FUn bin, bi, special cases for neutrons:

  • ZZAAA.99999 delayed particles from fission or residuals of ZZAAA

Tally tagging FUn bin, bi, special cases for photons:

00000.00001

bremsstrahlung from electrons

ZZ000.00003

fluorescence from nuclide ZZ

00000.00003

K X-rays from electrons

00000.00004

annihilation photons from e−

ZZ000.00005

Compton photons from nuclide ZZ

ZZAAA.00006

muonic X-rays from nuclide ZZAAA

00000.00007

Cerenkov photons

Tally tagging FUn bin, bi, special cases for electrons:

Electron special designations for ZZAAA.RRRRR:

ZZ000.00001

photoelectric from nuclide ZZ

ZZ000.00003

Compton recoil from nuclide ZZ

ZZ000.00004

pair production from nuclide ZZ

ZZ000.00005

Auger electron from nuclide ZZ

00000.00005

Auger electron from electrons

00000.00006

knock-on electrons

Tally tagging photon example:

FT5 TAG 3 FU5 -1.0 $ Source from photons 0000106012.00005 $ Compton from carbon in cell 1 0000106012.00000 $ Remaining photons from carbon in cell 1 0000026056.00102 $ Capture gammas from 56Fe in cell 1 0000026056.00000 $ Remaining photons / gammas from 56Fe 0000000000.00051 $ Remaining 1st level from (n,n’) gammas 10000000000.00000 $ Remaining gammas

2.4.4.12 LET: Tally Stopping Powers Instead of Energy

FTn LET

Linear energy transfer values provided in the energy bins of charged particle tallies are interpreted as stopping power values with units of MeV/cm.

2.4.4.13 ROC: Receiver-Operator Characterization

FTn ROC n TFn F D U S M C E T F D U S M C E T

The FT ROC receiver-operator separates the tally into signal and noise and runs histories in batches of n source particles, where the total number of NPS source particles should be 50–100 times n. The 1st eight entries on the TFn tally fluctuation bin card point to the signal bin of the tally; the 2nd eight entries point to the noise bin. For example, the signal may be the value in one energy, time, cosine, or other bin, and the noise may be the value in some other bin. The probability of detection and the probability of false alarm are then printed in Table 163 of the MCNP output file.

2.4.4.14 PDS: Point Detector Sampling

FTn PDS c

Point detectors (n ending in 5, or F5 tallies) suffer convergence problems when collision sampling has a large variation in probabilities. For example, photon coherent scattering in the forward direction has orders of magnitude greater tally contribution than in the backward direction. In these cases—particularly with photon point detectors, F5:p—the PDS special treatment changes sampling to more often sample the rare high contribution events.

c = 0 or c = −1 Next-event estimator sampling is performed post-collision; only a single reaction and isotope is sampled (historic MCNP4 and MCNP5 behavior) (DEFAULT).

c = 1 Next-event estimator sampling is performed using post-collision sampling of the collision isotope and pre-collision sampling of all reaction channels (recommended for photons).

c = 2 Next-event estimator sampling is performed using pre-collision sampling of all collision isotopes and pre-collision sampling of all reaction channels.

2.4.4.15 FFT: First Fission Tally

FTn FFT [LKJI] FUn z1 z2 z3 . . .

The first fission tally records the nuclide with the first fission or (n,xn) reaction. This special treatment was developed specifically for safeguards applications because the signal from fissions is roughly scaled by the number of fission neutrons emitted by the first fissioning nuclide. The optional LKJI values are

  • L = 0/1 Omit/include neutron and photon-induced fissions treated by model physics

  • K = 0/1 Omit/include neutron spontaneous fissions (PAR = SF source particles)

  • J = 0/1 Omit/include photon-induced fissions treated by library physics

  • I = 0/1 Omit/include neutron-induced fissions treated by library physics (E <~20 MeV)

The zi values are ZAID identifiers, 0 for everything else, 16 for (n,xn), and 18 for some other fission. For example,

FTn FFT 0001 FUn 92238 0 16 18 94241 92235 94239

The FFT 0001 is the default, including only fissions treated by library physics. The tally is subdivided into user bins for those particles having had their first multiplicity reaction from 235U fission (92235), 238U fission (92238), 239Pu fission (94239), 241Pu fission (94241), any other fission (18), (n,xn) (16), or no multiplicity reactions at all before scoring.

Note that the FFT special treatment causes cell flagging (CF) and tally segmenting (FS) to apply at the location where the first fission or (n,xn) event occurs and not where the tally score is made. The FFT tally does not apply to F8 tallies.

2.4.4.16 COM: Compton Image Tally

FT8 COM T A

The FT8 COM tally option produces a Compton image stored in an associated FIR radiography tally T using algorithm A (only A = 1 is currently available). The Compton image is formed from an FT8 PHL specification of dual-region coincidences of planar lattice tallies. At the end of each particle history, Compton/photoelectric energy deposition in the front/back of these dual-panel detectors is used to create a circular “image” of the incident photon on a specified image plane. For more details, see the MCNP code manual section 3.3.5.18.

2.4.4.17 CAP: Coincidence Capture

The coincidence capture tally is a key safeguards application that models multiplicity counting methods. The tally produces the factorial moments of the neutron multiplicity distribution. The use of this tally is described in Sect. 3.2. The FT8 CAP syntax is

$$ \textrm{FT}8\kern0.5em \textrm{CAP}\kern0.5em \left[-{m}_{\textrm{c}}\right]\left[-{m}_{\textrm{o}}\right]\kern0.5em {i}_1{i}_2\dots {i}_n\kern0.5em \left[\textrm{GATE}\kern0.5em {t}_{\textrm{d}}\kern0.5em {t}_{\textrm{w}}\right]\kern0.5em \left[\textrm{EDEP}\kern0.5em {t}_{\textrm{g}}\kern0.5em {t}_{\textrm{t}}\right] $$

where

  • mc = optional maximum number of captures (Default = 21),

  • mo = optional maximum number of moments (Default = 12),

  • in = capture nuclides such as 3006 or 5010 for 6Li or 10B,

  • td = predelay time (shakes),

  • tw = gate width (shakes),

  • tg = trigger tally number, and

  • tt = trigger tally threshold (MeV) (Default=0.0).

Coincidence and multiplicity calculations are described in greater detail in Sect. 3.1. The FT8 tallies from Example 3.2 are

fc8 Ungated Coincidence tally f8:n (10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27) ft8 cap 2003 c fc18 Gated Coincidence tally f18:n (10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27) ft18 cap 2003 gate 450 6400 $ 1us=10 shakes --> 4.5 us predelay and 64 us gate c FC108 Full gate to compare moments with ungated F8 F108:n (10 16I 27) FT108 CAP 2003 gate 0 1e20 c FC116 Energy deposition tally to demonstrate FT118 CAP EDEP F116:n (10 16I 27) c FC118 Energy deposition capture tally F118:n (10 16I 27) FT118 CAP EDEP 116 .5

The FT8 CAP tally counts the number of captures and records the multiplicity and calculates the moments. The data are stored in cosine bins that enable tally plotting because there are no multiplicity bins in the 8-dimensional tally array. In the standard tally output (in the cosine format), these data are hard to read, so they are also presented in more readable form in Print Table 118:

1 neutron captures, moments and multiplicity distributions. tally 8 print table 118 weight normalization by source fission neutrons = 2154543 cell a is (10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27) cell: a neutron captures on 3he captures captures multiplicity fractions histories by number by weight by number by weight error captures = 0 649116 0 0.00000E+00 6.49116E-01 3.01278E-01 0.0007 captures = 1 280151 280151 1.30028E-01 2.80151E-01 1.30028E-01 0.0016 captures = 2 58825 117650 5.46055E-02 5.88250E-02 2.73028E-02 0.0040 captures = 3 9573 28719 1.33295E-02 9.57300E-03 4.44317E-03 0.0102 captures = 4 1811 7244 3.36220E-03 1.81100E-03 8.40549E-04 0.0235 captures = 5 380 1900 8.81858E-04 3.80000E-04 1.76372E-04 0.0513 captures = 6 108 648 3.00760E-04 1.08000E-04 5.01266E-05 0.0962 captures = 7 29 203 9.42195E-05 2.90000E-05 1.34599E-05 0.1857 captures = 8 5 40 1.85654E-05 5.00000E-06 2.32068E-06 0.4472 captures = 9 2 18 8.35444E-06 2.00000E-06 9.28271E-07 0.7071 total 1000000 436573 2.02629E-01 1.00000E+00 4.64136E-01 0.0015 factorial moments by number by weight 3he 4.36573E-01 0.0015 2.02629E-01 0.0015 3he(3he-1)/2! 1.04651E-01 0.0051 4.85722E-02 0.0051 3he(3he-1)(3he-2)/3! 2.42400E-02 0.0177 1.12506E-02 0.0177 3he(3he-1) .... (3he-3)/4! 6.94800E-03 0.0514 3.22481E-03 0.0514 3he(3he-1) .... (3he-4)/5! 2.16900E-03 0.1171 1.00671E-03 0.1171 3he(3he-1) .... (3he-5)/6! 6.19000E-04 0.2259 2.87300E-04 0.2259 3he(3he-1) .... (3he-6)/7! 1.41000E-04 0.3846 6.54431E-05 0.3846 3he(3he-1) .... (3he-7)/8! 2.30000E-05 0.5619 1.06751E-05 0.5619 3he(3he-1) .... (3he-8)/9! 2.00000E-06 0.7071 9.28271E-07 0.7071

These results are slightly different (though statistically the same) from those in Fig. 3.11, which is part of the more detailed discussion of coincidence counting in Sect. 3.2.6.2. The difference arises from running the problem on different computer platforms and operating systems. Tally F8 is the simplest form, counting all captures in 3He in the cells listed on F8:n.

1 neutron captures, moments and multiplicity distributions. tally 18 print table 118 weight normalization by source fission neutrons = 2154543 cell a is (10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27) cell: a neutron captures on 3he time gate: predelay = 4.5000E+02 gate width = 6.4000E+03 pulses occurrences occurrences pulse fraction in gate histogram by number by weight by number by weight error captures = 0 372500 0 0.00000E+00 3.72500E-01 1.72890E-01 0.0014 captures = 1 55609 55609 2.58101E-02 5.56090E-02 2.58101E-02 0.0044 captures = 2 7100 14200 6.59072E-03 7.10000E-03 3.29536E-03 0.0129 captures = 3 1134 3402 1.57899E-03 1.13400E-03 5.26330E-04 0.0332 captures = 4 187 748 3.47173E-04 1.87000E-04 8.67933E-05 0.0827 captures = 5 36 180 8.35444E-05 3.60000E-05 1.67089E-05 0.1924 captures = 6 7 42 1.94937E-05 7.00000E-06 3.24895E-06 0.5533 total 436573 74181 3.44300E-02 4.36573E-01 2.02629E-01 0.0045 factorial moments by number by weight n 7.41810E-02 0.0055 3.44300E-02 0.0054 n(n-1)/2! 1.20890E-02 0.0201 5.61093E-03 0.0201 n(n-1)(n-2)/3! 2.38200E-03 0.0612 1.10557E-03 0.0612 n(n-1)(n-2) ... (n-3)/4! 4.72000E-04 0.1549 2.19072E-04 0.1549 n(n-1)(n-2) ... (n-4)/5! 7.80000E-05 0.3140 3.62026E-05 0.3140 n(n-1)(n-2) ... (n-5)/6! 7.00000E-06 0.5533 3.24895E-06 0.5533

Tally F18 is a gated capture tally. When a neutron is captured, a “gate” is opened. After a predelay of 450 shakes (4.5 μs), the number of captures is counted for the next 6400 shakes (64 μs).

1 neutron captures, moments and multiplicity distributions. tally 108 print table 118 weight normalization by source fission neutrons = 2154543 cell a is (10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27) cell: a neutron captures on 3he time gate: predelay = 0.0000E+00 gate width = 1.0000E+20 pulses occurrences occurrences pulse fraction in gate histogram by number by weight by number by weight error captures = 0 350884 0 0.00000E+00 3.50884E-01 1.62858E-01 0.0014 captures = 1 70733 70733 3.28297E-02 7.07330E-02 3.28297E-02 0.0036 captures = 2 11908 23816 1.10539E-02 1.19080E-02 5.52693E-03 0.0091 captures = 3 2335 7005 3.25127E-03 2.33500E-03 1.08376E-03 0.0207 captures = 4 524 2096 9.72828E-04 5.24000E-04 2.43207E-04 0.0437 captures = 5 144 720 3.34178E-04 1.44000E-04 6.68355E-05 0.0833 captures = 6 36 216 1.00253E-04 3.60000E-05 1.67089E-05 0.1667 captures = 7 7 49 2.27426E-05 7.00000E-06 3.24895E-06 0.3780 captures = 8 2 16 7.42617E-06 2.00000E-06 9.28271E-07 0.7071 total 436573 104651 4.85722E-02 4.36573E-01 2.02629E-01 0.0036 factorial moments by number by weight n 1.04651E-01 0.0051 4.85722E-02 0.0049 n(n-1)/2! 2.42400E-02 0.0177 1.12506E-02 0.0177 n(n-1)(n-2)/3! 6.94800E-03 0.0514 3.22481E-03 0.0513 n(n-1)(n-2) ... (n-3)/4! 2.16900E-03 0.1171 1.00671E-03 0.1171 n(n-1)(n-2) ... (n-4)/5! 6.19000E-04 0.2259 2.87300E-04 0.2259 n(n-1)(n-2) ... (n-5)/6! 1.41000E-04 0.3846 6.54431E-05 0.3846 n(n-1)(n-2) ... (n-6)/7! 2.30000E-05 0.5619 1.06751E-05 0.5619 n(n-1)(n-2) ... (n-7)/8! 2.00000E-06 0.7071 9.28271E-07 0.7071

If the gate is infinite, as in F108, then the 1st, 2nd, 3rd, … moments are the same as the ungated 2nd, 3rd, 4th, … moments of the F8 tally.

1 neutron captures, moments and multiplicity distributions. tally 118 print table 118 weight normalization by source fission neutrons = 2154543 cell a is (10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27) cell: a Energy losses > 5.0000E-01 MeV in tally F116 captures captures multiplicity fractions histories by number by weight by number by weight error captures = 0 791900 0 0.00000E+00 7.91900E-01 3.67549E-01 0.0005 captures = 1 183826 183826 8.53202E-02 1.83826E-01 8.53202E-02 0.0021 captures = 2 21834 43668 2.02679E-02 2.18340E-02 1.01339E-02 0.0067 captures = 3 2138 6414 2.97697E-03 2.13800E-03 9.92322E-04 0.0216 captures = 4 267 1068 4.95697E-04 2.67000E-04 1.23924E-04 0.0612 captures = 5 31 155 7.19410E-05 3.10000E-05 1.43882E-05 0.1796 captures = 6 3 18 8.35444E-06 3.00000E-06 1.39241E-06 0.5773 captures = 7 1 7 3.24895E-06 1.00000E-06 4.64136E-07 1.0000 total 1000000 235156 1.09144E-01 1.00000E+00 4.64136E-01 0.0021 factorial moments by number by weight n 2.35156E-01 0.0021 1.09144E-01 0.0021 n(n-1)/2! 3.02260E-02 0.0077 1.40290E-02 0.0077 n(n-1)(n-2)/3! 3.61100E-03 0.0302 1.67599E-03 0.0302 n(n-1)(n-2) ... (n-3)/4! 5.02000E-04 0.1080 2.32996E-04 0.1080 n(n-1)(n-2) ... (n-4)/5! 7.00000E-05 0.3440 3.24895E-05 0.3440 n(n-1)(n-2) ... (n-5)/6! 1.00000E-05 0.7211 4.64136E-06 0.7211 n(n-1)(n-2) ... (n-6)/7! 1.00000E-06 1.0000 4.64136E-07 1.0000

The F8 CAP tally scores when a capture occurs. This is a good approximation for detectors like 3He, in which almost 100% of captures lead to a detected pulse. The tally does, however, have the ability to count “captures” by energy deposition, which is useful for simulation of other detectors such as boron-lined counters. For example,

FT118 CAP EDEP 116 .5

counts a capture for every track in the specified cells that deposit greater than 0.5 MeV of energy in cell 116. In the case of a 3He detector, a closer agreement with the F8 tally can be achieved by using the pulse-height energy deposition tally with a lower threshold to count all of the events:

FC128 *F8 energy deposition tally to demonstrate FT118 CAP EDEP *F128:n (10 16I 27) E128 100 NT c FC118 Energy deposition capture tally F118:n (10 16I 27) FT118 CAP EDEP 128 1e-10

The cells specified on the F118:n card are required but ignored; the tally is based on the energy deposition in tally F128. The *F8 heating tally is a surface estimator that tallies the difference between particles that enter and exit the cells specified by *F128. Note that a threshold of zero will not work because MCNP deducts 10−12 of the energy that passes through each cell in order to distinguish between particles that never entered the cell and those that entered but did not lose any energy. An example of the use of the EDEP option is given in Sect. 4.4.2.

2.4.4.18 PHL: Pulse-Height Light Tally

The pulse-height light special treatment, FT8, is described with the other F8 pulse-height tallies in Sect. 2.4.5.

2.4.5 Pulse-Height Tallies

Safeguards detectors do not measure flux or energy deposition directly. They measure pulses or signals. The MCNP pulse-height tally and its various forms are used to model scintillator detectors.

The F8 pulse height looks like just another cell tally

f8:n c1 c2 c3 . . .

but it is unlike other MCNP tallies in that it does not model currents and fluxes or reaction rates that can be modeled with the Boltzmann transport equation; rather it is a non-linear tally. With the FT8 RES and FT8 CAP special treatments, it is a collision estimator that counts residuals and captures. In its simplest form, it is a surface tally, counting the energy that enters a cell and subtracting out what leaves to put a pulse of that energy in the tally energy bins.

The MCNP code’s F1, F2, F4, F5, F6, and F7 tallies score the energy bins as the energy at which the event happened, not the amount deposited. Consequently, the E0 default energy can apply to either those standard tallies or pulse-height tallies, but not both. Charged particles, on average, emit the correct number of secondary particles. However, any given track may sample the production of more secondary particles that carry more energy than the incident particle which leads to a negative pulse-height energy deposition. Negative energy production from sampling too many secondary particles is nonphysical but correct, on average, to get the total energy deposited. Thus, pulse-height tallies should have a zero energy bin to tally these rare nonphysical imbalance events with negative scores. Particles that do not enter the pulse-height tally cell make no score. Particles that enter the pulse-height tally cell but deposit no energy need to be distinguished, so the MCNP code gives them an energy deposition of 10−12 times the particle’s energy. Thus the energy bins for a pulse-height tally should be something like

e8 0 1e-6 .001 .01 .1 100ilog 101

The energy pulse is accumulated from all tracks of a particle’s history, including those of its created secondary particles. Consequently, variance reduction—which affects the particle weight—greatly complicates the accumulation of these tracks. Variance reduction can be used with pulse-height tallies in some cases, but the algorithm is complicated and fragile. Thus, variance reduction should be avoided. If used, it should be done with caution [7].

A *F8 tally multiplies the pulse by the energy of the track at each collision in the pulse-height cell.

Consider a photon that enters a cell at 10 MeV, then exits at 9 MeV, then re-enters at 7 MeV, and exits at 6 MeV. The F6 tally would score 1 MeV in the 10 MeV bin and 1 MeV in the 7 MeV bin.

The *F8 tally would have 2 MeV in the 2 MeV bin. The F8 tally would have 1 pulse in the 2 MeV bin.

Pulse-height tallies may also be used for the FT8 PHL pulse-height light special treatment. The PHL treatment enables more realistic detector modeling including coincidence and anticoincidence. Example 2.20 has simple F8, energy *F8, and light FT8 PHL pulse-height tallies:

Example 2.20 NE213 Scintillator

NE213 scintillator 11 1 -0.9 -11 imp:n=1 20 0 -100 11 imp:n=1 22 0 100 imp:n=0 11 RCC 0 0 0 0 0 3.81 2.54 100 SO 100 sdef pos 0 0 -1 erg 15 vec 0 0 1 dir 1 m1 1001 1.213 6000 1 nps 1e5 print -85 -86 -162 -70 -98 prdmp 2j -1 mode n h phys:n J 20 J J J J 1 $ analog < 20 MeV, ion elastic recoil cut:n 2J 0 0 cut:h j 0 0 0 FC6 Track length heating converted to light output f6:h 11 $ F6 is converted to light DE6 LIN 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 9 9.2 9.4 9.6 9.8 10 DF6 LIN 0.00 0.01 0.04 0.08 0.12 0.17 0.23 0.29 0.35 0.42 0.50 0.57 0.66 0.74 0.83 0.93 1.02 1.12 1.23 1.34 1.45 1.56 1.67 1.79 1.92 2.04 2.17 2.30 2.43 2.57 2.70 2.84 2.99 3.13 3.28 3.43 3.58 3.74 3.90 4.05 4.22 4.38 4.55 4.72 4.89 5.06 5.23 5.41 5.59 5.77 5.95 FC16 Track length MeV/g deposited f16:h 11 $ F16 is Mev/g deposited FC26 Track length MeV deposited f26:h 11 sd26 1 $ Divide by 1 rather than grams e0 0.2 48I 10 11 13 15 17 20 FC8 Total light output f8:h 11 $ F8 is total light output FT8 PHL 1 6 1 0 e8 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50 FC18 Total energy deposited F18:h 11 Ft18 PHL 1 16 1 0 $ 0.37 1.32 and 3.79 are expt thresholds e18 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50 FC28 Pulse height tally f28:h 11 e28 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50 FC38 Pulse height energy deposition *f38:h 11 e38 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50

Example 2.20 is an oversimplified NE-213 scintillator model that consists of a 3.81 cm high and 2.54 cm radius can with a mono-directional, 15 MeV neutron source directed up the z-axis of the can. This example is a coupled neutron-proton calculation: MODE N H. The 2nd PHYS:N entry turns on analog capture below 20 MeV. Analog capture is more evidently specified by setting the weight cutoffs, the 3rd and 4th entries of CUT: N 2J 0 0 to zero. The 7th PHYS:N entry turns on proton recoil when a neutron scatters off hydrogen, and the scattered hydrogen nucleus is transported as a proton. CUT:H J 0 0 0 turns on analog capture, and the 2nd entry lowers the proton energy cutoff to zero. The default for protons is a 1 MeV lower cutoff, but now the MCNP code will set it to 0.001 MeV (1 keV), the lowest possible for protons in the MCNP code.

The F16:H tally is a standard, simple, proton-heating or energy-deposition tally with units of MeV/g. The F26:H tally has a segment divisor, DS26, which overrides the divisor of grams, and thus it is the standard, simple, proton-heating or energy-deposition tally with units of MeV. The F6:H proton tally is converted to light output by the dose function, DE6 and DF6. The dose functions are obtained from experimental measurement as described in Chap. 8 of Ref. [8]. The tally is multiplied by the linearly interpolated point on DF6 at the linearly interpolated values between the energies on DE6. Because the tally is multiplied by as much as a factor of 6 at higher energies, the upper tally energy bin is exceeded, and the MCNP code issues the warning

warning. tally not scored beyond last energy bin. nps = 176 tal = 8 erg = 5.2820E+01

The energies of the F6 tallies are set by the E0 default energy. The pulse-height energy bins must be specified for each pulse-height tally because pulse-height tallies and non-pulse-height tallies cannot both use the same E0 default energy binning.

Tally F28 is a simple pulse-height tally with units of pulses in each energy bin.

Tally *F38 is an energy pulse-height tally with units of MeV in each energy bin. The total bin value is 0.909 MeV, which agrees well with the F26 total bin of 0.907 MeV. F38 uses a surface estimator and F26 uses a track-length estimator.

Tallies F8 and F18 are pulse-height light tallies. The FT8 PHL syntax is

FT8 PHL n ta1 ba1 ta2 ba2 ... m tb1 bb1 tb2 bb2 ...

  • n = number of F6 tallies for the first detector region

  • tai, bai = pairings of tally number and F-bin number for the n F6 tallies of the first detector region

  • m = number of F6 tallies for the second detector region

  • tbi, bbi = pairings of tally number and F-bin number for the m F6 tallies of the second detector region

The pulse-height light tally gets its energy deposition in the cell not from the surface estimator weight balance or a collision estimator at each collision but rather from another energy deposition tally. The energy deposition tally may be a type 6 tally, as illustrated here (F6), or a type 8 tally (*F38). Because these heating tallies can have different particle types, the FT8 PHL tally can mix and match contributions from different particles. For example,

mode n h d s t a phys:n 6J 4 cut:h,d,t,s,a J 0 0 0 e0 100 nt f116:h 3 4 f126:d 3 4 f136:t 3 4 f146:s 3 4 f156:a 3 4 fc208 Combined PHT for h,d,t,s,a f208:n 3 e208 0 1e-3 99ilog 10 100 ft208 phl 5 116 1 126 1 136 1 146 1 156 1 0 GEB .01 .02 .04

In this example, the 7th PHYS:N entry enables neutrons that strike any of the light nuclei to generate secondary light ions from both elastic and selected inelastic reactions. Analog capture and the lowest possible energy cutoff are invoked with the CUT entries. The pulse-height light tally, FT8 PHL has five heating tallies in its one region, each using the energy deposition in the first cell, cell 3, of the F6 tally.

In addition to the F8 pulse-height tally, the *F8 energy pulse-height tally, and the FT8 PHL pulse-height light tally, there is the possibility of a coincidence–anticoincidence pulse-height light tally by having up to four FT8 PHL regions. For example,

f116:h 3 4 fc308 PHT for h in cell 3 no/yes having cell 4 energy deposition f308:h (3 4) e308 0 1e-3 99ilog 10 100 fu308 0 100 fq308 e u ft308 phl 1 116 1 1 116 2 GEB .01 .02 .04

In this example, the FT308 PHL has 2 regions: tally F116 cell 3 (first cell) and tally F116 cell 4 (second cell). The first region, cell 3, uses energy bins E308. The second region, cell 4, uses user bins FU308 for energy bins. (A third region would use C cosine bins and a fourth region would use FS segment bins for energy.)

In this case (with two regions) the scores are made in a two-dimensional matrix, where one dimension represents the score in cell 4 (in user bin energies defined on the fu308 card) and the other dimension represents the score in cell 3 (in regular energy bins defined on e308 card). If there is no score in cell 4, the score recorded in cell 3 appears in the lowest energy bin of cell 4. If there is no score recorded in cell 3, the score in cell 4 appears in the lowest energy bin of cell 3. If there is no score in either cell there is no entry in the tally. The fq308 e u card causes the user bins (cell 4 score) to run horizontally and the regular energy bins (cell 3 score) to run vertically. The first row of data is the pulse height spectrum in cell 4 when there is no score in cell 3. The first column of data is the pulse height spectrum in cell 3 when there is no score in cell 4. The other bins show the energy relationship between coincident events.

Thus, coincidence/anticoincidence is achieved by seeing which parts of the pulse height were caused from contributions from which regions. As another example,

F8:N 5 FT8 PHL 1 6 1 1 16 1 E8 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 FU8 0 10 F6:E 5 F16:E 6

The FU8 0 10 card creates two user bins, one <0 and one 0-10. The first user bin of the pulse-height tally in cell 5 (F6 electron energy deposition) scores only when there is no energy deposition in cell 6 (F16); the second user bin records the pulse-height tally in cell 5 (from F6) only when there is energy deposition in cell 6 (F16) with 0 < E < 10.

More information on the use of pulse height tallies can be found in [9].

2.4.6 Point Detectors and Next-Event Estimators

Point detectors enable the estimation of flux at a point. The probability of a neutron or a photon having a random walk exactly to a point is zero, but the flux at the point can be estimated using a next-event estimator. The next-event estimator may also be used to deterministically transport neutrons or photons to a spherical region—the basis of the variance-reduction technique DXTRAN (deterministic transport).

The next-event estimator accumulates the flux, Φ, at a point by a contribution of

$$ \varPhi = Wp\left(\mu \right){e}^{-\sigma x}/2\pi {R}^2 $$

This estimate of Φ is accumulated at every source or collision event.

  • W = particle weight

  • p(μ) = p(Ω)/2 = scattering probability density function

  • R = distance between event and detector point

  • eσx = attenuation through all cells/materials along path R

Although the probability of scattering exactly in the direction, Ω, from the event to the detector is zero, the probability density function of the cosine, μ, assuming azimuthal symmetry, p(μ) = p(Ω)/2, is finite.

0 ≤ p(μ) < ∞ -1 ≤ μ ≤ 1

For photon coherent scattering with form factors, p(μ) is nearly a delta function at μ = 1.

p(1) ≈ ∞

Consequently, convergence is poor unless coherent scatter is turned off (PHYS:P 2J 1) or the special treatment, FT5 PDS, is applied to sample collision outcomes rather than incident reaction probabilities.

Convergence is also poor if there are large estimates from collisions close to the detector, in which case R approaches zero. The MCNP input for a point detector allows specification of a radius of exclusion, R0. If the detector is in a vacuum, then no collisions will occur close to the detector, and R0 can safely be set to zero. R0 can also be set to zero if there are no large estimates from collisions close to the detector as is often the case in air. If there are large collisions close to the detector, then R0 can be specified as about a mean free path; all collisions inside this radius assume a flat flux uniform distribution, which is a good approximation if R0 is small. In our experience a better approach is to use the nested DXTRAN spheres variance reduction technique to better sample collisions near the point detector.

The attenuation term, eσx, slows calculations so much that a Russian roulette variance reduction technique is automatically turned on to prevent prohibitive running times. This attenuation takes a lot of computer time because the transmission from the event to the detector requires transport through each intervening cell and material. The non-attenuation terms of the next-event estimator are known at the event before transmission to the detector

$$ Wp\left(\mu \right)/2\pi {R}^2. $$

If these non-attenuation terms are less than a specified value, k, then Russian roulette is played on the transmission to the detector. That is, if

$$ X= Wp\left(\mu \right)/2\pi {R}^2<k, $$

then the transmission is rouletted with probability

1 – X / k

The transmission is continued with weight W′ = W k/X with probability X/k.

The value of k is specified on the DD input card along with a parameter, M, controlling diagnostic prints (MCNP manual section 3.3.6.11). The default for k is k = 0.1Φ; that is, the default Russian roulette cutoff value is a tenth of the average value of all previous Φ estimates. This default usually results in one transmission to the detector for each source history. If k = 0, then no Russian roulette is played, and all source and collision events cause a time-consuming transmission to the detector. If k < 0, then Russian roulette is based on the constant k. If k > 0, then Russian roulette is based on the current value of kΦ.

Whereas the value of Φ changes with each history, the Russian roulette—and hence, transport—are not identical sampling schemes, and the Central Limit Theorem of statistics is violated. Though a mathematician would argue that the estimated relative errors are then wrong, the effect is negligible. Nevertheless, it is highly recommended that a short calculation be made to determine the average value of Φ and then a tenth of this value be used as the negative value of k. Setting k negative no longer violates the Central Limit Theorem.

Point detectors are specified as

F5:p X Y Z R0

  • X, Y, Z = point detector location coordinates

  • R0 = radius of exclusion for collisions close to detector

The MCNP code also allows ring detectors, which are more efficient but useful only for geometries rotationally symmetric about a major axis. The code also allows radiography tallies, which use the point detector repeated over a grid to model radiographical images. Also the MCNP code uses the same next-event estimator point detector mechanics for DXTRAN spheres. These are illustrated in Example 2.22.

Example 2.21 illustrates many of the above issues.

Example 2.21 Point Detector Example Using Scintillator Geometry of Example 2.20

NE213 scintillator 11 1 -0.9 -11 imp:n=1 20 0 -100 11 imp:n=1 22 0 100 imp:n=0 11 RCC 0 0 0 0 0 3.81 2.54 100 SO 100 sdef pos 0 0 -1 erg .5 vec 0 0 1 dir 1 ARA=1 m1 1001 1.213 6000 1 nps 1e5 print -85 -86 -162 -70 -98 prdmp 2j -1 MODE p PHYS:p 2J 1 FC15 Photon point detector F15:P 0 0 50 0 DD15 -3e-7 500

Example 2.21 uses the scintillator geometry and materials of Example 2.20 and converts to a photon-only problem (MODE p). The source has been modified to be 500 keV, ERG = 0.5 MeV—in the range of many coherent scatters. Because the source is mono-directional (vec 0 0 1 dir 1), the ARA area parameter must be added, or the MCNP code terminates with the following fatal error:

fatal error. sdef requires ara if dir=constant with detectors or dxtran.

A point detector 51 cm in the z-direction from the source point at 0,0,-1 has been added.

F15:P 0 0 50 0

If the PHYS:P and DD15 cards are omitted, then the coherent scattering causes unacceptable convergence with large relative errors. The first 600 large detector scores print a diagnostic. A large detector score is > MΦ or > −Mk, which is input on the DD card. These messages appear between the 1st 50 particle histories and the problem summary table in the output file. When the PHYS:P and DD cards are omitted in the problem of Example 2.21, the first two lines are

det t wgt psc amfp ddetx 1 4.1480E-01 1.0000E+00 5.8540E+03 3.4184E-02 4.6590E+01 radius erg cell nps nch p nrn ipsc 0.0000E+00 5.0000E-01 11 1126 1 7 6

where

  • det = 1, the 1st detector (F15) in the problem;

  • t = 4.1480E-01 = detector tally estimate;

  • wgt = 1.0 = W = particle weight at event;

  • psc = 5.8540E + 03 = p(μ) density function value >> 1.0;

  • amfp = σt*x = sum of macroscopic cross section times transmission distance through each cell of geometry between event and detector;

  • ddetx = 4.6590E+01= R = distance from event to detector;

  • radius = R0 = exclusion radius for collisions close to detector;

  • erg = 5.0000E-01 = particle energy (500 keV);

  • cell = 11 = cell where event happened;

  • nps = 1126 = history number of this large score;

  • nch = 1 = number of collisions so far in this history;

  • p = particle type (photon);

  • nrn = 7 = random number in this history of the event; and

  • ipsc = 6 = scattering law of this collision (coherent).

When the PHYS:p 2J 1 is added to the input file, coherent photon scatter is turned off. The large score diagnostics for scores 1000 times the average do not appear because the largest scores are only 100 times the average. The average tally per history is 2.85337E-06. Convergence is acceptable.

The detector diagnostics card controls both the detector Russian roulette and the large scores diagnostics:

DD15 -3e-7 500

The DD card plays Russian roulette on all scores < 3.e−7, which was chosen to be 0.1 times the average score of 1.85337E−06 of the run without the DD card. The large score diagnostics are printed for any scores >500 × 3e−7 = 1.5e−4.

The tally 15 printout and diagnostics are

1tally 15 nps = 100000 Photon point detector tally type 5 particle flux at a point detector. units 1/cm**2 particle(s): photons detector located at x,y,z = 0.00000E+00 0.00000E+00 5.00000E+01 2.84992E-06 0.0158 detector located at x,y,z = 0.00000E+00 0.00000E+00 5.00000E+01 uncollided photon flux without coherent scattering 1.06682E-07 0.1217 detector score diagnostics cumulative tally cumulative fraction of per fraction of maximum score transmissions transmissions history total tally 3.00000E-07 19229 0.71919 6.00000E-12 0.00000 6.00000E-07 4 0.71934 1.60652E-11 0.00001 1.50000E-06 6 0.71956 6.06829E-11 0.00003 3.00000E-06 1 0.71960 1.65946E-11 0.00003 6.00000E-06 6 0.71983 2.57346E-10 0.00013 1.50000E-05 453 0.73677 6.36520E-08 0.02246 3.00000E-05 3851 0.88080 8.08683E-07 0.30622 3.00000E-04 3187 1.00000 1.97722E-06 1.00000 1.00000E+36 0 1.00000 0.00000E+00 1.00000 average tally per history = 2.84992E-06 largest score = 2.86963E-04 (largest score)/(average tally) = 1.00692E+02 nps of largest score = 53475 score contributions by cell cell misses hits tally per history weight per hit 1 11 8707 26737 2.84992E-06 1.06591E-05 2 20 100000 0 0.00000E+00 0.00000E+00 total 108707 26737 2.84992E-06 1.06591E-05 score misses russian roulette on pd 0 psc=0. 103806 russian roulette in transmission 4798 underflow in transmission 103 hit a zero-importance cell 0 energy cutoff 0

Note that only 0.00013 of the total score comes from the 71.983% of the transmissions with a score less than 6e−6. Thus, a better DD choice would be

DD15 -6e-6 20

Whereas the largest score is 2.86963E−04, it is useful to have the large score diagnostics for scores larger than 20 × 6e−6 = 1.2e−4. The largest/average score is 100, which is good. If this ratio is much larger, there are usually convergence problems.

Whereas the Russian roulette criteria are identical for all histories, the history of the largest score, nps = 53475, can be rerun to see what made the score large. If the default DD is used, then the criteria for the Russian roulette change for each history and are unknown and thus histories cannot be rerun.

Note that “score contributions by cell” can be used to play Russian roulette on the contributions from different cells using a PD card; cells with a low weight per hit often should have PD values <1.

One might ask how having each source or collision event score directly to the point detector gives a believable flux estimate. The next-event estimator is robust and accurate and any claim to the contrary is probably user error. Any score that is made in any tally comes from the last event preceding it. The next-event estimator is justified by calculating the next event from a source or collision that is the last event before making a score.

2.5 Plotting

The MCNP code supports plotting with the X Windows graphics system. The window manager on the display device must support the X protocol with software such as Xming. (Linux supports the X protocol natively).

The MCNP code makes geometry, cross-section, and tally plots. The geometry plots can be interactive or command driven. The cross-section and tally plots are command driven. Mesh tally plots overlay a mesh of tallies over the geometry plot. All MCNP plots can be made from the RUNTPE continuation file or during the execution of a problem. Tally plots can also be made from MCTAL files, which are created by the PRDMP 3rd entry. (See MCNP Manual [2] section 3.3.7.2.3).

2.5.1 Geometry Plotting and Command Files

All MCNP plots can use command files. Whenever an the MCNP code is run in the plotting mode, a COMOUT file (with an appended c if the NAME option is used) is created. This file can then be used or edited to make a series of plots. For example, consider the MCNP input file “ex11_6a” of Example 2.6. The geometry plots of Fig. 2.9 through Fig. 2.13 through can be reproduced with the following MCNP execution command:

MCNP6 i=ex11_6a n=plot11_6a. com=com11_6a ip

The command file, “com11_6a” is

py 6 or 10 6 0 ex 20 la 0 0 pause lev 0 la 2 2 cel mbody off pause lev 1 la 0 2 cel pause lev 2 la 0 0 pause OR 10 6 -5 PY 6 EX 7 LEV 2 LA 0 1 IJK pause pause end end

All plot commands can be shortened to two or three characters (until non-ambiguous). The commands PY, ORigin, Extent, Label, and LEVel are explained in the MCNP manual and also listed in the plot help package by typing “?”, “help,” or “option” when in the interactive plotter. The PAUSE command stops the plot and sends the following message to the plot terminal:

strike any key to continue. quit returns to interactive plot.

Striking any key (we prefer “enter”) makes the next picture up to the next pause command. Typing “quit” returns to the interactive plot mode.

For geometry plots, the help command lists only possible commands; the same plots may be made interactively. For cross-section and tally plots, the MCPLOT help command provides a complete description of each command and its options.

For geometry plots, the command file commands may be entered in the command window, in the lower right of the plot window where it says “click here or picture or menu,” or by clicking the interactive buttons. The equivalent interactive buttons for the com11_6a file are: ZX, OR, ZOOM 5, L1; then LEV, MBODY, L1, L2; LEV, L1; then LEV, L2; then IJK (on right margin), L2, ZOOM 3, OR.

Note that the NAME=plot11_6a command will cause a new COMOUT command file, plot11_6a.c, to be made. The COMOUT file (or .c file) can then be edited to reproduce or make even fancier plots with labels, titles, scales, etc.

While the MCNP code is running transport, the geometry plotter can be called by <ctrl-c> M to get the MCPLOT tally plot and then “PLOT” to get the geometry plotter. The PAC and N buttons (lower right), followed by L2 for cell quantities, can be used to overlay the geometry with the number of particles that enter, population, number of collisions, etc., of each cell in the geometry plot. N refers to the column number in the MCNP output file particle activity by cell (PAC) printout.

The command “LA 2 2” doubles the size of cell and surface labels. Other geometry plot buttons are

  • UP, RT, DN, LF (on top menu) to go up, right, down, or left one frame;

  • REDRAW on the bottom menu is needed to refresh the picture when the refresh is not done automatically;

  • CURSOR to enclose a region on which to zoom in;

  • CellLine to not outline cells by drawing surfaces;

  • COLOR to toggle surface-only plots without color and to color by other cell quantities first clicked on the right side buttons. For example, “VOL” then “COLOR” will make the cells colored proportionately to cell volume;

  • SCALES to get a grid of dimensions in centimeters;

  • LEVEL to toggle repeated structures levels;

  • XY, YZ, ZX to get PZ, PX, and PY plots;

  • L1, L2 to display surface and cell labels;

  • MBODY to toggle macrobody facets on and off; and

  • LEGEND to get color scaling.

If the plot is moved or otherwise goes blank, clicking in the upper left corner will restore the plot. When entering commands from the terminal, the command “INTERACT” will return to the interactive plot mode. In the interactive plot mode, clicking the “PLOT>” button on the bottom returns to the command mode. Terminating the plotting session should be done by clicking “END” on the bottom menu because closing the plot window will cause the MCNP code to crash and lose files.

One way to print a plot to another document without using the interactive buttons, first go to the command mode on the terminal. If in the interactive mode, click “PLOT>” on the bottom menu. From the command mode, enter the command “VIEWPORT SQUARE” and only the square plot will appear. Click inside the picture and then take a screenshot (<ctrl-alt><print screen> on Windows systems) and then click inside the document <ctrl-v> (on Windows). To exit the command prompt mode and return to the interactive plotter, enter “INTERACT.” Alternatively the “file” command can be used to produce a postscript format file.

2.5.2 Cross-Section Plotting

MCNP cross sections can be plotted for neutrons, photons, photonuclear data, protons, and other cross-section libraries in ACE (A Compact ENDF) format [10]. MCNP stopping powers for electrons can also be plotted. The cross-section plotter is command driven, not interactive. The MCNP code is run in the IXZ mode.

To reproduce all the figures in Sect. 2.2, “Materials and Cross Sections,” the MCNP code is run as

MCNP6 IXZ i=EX152 N=J152. COM=COM152

In this example, EX152 is the input file of Example 2.11 with MODE P E changed to MODE N P E. The MCNP code will produce files with names J152.o, J152.c, etc. The plot command file is COM152:

xs m313 cop xs 92238.80c cop xs 92235.80c pause xlim 1e-11 1e-5 xs 6000.80c mt 1 cop xs grph.10t mt 1 cop mt 2 cop mt 4 pause xlim .001 100 ylim .01 1e7 par p mt -5 xs 1000.04p cop xs 8000.04p & cop xs 40000.04p cop xs 92000.04p pause ylim .01 1e7 par p xs 92000.04p mt -5 cop mt -1 cop mt -2 cop mt -3 cop mt -4 pause ylim 1 200 par e xs m609 cop xs m201 cop xs m313 pause par e xs m313 ylim .5 30 mt 3 cop mt 1 cop mt 2 pause end end

The cross section plot commands are largely the same as for the tally plotter (see MCNP manual [2] section 5.3). The commands relating to the picture such as XLIM, YLIM, etc., are the same; only commands unique to tallies or data are different. Once a command is set, it holds for all subsequent plots until specifically changed. RESET ALL resets to the default commands. XS, PAR, and MT are specific for the cross-section plotter. PAR defaults to neutrons, and if neutrons are not in the MCNP input file, nothing will be plotted until PAR is set to a particle type on the MODE card of the input file. All the possible plot and tally commands are listed when “?”, “help,” or “options” is entered. Each command is described in detail by entering the command name after the “?”, “help,” or “options.” For example, “help XS” prints on the terminal:

xs > Syntax: xs m Plot a cross-section according to the value of m. Option 1: m=Mn, a material card in the INP file. Example: XS M15. The available materials will be listed if a material is requested that does not exist in the INP file. Option 2: m=z, a nuclide ZAID. Example: XS 92235.50c. The full ZAID must be provided. The available nuclides will will be listed if a nuclide is requested that does not exist in the INP file. See "mt", "par".

“Help MT” describes how to get specific reactions in the cross-section plotter:

mt > Syntax: mt m Plot reaction n of material XS m. The default is the total cross section. The available reaction numbers will be listed if one enters a reaction number that does not exist (e.g., 999). See "par", "xs".

Note that only cross sections and reactions specified in the MCNP input file may be plotted and they must be fully specified. XS=1000 fails; XS=1000.04p is required in this example. If a cross section is requested that is not available, then the plotter will list the ones available:

xs 100 can't find zaid 100 available zaids: 1001.80c 6000.80c 7014.80c 8016.80c 40090.80c 40091.80c 40092.80c 40094.80c 40096.80c 92235.80c 92238.80c benz.10t grph.10t lwtr.10t poly.10t 1000.04p 6000.04p 7000.04p 8000.04p 40000.04p 92000.04p

If electrons have been loaded, they will have a similar ZAID set as the photons, with p changed to e. (The electron libraries are 03e and 01e.)

Similarly, when plotting a nuclide like xs=92235.80c (but not a material like xs=m311), specifying an unavailable MT reaction number causes the MCNP cross-section plotter to print the available reactions:

mt -12 available (not expunged) endf reactions: for isotope 92235.80c zaid 92235.80c has : mt = 16. mt = 17. mt = 18. mt = 37. mt = 51. mt = 52. mt = 53. . . . mt = 89. mt = 90. mt = 91. mt =102. mt = 4. some mts may have been expunged or may not exist in the material or zaid cross sections given. if an mt exists but is not listed here, call it out with an fm tally card.

The principal reaction MT values are shown in Table 2.1.

The S(α,β) reactions are 1 = total, 2 = elastic, 4 = inelastic.

The neutron cross sections do not include thermal effects. The photon cross sections do not include form factors for coherent and incoherent scattering.

In command file COM152, the other commands are:

  • COPlot = overlay the next plot on the present one to plot multiple curves;

  • PAUSE = do not advance to the next plot until a key is struck. If a number, N, follows, then hold the plot for N seconds;

  • XLIMs = set the range of x-axis values; and

  • YLIMs = set the range of y-axis values.

Model cross sections cannot be plotted. Secondary energy and angle distributions cannot be plotted.

2.5.3 Tally Plotting

All MCNP tallies, most KCODE criticality quantities, and statistical information can be plotted with the MCNP tally plotter. The plots can be made during the course of a calculation either by interrupt, <ctrl-c><m>, or periodically with the MPLOT run-time plotter (MCNP manual [2] section 3.3.7.2.5), which sends a tally plot to the terminal after specified numbers of histories or KCODE cycles. After a calculation, the RUNTPE file can be used to produce geometry, cross section, or tally plots and the MCTAL file can be used for all tally plots. These capabilities also apply to MESH tallies described in the next section.

RUNTPE continue files can be read only by the version and operating system that wrote them. MCTAL files are ASCII files in a backward-compatible format, so MCTAL files written by different versions of the MCNP code on different operating systems many years ago may still be used to plot tallies with the latest MCNP version. Furthermore, MCTAL files are much shorter and easier to store. Different MCNP runs may be plotted against each other from any combination of MCTAL and RUNTPE files.

Tally plots are made with the command

MCNP6 Z

after which the MCNP code will ask for a RUNTPE file name (RUNtpe =) or an MCTAL file name (RMCtal =).

The following plot command file, COM153, plots all tallies of Example 2.20 and Example 2.21 with MCTAL files j142.m and j143.m:

rmc=j142.m tal 6 xlim .2 20 loglog la "light output" & title 1 "Proton Heating Tally with DE/DF dose function" pause tal 26 xlim 0 10 ylim .03 .2 linlog la "energy deposition" & Title 1 "Proton Energy Deposition" pause tal 8 xlim .01 10 ylim .001 .2 linlog la "light output pulses" & title 1 "Based upon light output F6: FT8 PHL 1 6 1" pause tal 18 ylim .005 .04 la "energy deposition pulses" & title 1 "Based upon energy deposition F16: FT8 PHL 1 16 1" pause tal 28 xlim 0 10 ylim .005 .03 title 1 "Simple pulse-height F8 tally" pause tal 18 Title 1 "Comparison of F8 and FT8 PHL 16" la "FT8 PHL 1 16 1" & cop tal 28 la "Simple F8" pause tal 38 ylim .001 .1 title 1 "*F8 Energy-Deposition" la "energy deposition" pause reset all rmc j143.m tfc m title 1 "TFC M --- Mean Value Convergence" la "mean" pause end end

Figure 2.25 shows the plot of the F6 tally in Example 2.20:

Fig. 2.25
figure 25

Tally plot of F6 tally in Example 2.20

The F6 tally of Example 2.20 is repeated here to accompany the plot of Fig. 2.26.

FC6 Track length heating converted to light output f6:h 11 $ F6 is converted to light DE6 LIN 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 9 9.2 9.4 9.6 9.8 10 DF6 LIN 0.00 0.01 0.04 0.08 0.12 0.17 0.23 0.29 0.35 0.42 0.50 0.57 0.66 0.74 0.83 0.93 1.02 1.12 1.23 1.34 1.45 1.56 1.67 1.79 1.92 2.04 2.17 2.30 2.43 2.57 2.70 2.84 2.99 3.13 3.28 3.43 3.58 3.74 3.90 4.05 4.22 4.38 4.55 4.72 4.89 5.06 5.23 5.41 5.59 5.77 5.95 e0 0.2 48I 10 11 13 15 17 20

Fig. 2.26
figure 26

Tally plot of F26 tally in Example 2.20. Note how the eight tally dimensions—F, D, U, S, M, C, E, T—are listed on the right

The tally plot commands are

rmc=j142.m tal 6 xlim .2 20 loglog la "light output" & title 1 "Proton Heating Tally with DE/DF dose function"

Figure 2.26 shows the plot of the F26 tally in Example 2.20.

FC26 Track length MeV deposited f26:h 11 sd26 1 $ Divide by 1 rather than grams e0 0.2 48I 10 11 13 15 17 20

The tally plot commands are

rmc=j142.m tal 26 xlim 0 10 ylim .03 .2 linlog & la "energy deposition" Title 1 "Proton Energy Deposition"

Figure 2.27 shows the plot of the F8 tally in Example 2.20.

Fig. 2.27
figure 27

Tally plot of F8 tally in Example 2.20. Note how the eight tally dimensions—F, D, U, S, M, C, E, T—are listed on the right

FC8 Total light output f8:h 11 $ F8 is total light output FT8 PHL 1 6 1 0 e8 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50

Note that the energy deposition for this F8 tally is from the heating tally F6 and does not require energy bins for the F6 tally.

The tally plot commands are

rmc=j142.m tal 8 xlim .01 10 ylim .001 .2 linlog & la "light output pulses" & title 1 "Based upon light output F6: FT8 PHL 1 6 1"

Figure 2.28 shows the plot of the F18 tally in Example 2.20.

Fig. 2.28
figure 28

Tally plot of F18 tally in Example 2.20. Note how the eight tally dimensions—F, D, U, S, M, C, E, T—are listed on the right

FC18 Total energy deposited F18:h 11 Ft18 PHL 1 16 1 0 $ 0.37 1.32 and 3.79 are expt thresholds e18 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50

Note that the energy deposition for this F8 tally is from the heating tally F16 and does not require energy bins for the F16 tally.

The tally plot commands are

rmc=j142.m tal 18 ylim .005 .04 la "energy deposition pulses" & title 1 "Based upon energy deposition F16: FT8 PHL 1 16 1"

Figure 2.29 shows the plot of the simple (no FT special treatment) F28 tally in Example 2.20.

Fig. 2.29
figure 29

Tally plot of F28 tally in Example 2.20. Note how the eight tally dimensions—F, D, U, S, M, C, E, T—are listed on the right

FC28 Pulse height tally f28:h 11 e28 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50

Note that the energy deposition for this F8 tally is from a surface/collision estimator, which adds the energy entering the cell and produced by secondary particles in the cell and then subtracts out the energy of particles leaving the cell.

The tally plot commands are

rmc=j142.m tal 28 xlim 0 10 ylim .005 .03 & title 1 "Simple pulse-height F8 tally"

The name of the file is “j142.m,” the name of the MCTAL file being read. The LAbel command should be used to get a better label.

Figure 2.30 shows a coplot of tallies F18 and F28 and demonstrates the capability to plot multiple curves with the MCNP plotter.

Fig. 2.30
figure 30

Plot of Example 2.20 F18 and F28 tallies plotted together. Note how the eight tally dimensions—F, D, U, S, M, C, E, T—are listed on the right

FC18 Total energy deposited F18:h 11 Ft18 PHL 1 16 1 0 $ 0.37 1.32 and 3.79 are expt thresholds e18 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50 FC28 Pulse height tally f28:h 11 e28 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50

Note that the pulses of the F18 pulse-height light tally energy deposition are from the F6 track length heating tally. The pulses of the F28 simple pulse-height tally are from the energy balance of surface and collision estimators. The two tallies agree within statistics despite being generated from totally independent estimators.

The tally plot commands are

tal 18 Title 1 "Comparison of F8 and FT8 PHL 16" & la "FT8 PHL 1 16 1" cop tal 28 la "Simple F8"

Figure 2.31 shows tally *F38, an energy deposition pulse-height tally:

Fig. 2.31
figure 31

Plot of Example 2.20 *F38 energy deposition pulse-height tally. Note how the eight tally dimensions—F, D, U, S, M, C, E, T—are listed on the right

FC38 Pulse height energy deposition *f38:h 11 e38 0 1e-6 .001 .05 5I .35 .37 .4 9I 1.0 1.2 1.32 1.4 10I 3.6 3.79 3.8 30I 10 6I 50

The total tally bin is 0.908 ± .01 MeV, which is the same (within statistical uncertainty) energy deposition, 0.905 ± .01 MeV, as the F26 track length tally with the SD divisor to get units of MeV:

FC26 Track length MeV deposited f26:h 11 sd26 1 $ Divide by 1 rather than grams

The tally plot commands are

rmc=j142.m tal 38 xlim 0 10 ylim .001 .1 & title 1 "*F8 Energy-Deposition" la "energy deposition"

Figure 2.32 shows how it is possible to plot tallies with no energy, time, cosine, etc., bins—that is, tallies with a single number as the total tally.

Fig. 2.32
figure 32

Plot of convergence of Example 2.21 point detector tally

It also shows tally convergence of the point detector tally in Example 2.21:

FC15 Photon point detector F15:P 0 0 50 0 DD15 -3e-7 500

The convergence and all convergence criteria of the 10 statistical tests can be plotted with the TFC plot command:

rmc j143.m tfc m title 1 "TFC M --- Mean Value Convergence" la "mean"

Note that because “TITLE 2” is not set, the second title line is carried forward from a previous plot.

Figure 2.33 shows the KCODE criticality plot of the principal estimators vs the cycle number using the MCNP tally plotter for the problem of Example 2.14.

Fig. 2.33
figure 33

KCODE Criticality plot

The plot COM file is

run=j133.r ylim .794 .8 kc 16 la "CAT" & title 1 "Criticality Convergence" title 2 "Estimators vs Cycles" & cop kc 11 la "Collision" cop kc 12 la "Absorption" & cop kc 13 la "Track Length"

The various KCODE commands are listed when running the MCNP code in the tally plot mode and entering “HELP KCODE.” Here, the multiple curves KC = 16, 11, 12, and 13 are plotted together. KC = 16 is the cumulated covariance-weighted collision-absorption-track length (CAT) estimator of keff; KC = 11 is the cumulative collision estimator of keff; KC = 12 is the cumulative absorption estimator of keff; and KC = 13 is the track length estimate of keff. Because Example 2.14 did not have a PRDMP 3rd entry to create a MCPLOT file, the RUNTPE file, RUN=j133.r, has to be used.

Run-Time Plotting: MCNP plots may also be generated during the course of a run with the MPLOT input card. For Example 2.14, add the following to the MCNP input file:

MPLOT freq 5 kc 16 cop kc 11 cop kc 12 cop kc 13

FREQ=5 puts up a plot of the calculation every 5 KCODE cycles. For non-KCODE calculations, FREQ is the number of histories between which plots will be sent to the computer terminal. MPLOT arguments are MCPLOT commands. A user can visualize source convergence with an FMESH plot (see next section) called on MPLOT, for example.

2.5.4 Mesh, Radiography, and Ring Tallies

Example 2.22 illustrates a number of tally capabilities. It also uses the vertical input format for the source definition with the # operator (MCNP manual [2] section 2.8.2).

Example 2.22 NE-213 Detector with Point, Ring, Radiography, and Mesh Tallies

NE213 scintillator 11 1 -0.9 -11 imp:n=1 20 0 -100 11 imp:n=1 22 0 100 imp:n=0 11 RCC 0 0 0 0 0 3.81 2.54 100 SO 100 MODE n h m1 1001 1.213 6000 1 mx1:h J 6012 nps 1e5 print -85 -86 -162 -70 -98 prdmp 2j 1 c SDEF PAR N x=d21 y=d22 z -50 erg 100 vec 0 0 1 dir 1 ARA=56 # SI21 SP21 SI22 SP22 -4 0 -3.5 0 4 1 3.5 1 c TALNP 55 65 FC15:n Point detector at side F15:n 4 0 2 0 dd15 -3e-4 FC25:n Ring detector at side FZ25:n 2 4 0 dd25 -3e-4 FC55 Pinhole image projection radiograph at side FIP55:n 0 50 2 0 0 0 2 5 0 40 c55 -3 99i 3 FS55 -3 99i 3 FC65 Transmitted image projection radiograph on axis FIR65:n 0 0 40 0 0 0 0 0 0 0 c65 -20 99i 20 FS65 -20 99i 20 FC75 Transmitted image projection radiograph on side FIR75:n 4 0 0 0 0 0 0 0 0 0 c75 -10 99i 30 FS75 -20 99i 20 c FC104 MCNP5 style neutron FMESH tally FMESH104:n geom=rec origin -5 -5 -30 IMESH 5 IINTS 100 JMESH 5 JINTS 100 KMESH 30 31 KINTS 1 1 FC114 MCNP5 style proton FMESH tally FMESH114:h geom=rec origin -5 -5 -30 IMESH 5 IINTS 100 JMESH 5 JINTS 100 KMESH 30 31 KINTS 1 1 FC124 MCNP5 style neutron source FMESH tally FMESH124:n geom=rec origin -5 -5 -60 type source IMESH 5 IINTS 100 JMESH 5 JINTS 100 KMESH 60 KINTS 1 c TMESH RMESH201:n FLUX TRAKS CORA201 -5 99I 5 CORB201 -5 99I 5 CORC201 -30 30 31 RMESH211:h FLUX TRAKS CORA211 -5 99I 5 CORB211 -5 99I 5 CORC211 -30 30 31 RMESH232 N CORA232 -5 99I 5 CORB232 -5 5 CORC232 -60 30i -.1 60i 4 10i 20 ENDMD MPLOT FREQ 1e4 plot pz 1 ex 5.5 la 0 1 tal201.1 col on

The following plot COMmand file shows how to plot the F15 point, FZ25 ring, F55, F65, F75 radiography tallies and FMESH104, FMESH114, and FMESH124 tallies. Plotting the TMESH tallies is problematic and is dealt with thereafter.

rmc j161.m tal 15 la "F15 point" tfc m cop tal 25 la "FZ25 Ring" pause tal 55 free sc pause tal 65 free sc pause tal 75 free sc pause run j161.r plot fmesh 104 pause fmesh 114 pause fmesh 124 pause pause end end

Figure 2.34 is a coplot of the F:15 point detector tally and FZ25 ring detector tally from Example 2.22.

Fig. 2.34
figure 34

Coplot of F15 point detector and equivalent FZ25 ring detector tallies. The plot command is “rmc j161k.m tal 15 la “F15 point” tfc m cop tal 25 la “FZ25 Ring””

FC15:n Point detector at side F15:n 4 0 2 0 dd15 -3e-4 FC25:n Ring detector at side FZ25:n 2 4 0 dd25 -3e-4

The point and ring detectors converge within statistical uncertainty because of the cylindrical symmetry of the problem. In such problems, the ring detector is more efficient and converges with a higher FOM. Note that the SDEF source requires the ARA entry to give the area of the plane-wave source and enable next-event estimators such as the point detector, ring detector, and radiography tallies.

The pinhole radiograph tally (FIP) is

FC55 Pinhole image projection radiograph at side FIP55:n 0 50 2 0 0 0 2 5 0 40 c55 -3 99i 3 FS55 -3 99i 3

The FIP55 tally is a pinhole radiography tally. The pinhole is like the focal point in a camera. The pinhole radiography tally captures the image of a particles emitted from or reflected by a region. The first three entries are the coordinates of the pinhole (0,50,2). The 4th entry is unused. The 5th–7th entries are the center of the object (0,0,2). The tally axis is the vector from the object (0,0,2) to the pinhole, (0,50,2)—namely the y-axis—normal to the incident beam axis. No collision or source points outside a 5 cm radial collimator (8th entry) contribute to the radiograph, saving computational time. The 9th entry is the radius of the pinhole (0=perfect pinhole). The image plane is 40 cm (10th entry) beyond the pinhole in the tally axis (y-axis) direction. The s- and t-axes are x and y lengths of the radiograph plane and are specified on the FS (segment) and C (cosine) cards as dimensions in cm. TALNP prevents excessive printout to the OUTP file. The plot is illustrated in Fig. 2.35.

Fig. 2.35
figure 35

Pinhole radiography tally F55 of Example 2.22. The image plane is a side (y-axis) view of NE213 detector secondary and scattered neutrons. The plot command is “tal 55 free sc”

Any tally can be plotted as a two-dimensional array of any two of the 8 tally dimensions. For radiography tallies, it is convenient to plot the cosine and segment dimension. The “free” tally plot command can be in one dimension as in “free s” or “free c,” or in two dimensions as in “free sc.”

Example 2.22 has two transmitted image (FIR) rectangular coordinate radiography tallies. Cylindrical transmitted image radiography tallies, FIC, are available but not illustrated.

FC65 Transmitted image projection radiograph on axis FIR65:n 0 0 40 0 0 0 0 0 0 0 c65 -20 99i 20 FS65 -20 99i 20 FC75 Transmitted image projection radiograph on side FIR75:n 4 0 0 0 0 0 0 0 0 0 c75 -10 99i 30 FS75 -20 99i 20

The first three entries are the coordinates of the center of the image plane, which is normal to the tally axis. The 4th entry is unused. The tally axis is the vector from the image center (5th–7th entries) (0,0,0) to the image plane center [(0,0,40) for FIR65 and (4,0,0) for FIR75]. The last 8th entry (0) tallies both direct and indirect (scattered) contributions; = 1 for scattered only. The 9th entry is the radial field of view (= 0 for infinite). The 10th entry scores to the center of each s- and t-axis grid point rather than to a random point in the grid.

The transmitted image radiography tallies of Example 2.22 are plotted below: Fig. 2.36 is the image from the F65 on-axis tally; Fig. 2.37 is the image from the F75 side-view tally.

Fig. 2.36
figure 36

Transmitted image radiography tally F65 of Example 2.22. The image plane is in the SDEF source direction (on z-axis) and behind the object. The tally plot command is “tal 65 free sc”

Fig. 2.37
figure 37

Transmitted image radiography tally F75 of Example 2.22. The image plane is a side view normal to the x-axis, normal to the source direction z axis direction, and shows a side view of neutron emission from the object. The tally plot command is “tal 75 free sc”

The MCNP5 style mesh tallies in Example 2.22 show the neutron flux (FMESH104), proton flux (FMESH114), and neutron source (FMESH124) integrated over the z-axis source direction and in the x–y plane:

FC104 MCNP5 style neutron FMESH tally FMESH104:n geom=rec origin -5 -5 -30 IMESH 5 IINTS 100 JMESH 5 JINTS 100 KMESH 30 31 KINTS 1 1 FC114 MCNP5 style proton FMESH tally FMESH114:h geom=rec origin -5 -5 -30 IMESH 5 IINTS 100 JMESH 5 JINTS 100 KMESH 30 31 KINTS 1 1 FC124 MCNP5 style neutron source FMESH tally FMESH124:n geom=rec origin -5 -5 -60 type source IMESH 5 IINTS 100 JMESH 5 JINTS 100 KMESH 60 KINTS 1

These tallies can be plotted from the MCNP RUNTPE continuation file in the geometry plot mode simply by clicking the FMESH button in the lower left of the MCNP geometry plot or by entering commands at the command prompt or in a COMmand file. The MCNP code automatically determines the view as through the middle of the largest mesh span. Figure 2.38 shows the FMESH104 neutron mesh tally. Figure 2.39 shows FMESH114 proton mesh tally. Figure 2.40 shows the FMESH124 neutron source mesh tally. In each of these plots, the PLOT> command “VIEWPORT SQUARE” from the terminal was used to get the plot without the interactive geometry plot buttons showing. The command “VIEWPORT RECT” restores the rectangular plot plane and “INTERACT” switches from PLOT> command mode back to interactive mode.

Fig. 2.38
figure 38

FMESH neutron flux tally F104, which excludes the source beam and consists of secondary neutrons. The plot is made by clicking the FMESH button on the lower left of the geometry plot or by entering the command “FMESH 104”

Fig. 2.39
figure 39

FMESH proton flux tally F114. The protons are secondary particles from neutron reactions. The plot is made by clicking the FMESH button on the lower left of the geometry plot or by entering the command “FMESH 114”

Fig. 2.40
figure 40

FMESH neutron source tally F124. The neutrons are from proton collisions. The z = 0 plot plane is coincident with the bottom of the NE213 detector object, resulting in the semi-circle surface boundary. The plot is made by clicking the FMESH button on the lower left of the geometry plot or by entering the command “FMESH 124”

The MCNPX style TMESH tallies in Example 2.22 enable more plotting capabilities than the FMESH tallies. Tracks, populations, doses, DXTRAN, point detector, multiple particle energy deposition, and more can be plotted in addition to fluxes and source points. Unfortunately, the interface is inconsistent with the FMESH interface and harder to use. The TMESH tallies start with “TMESH” and end with “ENDMD” and must be all together:

TMESH RMESH201:n FLUX TRAKS CORA201 -5 99I 5 CORB201 -5 99I 5 CORC201 -30 30 31 RMESH211:h FLUX TRAKS CORA211 -5 99I 5 CORB211 -5 99I 5 CORC211 -30 30 31 RMESH232 N CORA232 -5 99I 5 CORB232 -5 5 CORC232 -60 30i -.1 60i 4 10i 20 ENDMD

There are also bugs that make entering PLOT> commands ineffective. Entering commands from a COMmand file, from the terminal, or in the lower left of the interactive geometry plot where it says “click here or picture or menu” does not work unless the lower right button “TAL” in the interactive geometry plot reading from the RUNTPE is first toggled. After that, the “TAL,” “PAR,” and “N” buttons can be used or the TMESH plot commands may be entered. Note that the FMESH button must be clicked off. (There are TMESH plotting examples in the MCNP manual [2] Sect. 6.4.4.)

The neutron flux and particle track rectangular RMESH201 mesh tallies are shown in the TMESH plots of Figs. 2.41 and 2.42. The XY view is normal to the source z-axis. The proton flux and particle track rectangular RMESH211 mesh tallies are shown in the TMESH plots of Figs. 2.43 and 2.44. Again, the XY view is normal to the source z-axis. The neutron source rectangular RMESH232 mesh tally is shown in the TMESH plot of Fig. 2.45. The XZ view is normal to y-axis with the source z-axis in the vertical direction.

Fig. 2.41
figure 41

TMESH plot of rectangular mesh tally RMESH201 neutron flux made from RUNTPE in the geometry plot mode. The view is down the source z-axis, “pz 1 ex 5.5 or 0 0 1.” The plot command after clicking off FMESH and cell labels and entering VIEWPORT SQUARE is “la 0 1 tal201.1 la 0 0.” The 0.1 gets the first PAR entry, namely flux. The first LA command gets the tally for coloring and the second turns off cell labels

Fig. 2.42
figure 42

TMESH plot of rectangular mesh tally RMESH201 neutron tracks made from RUNTPE in the geometry plot mode. The plot view is down the source z-axis with geometry plot command “pz 1 ex 5.5 or 0 0 1.” The plot command after clicking off FMESH and cell labels and entering “VIEWPORT SQUARE” is “la 0 1 tal201.2 la 0 0.” The 0.2 gets the second PAR entry, namely TRAKS. The first LA command gets the tally for coloring and the second turns off cell labels

Fig. 2.43
figure 43

TMESH plot of rectangular mesh tally RMESH211 proton flux made from RUNTPE in the geometry plot mode. The plot view is “pz 1 ex 5.5 or 0 0 1,” namely down the source z-axis. In the interactive mode, FMESH must be off, cell labels L2 off, and the lower right “TAL” and “N” buttons must be cycled through until tally 211 is the edit quantity (middle left). Then click PLOT> and enter “VIEWPORT SQUARE” and “la 0 1 tal211.1 la 0 0.” The 0.1 gets the first PAR entry, namely flux. The first LA command gets the tally for coloring and the second turns off cell labels

Fig. 2.44
figure 44

TMESH plot of rectangular mesh tally RMESH211 proton tracks made from RUNTPE in the geometry plot mode.The plot view is “pz 1 ex 5.5 or 0 0 1,” namely down the source z-axis. In the interactive mode, FMESH must be off, cell labels L2 off, and the lower right “TAL” and “N” buttons must be cycled through until tally 211 is the edit quantity (middle left). Then click PLOT> and enter “VIEWPORT SQUARE” and “la 0 1 tal211.2 la 0 0.” The 0.2 gets the first PAR entry, namely flux. The first LA command gets the tally for coloring and the second turns off cell labels

Fig. 2.45
figure 45

TMESH plot of rectangular mesh tally RMESH232 neutron source made from RUNTPE in the geometry plot mode. The plot view is “py 0 ex 34 or 0 0 -25.” Click off FMESH and cell labels. Click on TAL and N until tal232.1 becomes the edit quantity. Click PLOT> and enter “VIEWPORT SQUARE” to get the picture without buttons. The SDEF source is seen in the –Z mesh region at the bottom of the plot. The source plane is a line in this view, but appears like a rectangle with a thickness of the tally mesh. The red color is the maximum source strength. Secondary neutrons from protons are seen at the NE213 detector at the top of the plot. These secondary neutrons appear as blue dots because of their lower intensity

The run-time plotting for TMESH tallies utilizes the MPLOT command in the input file. The plot picture will flash on the computer terminal with a frequency of every FREQ = 1.E–4 histories. The remaining commands specify the plot view and which TMESH tally to plot.

MPLOT FREQ 1e4 plot pz 1 ex 5.5 la 0 1 tal201.1 col on

2.6 Statistics and Convergence

Example 2.23 is a thin, 2 cm-thick disc with a huge 1×1010 cm radius.

Example 2.23 Demonstration of False Convergence and Ten Statistical Tests

Isotropic source in RCC disc 1 0 -1 imp:n=1 2 0 1 imp:n=0 1 RCC 0 0 -1 0 0 2 1e10 nps 2e4 prdmp 2j 1 print c fc11 Neutrons crossing surface f11:n (1.1 1.2 1.3) 1.1 1.2 1.3 ft11 frv 0 0 1 c11 -.95 38i 1 T fq11 C F c fc4 Neutron flux in RCC disc f4:n 1 sd4 1 c SDEF vec 0 0 1 dir=d1 si1 -1 9ilog -1e-10 1e-10 9ilog 1 sp1 0 .9 8ilog 9e-10 2e-10 9e-10 8ilog .9

The problem is contrived to demonstrate the ten MCNP statistical tests that can help identify false convergence. Tally F4 of this problem passes nine of the ten MCNP statistical tests although it is falsely converged after 20000 histories to an estimate of 9.19 ± 5.9%. The correct answer is 24.06 neutrons, more than a factor of two different. Sometimes false convergence gives estimates that are orders of magnitude wrong.

False convergence is caused by undersampling important regions of problem space. Undersampling can occur if an important region of the geometry is not sufficiently sampled. Undersampling also occurs if important physics, such as neutron resonances or photon coherent scatter, are improperly sampled. In Example 2.23, undersampling is caused by insufficiently sampling source angle.

The source in Example 2.23 is an isotropic source because the number of neutrons started in each cosine bin from −1 < μ < 1 is proportional to the size of the cosine bin. This is equivalent to these other forms of an isotropic source:

SDEF

or

SDEF vec 0 0 1 dir=d1 si1 -1 1 sp1 0 1

Figure 2.46 shows that the source is indeed isotropic. It is a plot of the tally F11 neutron current in cosine bins.

Fig. 2.46
figure 46

Tally plot of Example 2.23 tally F11 neutron current in cosine bins; number of particles crossing all surfaces of the disc relative to the z-axis

The tally plot command is

rmc j20k.m tal 11 free c la "20k" linlin cop rmc j8.m tal 11 la "1e8"

The MCTAL file j20k.m was created with NPS 2e4 and the MCTAL file j8.m was created with NPS 1e8. The number of neutrons in each of the 40 equal cosine bins is the same within statistical uncertainty, which is an isotropic distribution.

The output file tally fluctuation chart for the F4 flux tally shows the tally results as a function of the number of histories:

tally 4 nps mean error vov slope fom 2000 8.2767E+00 0.1805 0.5439 1.9 117916 4000 7.8166E+00 0.1082 0.3385 2.0 163975 6000 8.4977E+00 0.0983 0.1613 2.0 132371 8000 8.7945E+00 0.0949 0.1873 2.3 106554 10000 8.8567E+00 0.0820 0.1405 2.2 114229 12000 8.8204E+00 0.0759 0.1160 2.2 133236 14000 8.7715E+00 0.0679 0.1014 2.5 138769 16000 8.7582E+00 0.0617 0.0886 2.4 143993 18000 9.0119E+00 0.0592 0.0715 2.3 136872 20000 9.1865E+00 0.0589 0.0726 2.3 122876

The “mean” is the flux estimate of the answer. The “error” is the relative error, which is the standard deviation of the mean divided by the mean. The “vov” is the variance of the variance, or the relative error of the relative error. The “slope” is the rate of convergence. The “fom” is the Figure of Merit.

The tally printed in the output file is

1tally 4 nps = 20000 Neutron flux in RCC disc tally type 4 track length estimate of particle flux. units 1/cm** cell 1 9.18646E+00 0.0589 results of 10 statistical checks for the estimated answer for the tally fluctuation chart (tfc) bin of tally 4 tfc bin --mean-- ---------relative error--------- behavior behavior value decrease decrease rate desired random <0.10 yes 1/sqrt(nps) observed random 0.06 yes yes passed? yes yes yes yes ============================================================== ----variance of the variance---- --figure of merit-- -pdf- value decrease decrease rate value behavior slope <0.10 yes 1/nps constant random >3.00 0.07 yes yes constant random 2.32 yes yes yes yes yes no number of nonzero history tallies = 20000 efficiency for the nonzero tallies = 1.0000 history number of largest tally = 18954 largest unnormalized history tally = 4.42936E+03 (largest tally)/(average tally) = 4.82162E+02 (largest tally)/(avg nonzero tally)= 4.82162E+02

As shown in the tally fluctuation chart at the end of the output file (above), the mean value, 9.18646E+00, is constant within statistics—or “random”—passing the first of the ten statistical tests. The relative error, E = 0.0589, is <10%, decreasing, and sufficiently decreasing as 1/sqrt(NPS), thus passing the 2nd–4th statistical tests. For next-event estimator tallies, F5 point detectors, the passing criteria is a relative error <5%. The variance of the variance, VOV = 0.0726, is less than 10%, decreasing, and sufficiently decreasing as 1/NPS, passing the 5th–7th statistical tests. For convergence, the VOV test is more important and more sensitive than the relative error test. The Figure of Merit, FOM = 122876, is 1/(TE2), where T is the computation time and E is the relative error. Whereas E2 is proportional to 1/NPS and T is proportional to NPS, the FOM should be a constant. Note that the FOM will vary from one computer to another because T is machine dependent, so FOM can be compared only on the same computer, MCNP version, and installation. Because the FOM is sufficiently constant, the 8th and 9th convergence tests are passed. Only the 10th test—the slope test—fails, as indicated by “no.”

The slope test estimates the rate of large scores that contribute to a tally. If a region of phase space is undersampled, then its score to the tally will be large when it finally is sampled because the Monte Carlo calculation is unbiased and the frequency of large scores times the magnitude of those scores should be the expectation value—or mean score—of the problem.

Figure 2.47 in the output shows how the slope is calculated.

Fig. 2.47
figure 47

Print Table 161 is the unnormed tally density function plot. The abscissa is the number of unnormalized nonzero scores. Following columns are the number of scores, the number density (the number of scores divided by the width of tally bins in column 1), the log of that number, and then asterisks proportionate to the number of scores. The slope fit is denoted by “s”

All scores from 1e−30 to 1e30 are put in 600 logarithmic bins according to score size. Then the first nonzero score bin, in this case 1.26+00, to the highest, in this case 5.01+03, is printed vertically with the number of scores in each bin. The “num den” is the number density, or the number of scores, divided by the size of the bin. Then the logarithm of this density is taken and “*” (asterisks) are printed to represent the “log den” magnitude. The slope is the curve fitted to the 200 highest scores, indicated by the “s.” Print Table 161 shows 6 scores >2.00+03, which flatten the slope to its value of 2.32 and are 200 times the average score of 9.18. The largest unnormalized history tally is 4.42936E+03 (from print table 160), which is 4.82162E+02 times larger than the average tally.

To pass all convergence tests, the slope should be >3 [11]. To increase the slope, the parts of the problem that cause the large scores need to be sampled more often. To determine the cause of the largest score, one can rerun the largest history (the relevant history number is found in print table 160):

NPS 1 RAND HIST=18954 DBCN 2J 1 1

NPS runs a single history. The random number sequence of history 18954 is specified via the RAND card. An event log for the 1st through the 1st history is specified by the DBCN card. Thus, History 18954 is repeated, and from the 1st 50-history Print Table 110 and the event log, the cause of the large score becomes apparent. The cosine of the source in the z-direction, W, is very small and leads to a track length of λ = 4.42936E+03 cm before the particle escapes. The F4 flux tally is F4 = λ/V and V = 1 because of the SD card. In fact, a track length comparable to the disc radius, 1e8, is possible. Thus z-axis cosine bins near zero need to be sampled better. Setting a source bias of

sb1 0 1 20r

causes the SDEF source to start more particles in the small z-angle directions, as seen in Fig. 2.48 from Print Table 10 after the SB bias is added. This source direction bias enables the problem to pass all statistical checks and converge to the correct answer very quickly. (Note the weight multiplier keeps this a fair physical representation of an isotropic source.)

Fig. 2.48
figure 48

The source distribution is shown in optional (requires PRINT card) Print Table 10 of the output file. Each source bin is sampled with the same probability, but because the bin widths are different, the weights vary by nine orders of magnitude. Consequently, the central, very narrow bin—normal to the source direction axis—is sampled nine orders of magnitude more to sufficiently achieve convergence

Thus, the ten statistical tests identify false or poor convergence. The solution to false or poor convergence is variance reduction, which results in sampling the more important parts of the problem more frequently at the expense of less important parts.

What does one do when some tests fail? The ten statistical tests are themselves estimated quantities and can sometimes be misleading:

  • If the first test of the mean having a constant value fails, check the tally fluctuation charts at the end of the output file to see if it is really constant. In rare cases, the MCNP fitting routine fails.

  • If the 2nd test of the relative error >10% (5% for point detectors/next-event estimators) fails, run the problem for more particles or computer time. Note that if the slope <3, the relative error does not exist, and running longer may not help.

  • If the 3rd and 4th tests of the relative error behaving as 1 / SQRT (NPS) fail, in rare cases, the MCNP code is unable to fit the curve properly.

  • If the 5th–7th tests of the variance of the variance, VOV, fail, be aware that VOV does not exist if the slope is <3 and can be ignored in those cases. Otherwise, run longer, or check the tally fluctuation chart for a linearly decreasing VOV.

  • If the 8th and 9th tests of the FOM fail, check the tally fluctuation chart to see if it really is not constant. Note that—in the early stages of a problem—the MCNP code does extra checking for the first 50 to a few thousand histories. Thus, on short runs of merely thousands of histories, the FOM will increase as the checks finish. Also, the FOM is based on computer time so that very short problems of a few seconds often have difficulty passing the FOM tests.

  • If the 10th slope test fails, then the important parts of the problem have probably not been sampled well enough. Note that the slope test requires at least 600 histories and the ability to fit the slope to the 200 largest scores, so a slope of zero is observed for short problems. These short problems need to be run longer for more histories. The user should be aware of what the largest possible score to any tally could be and where large scores are likely. Then, variance reduction should be used to sample whatever causes the largest score and other likely large scores. Note that improper use of variance-reduction methods can also lead to failing the slope test, which often indicates poor convergence right from the start of a problem before the other tests are passed. The slope test is the most challenging test to understand and use to improve convergence.