Output Analysis

A deterministic code either converges or it does not, and you can usually tell which by looking. A Monte Carlo code always produces a number. Run the pin cell for ten histories and MCNP will report a flux with six significant figures. The figures are real; the flux is not. Everything on this page exists to answer one question about a printed tally: is this number an estimate of the answer, or an accident of the particular random walks that happened to be sampled?

MCNP tracks four quantities alongside every mean, and each answers a different part of that question. Taken separately, any one of them can be satisfied by a result that is wrong.

The relative error, printed as error and written R, is the estimated standard deviation of the mean divided by the mean: R = S_x̄ / x̄. It is the headline number, and it is the one most often over-trusted. R measures the spread of the histories that were sampled. It cannot measure the histories that were not — if an important contribution has never been sampled at all, R is small and the mean is wrong, and nothing in R will hint at it.

Because the histories are independent, R falls as 1/√N. That relationship sets the price of precision, and it is steep: cutting R in half costs four times the histories, and a run sitting at R = 0.30 needs nine times its current length to reach 0.10. Whenever a reduction in R looks free, something is wrong with the run, not with the arithmetic.

The variance of the variance, printed as vov, is the estimated relative variance of R itself. It depends on the third and fourth moments of the score distribution, which makes it far more sensitive to a long tail than R is. A rare, enormous score barely moves R and moves the VOV sharply, so the VOV is the earliest warning that the tally is being carried by a handful of histories. It should be below 0.10 and falling as 1/N.

The figure of merit, printed as fom, is FOM = 1/(R²T), where T is the elapsed minutes. Since R² falls as 1/N and T rises as N, the two cancel and the FOM should settle to a constant. That makes it a diagnostic of a different kind: its absolute value measures the efficiency of the calculation and is the right way to compare two variance-reduction schemes on the same problem, while its behavior over the run reports whether the sampling is stable. A FOM that drifts upward or downward is describing a problem whose character is still changing as it runs.

The slope is the one that most directly addresses the histories you have not seen. MCNP fits a Pareto distribution to the largest history scores and reports the exponent of the tail. The Central Limit Theorem — the entire justification for reporting a mean plus an error bar — requires the second moment of the score distribution to exist, and for the fitted tail that requires a slope greater than 3. A slope of 10 is reported as a perfect score, meaning the tail has a definite upper limit. A slope of 0 means MCNP could not fit the tail at all, which needs at least 500 nonzero scores in the bin.

The manual's guidance on R is empirical, and worth memorizing because it is the fastest way to know whether a result deserves any further attention.

Point detectors get a stricter threshold because their scores are dominated by contributions from close to the detector point, which are both the most important and the hardest to sample. At R = 0.10 a point detector may be known only to within a factor of a few. The thresholds also assume the problem is being sampled well everywhere — they are a statement about precision, and they say nothing about whether the model is right.

MCNP automates the reasoning above into ten checks and prints the result after each tally. Before reading them, note the scope, because it is the single most common misunderstanding about MCNP output: the ten checks are applied to the tally fluctuation chart bin only — one bin per tally, chosen with a TFn card, defaulting to the last bin. A tally with fifty energy bins gets ten checks on one of them. MCNP does compare R against the table above for every bin, and summarizes that separately, but the mean-behavior, VOV, FOM, and slope tests exist for the TFC bin alone.

Seven of the ten examine behavior over the last half of the run, which is why a longer run is not merely more precise but better diagnosed. The checks group naturally by the quantity they interrogate. On the mean, one check: no upward or downward trend over the last half. On R, three: an acceptable magnitude, monotonic decrease, and decrease at the 1/√N rate the Central Limit Theorem requires. On the VOV, three of the same shape: magnitude below 0.10, monotonic decrease, and a 1/N rate. On the FOM, two: statistically constant, and free of monotonic trend. On f(x), one: slope above 3.

The checks are nested where nesting is logical. An R that fails the monotonic-decrease test cannot pass the rate test, by definition, and the same holds for the VOV pair — so failures tend to arrive at least two at a time, and the count of failed checks overstates the number of independent problems.

The summary table reports each check three ways: what is desired, what was observed, and whether it passed.

results of 10 statistical checks for the estimated answer for the
tally fluctuation chart (tfc) bin of tally 4

 tfc bin    --mean--    ---------relative error---------   ----variance of the variance----   --figure of merit--  -pdf-
 behavior   behavior     value   decrease   decrease rate   value   decrease   decrease rate    value    behavior   slope

 desired     random      <0.10      yes        1/sqrt(nps)  <0.10      yes         1/nps        constant   random    >3.00
 observed    random       0.01      yes         1.02         0.00      yes         1.00         constant   random    10.00
 passed?      yes          yes      yes          yes          yes       yes         yes            yes       yes      yes

 ==============================================================================================================================

 this tally meets the statistical criteria used to form confidence intervals:
 check the tally fluctuation chart to verify.

Passing all ten is not a guarantee. The manual is explicit about this, and the closing line of the table says so in as many words: the checks improve the probability that a confidence interval covers the true answer, and some as-yet-unsampled part of the problem could still move it. When one or more checks fail, MCNP prints a warning and a page of printed-plot information about f(x), which is there so you can see whether the Pareto fit was reasonable rather than take the slope on faith.

The TFC is the history of the tally rather than its final state, printed at the end of the output file as one row per dump. It is the most informative few lines MCNP produces, and reading it takes about ten seconds once you know the two shapes to look for.

1tally fluctuation charts

                            tally  4
          nps      mean     error   vov  slope    fom

        64000   2.4265E-01 0.0448 0.0021  0.0    1462
       128000   2.4013E-01 0.0316 0.0010  0.0    1465
       256000   2.4108E-01 0.0223 0.0005 10.0    1461
       512000   2.4141E-01 0.0158 0.0002 10.0    1464
      1024000   2.4156E-01 0.0112 0.0001 10.0    1463

This is what convergence looks like. The mean moves within its own error bar and does not drift. The error falls by a factor near 1.41 — that is √2 — every time the history count doubles, which is the 1/√N law in the only place you will ever see it directly. The FOM sits on one value. The slope starts at 0 because there were not yet 500 nonzero scores to fit, then reaches its ceiling.

Two departures matter. The first is a mean that walks steadily in one direction while the error falls, which says the tally is still finding new contributions and the run has not converged no matter how small the error has become. The second is a FOM that steps down each dump, which usually means a rare high-weight event is occasionally entering the tally and the earlier, more optimistic error estimates were never trustworthy. A single jump in any column is not a concern; the manual is clear that small jumps in R, VOV, and FOM are not threatening. A trend is.

For a kcode calculation, this reasoning applies to the tallies but not to k_eff itself, whose convergence is a question about the fission source distribution rather than about tally statistics. That is covered in the criticality example.

For a single answer, read the output file. The lines worth finding are few enough to name: the tally block itself, the statistical-checks table, and the TFC.

# The final tally value and its error
grep -A 6 "^1tally *4" outp

# Did the ten checks pass, for every tally in the deck?
grep -i "statistical criteria" outp

# Any fatal error or warning MCNP wanted you to see
grep -iE "fatal|warning" outp

# Where did the particles go? Losses should be to the outer boundary, not to
# geometry errors or weight cutoffs you did not intend
grep -A 20 "neutron creation" outp

For anything repeated — a parameter sweep, a convergence study, a plot — parse the TFC rather than the tally block. It is fixed-width, has one row per dump, and gives you the whole history in one pass.

import re
import numpy as np
import matplotlib.pyplot as plt


def read_tfc(path, tally):
    """Return the tally fluctuation chart for one tally as an array of
    (nps, mean, error, vov, slope, fom) rows."""
    with open(path) as f:
        lines = f.readlines()

    # Find the TFC block for this tally, then the header row inside it
    start = None
    for i, line in enumerate(lines):
        if re.match(rf"\s*tally\s+{tally}\s*$", line):
            start = i
        elif start is not None and "nps" in line and "fom" in line:
            start = i + 1
            break
    if start is None:
        raise ValueError(f"no tally fluctuation chart for tally {tally}")

    rows = []
    for line in lines[start:]:
        fields = line.split()
        if len(fields) != 6:
            if rows:          # the block has ended
                break
            continue          # still in the blank lines under the header
        try:
            rows.append([float(v) for v in fields])
        except ValueError:
            break

    return np.array(rows)


tfc = read_tfc("outp", tally=4)
nps, mean, error = tfc[:, 0], tfc[:, 1], tfc[:, 2]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

# Is the mean settled? Plot it with its own error bar, not on its own.
ax1.errorbar(nps, mean, yerr=mean * error, marker="o")
ax1.set(xscale="log", xlabel="histories", ylabel="tally 4 mean")

# Is the error falling at the required rate? On log-log axes, 1/sqrt(N) is a
# straight line of slope -1/2, so a reference line makes the answer obvious.
ax2.loglog(nps, error, marker="o", label="observed")
ax2.loglog(nps, error[0] * (nps[0] / nps) ** 0.5, "--", label=r"$1/\sqrt{N}$")
ax2.set(xlabel="histories", ylabel="relative error")
ax2.legend()

fig.tight_layout()
fig.savefig("convergence.png", dpi=150)

Plotting the error against a 1/√N reference line is worth the extra three lines. The statistical checks answer the same question, but a run that is drifting away from that line shows it here several dumps before check 4 fails.

Mesh tallies are a different problem: meshtal is a separate file, and its layout depends on the mesh you asked for, so parsing it by hand is rarely worth the afternoon. Established tools already read it — PyNE and the mc-tools scripts both convert meshtal to formats that ParaView and VisIt can render as a volume. Remember that an FMESH without an FM card scores flux, not power, and that a multiplier on a mesh applies one material everywhere, so it cannot fix that.

A failed run poses one question: is the problem that too few histories were run, that the histories were spent in the wrong places, or that the model does not describe the thing you meant? The three have completely different remedies, and the FOM is what separates them.

If the FOM is constant and the error is simply larger than you need, the calculation is healthy and merely short. More histories will fix it, and the 1/√N law tells you exactly how many: the ratio of errors, squared. This is the only one of the three cases where a longer run is the right answer, and it is the one people reach for in the other two.

If the FOM is very low and constant, the run is stable but inefficient — the sampling is spending its time where the tally is not. Ten times the histories buys a factor of about three in error, which is usually unaffordable. This is what variance reduction is for: importance weighting, weight windows, or a source biased toward the part of phase space that reaches the tally. Judge whether it worked by comparing the FOM before and after, not the error, since the error will improve simply because the run took longer.

If the FOM drifts or the slope stays below 3, neither remedy applies yet. The tally is being carried by rare high-weight histories, and a longer run mostly buys more chances to be surprised. Find out which histories are scoring so heavily. Aggressive splitting or a weight window whose bounds do not match the physical attenuation both produce this signature, and in that case the variance reduction is the cause rather than the cure.

And if the answer is precise but physically implausible, the statistics are not the subject. Check the geometry with a plot, check that every cell has the material and density you intended and the sign convention you intended, check the source definition, and check the particle-loss summary for losses to anything other than the outer boundary. Troubleshooting works through the specific failure modes.

A tally value is not a result on its own. It is a result together with the model that produced it, and reproducing it a year later requires the same MCNP version, the same cross-section libraries, the same input, and the same statistical caveats. Version and library matter more than people expect: the same deck run against .70c and .80c data will not give the same k_eff, and the difference can exceed the reported uncertainty.

What is worth writing down, then, is whatever a reader would have to guess: which MCNP version and build, which library each material resolved to, which assumptions and simplifications went into the geometry, the statistical quality of the tallies rather than just their means, and any warning MCNP printed that you decided to accept. The input file and the analysis script belong in the record too — a plot that cannot be regenerated is a claim, not evidence.