- Email: [email protected]
A generalized approach to designing specific-purpose Monte Carlo programs: a powerful simulation tool for radiation applications
Nuclear Instruments and Methods in Physics Research A299 (1990) 1-9 North-Holland
1
Section I. Detectors
A generalized approach to designing specific-purpose Monte Carlo programs : a powerful simulation tool for radiation applications Kuruvilla Verghese and Robin P. Gardner Center for Engineering Applications of Radioisotopes, North Carolina State University, Raleigh, NC 27695-7909, USA
Well-designed Monte Carlo simulation programs provide the most accurate means of predicting the response of radiation measurement systems. Such simulations can replace large sets of experiments for the purpose of design optimization, calibration and analysis of sensitivity to measurement interferences. The use of detector response functions can extend Monte Carlo simulation to systems involving spectral responses. Compared to the use of general-purpose Monte Carlo codes, our work on developing specific-purpose codes has shown various advantages : greater "user-friendliness", use of an optimum set of variance reduction methods, superior code performance in most cases, and better understanding of the basic physics of the problem. The specific-purpose codes are relatively easy to develop using modules that we have developed for geometry, cross sections, detector response functions, standard as well as new and/or improved variance reduction schemes, etc. The measurement applications for which we have designed and used specific-purpose Monte Carlo simulation codes in recent years include on-line capture gamma-ray analysis of coal ; X-ray fluorescence analyzers; nuclear tools for oil-well logging by neutron moderation, pulsed neutron lifetime analysis and by gamma-ray backscatter; and industrial tomograpluc scanners. 1. Introduction Systems for radioisotope measurement applications are commonly designed by converting a measurement principle into a prototype measurement system, carrying out design calculations based on simplified models in order to approximately optimize the design, verifying the adequacy of the prototype design using an extensive set of laboratory measurements, refining the prototype system for commercial application, and calibrating the commercial device using sets of well-characterized calibration standards. These devices range in complexity of design and of the measurement principle all the way from relatively simple nuclear gauges to radiation imaging systems such as tomography scanners . As the complexity of the measurement application increases, it becomes essential to increase the sophistication of the design calculation model . One goal of this paper is to convincingly show that Monte Carlo simulation models present the best option for design calculations and that they are often sufficiently accurate and versatile for sensitivity evaluations and even for the measurement system calibration . With the recent revolution in computer capabilities, Monte Carlo simulation has become feasible without requiring access to a large-scale computational facility . It is possible to develop Monte Carlo programs which can run on relatively small computers and can simulate in reasonable time very complex geometries using detailed treatment of radiation-interaction physics and
microscopic data. What is required, however, is to tailor the program for the specific type of application . General-purpose codes such as MORSE [11 and MCNP [2] cover too wide a range of applications to be efficient and easy to use for radioisotope measurement applications. Another goal of this publication is to encourage wider use of Monte Carlo simulation within the radiation measurement community through a discussion of examples of successful specific-purpose codes which the authors and their group have developed over the last several years. Most of our Monte Carlo studies have been on systems that use neutrons, gamma-rays, or their combinations. Charged-particle Monte Carlo methods and examples will not be discussed here . In 1983, the authors published a paper [3] reviewing the development of specific-purpose Monte Carlo codes for radiation gauges and analyzers. It contrasted our experiences with general-purpose neutron-photon codes against a few of our specific-purpose codes (neutron moisture gauge, gamma-ray density gauge, X-ray fluorescence analyzers and a neutron-capture gamma-ray analyzer) and concluded that the latter are cost-effective and accurate enough to provide a nonexperimental method for calibration . Since that time, our applications have been extended to include nuclear geophyics problems, computed tomography, and the inverse problem (determination of elemental amounts from pulse-height spectra) of nuclear analyzers . The scheme for developing Monte Carlo programs has been generalized through a modular approach, and several improved methods of
O16ß-9002/90/$03 .50 © 1990 - Elsevier Science Publishers B V. (North-Holland)
I. DETECTORS
2
K Verghese, R. P. Gardner / Designing specific-purpose Monte Carloprograms
variance reduction have been formulated and implemented for deep-penetration neutron and photon problems . Specific-purpose Monte Carlo programs for radiation measurement applications are optimally formulated in a modular manner using subroutines which are common to many of the systems. Basically, these modules cover the following aspects of Monte Carlo: Generalized geometry and tracking of particle paths through the various zones, particle-collision physics, microscopic data pertinent to the radiation interactions, response functions of detectors, and variance reduction methods. Our approach to the development of these modules is briefly discussed. It is hoped that this discussion will help to remove the mystery of Monte Carlo techniques which has been perpetuated through the use of large general-purpose codes. Because of space constraints, it will be assumed that the reader is familiar with the fundamental methods of Monte Carlo for radiation transport simulation as outlined m refs . [4] and [5].
2. Generalized approach to specific-purpose Monte Carlo programs A Monte Carlo program for radiation measurement applications basically must have the following features : (1) it must be flexible enough for the user to easily make changes in the geometry of the measurement system and to track paths of particles through the various zones ; (2) the physics of all of the relevant radiation interactions must be included and sampling schemes for the interaction probability density functions must be sampled efficiently ; the microscopic physics data, such as cross sections, secondary-particle yields, etc., must be stored and accessed by the program in the most efficient manner ; (4) the response of the radiation detector (such as efficiencies, pulse-height spectral response, etc.) must be included and used in the most efficient manner ; and the variance reduction methods must be optimum for the type of measurement system, and it must be easy for the user to specify the parameters, if any, of the variance reduction schemes . It is in the last two items that specific-purpose codes differ significantly from the general-purpose codes. Spectral responses of detectors are not easily incorporated into general-purpose codes. Also, the user of a general-purpose code is usually offered a menu of variance reduction methods to choose from with the recommended method being the concept of splitting and
Russian Roulette [5] based on user-supplied importance regions and associated importance values. This method is generally not the optimum for radiation measurement problems and requires considerable experience to set up . In contrast, specific-purpose codes can be hard-wired with the optimum set of variance reduction methods for the type of problem. We will now proceed to discuss our approach to each of the above features. 2.1 . General geometric modelling
Radiation measurement system geometry can get quite complex m many applications and it is essential to have a general geometry modelling and particle-tracking program. T.H . Prettyman of our group has developed such a program called HERMETOR [6] which is being used to some of our Monte Carlo codes. HERMETOR allows the user to describe heterogeneous systems made up of homogeneous zones whose boundaries can be described by second-order surfaces . It uses a combination of boundary representation to construct zones and constructive solid geometry logic to combine the zone surfaces to form the desired geometric system . The input file for HERMETOR is simple and a 2D plotting package is used to verify the accuracy of the input. HERMETOR has proven invaluable in modelling complex systems such as nuclear well-logging tools [6] and tomographic scanners [6] which contain multiple zones of shielding, collimators and detectors. 2.2 . Neutron and photon interaction physics
The probability distributions and the appropriate methods for sampling these distributions for neutronphoton Monte Carlo simulations are well-documented in the literature (for example, ref. [5]) . In most radiation measurement applications, the source can be assumed to be a point emitting isotropically and the source energy spectra are well known either as empirical for252Cf mulae (for example, the fission spectrum) or as pointwtse tabulations as in the case of most (a, n) neutron sources . The major interactions for neutrons are elastic- and inelastic scattering and neutron-capture reactions. Secondary neutron emission reactions such as fission are rarely encountered. The angular distribution of elastically scattered fast neutrons can be treated in many cases with sufficient accuracy using the isotropic scattering assumption in the center-of-mass system . Thermal-neutron scattering is often treated adequately by the free-gas model . However, one word of caution in this regard is that the degree of sophistication required in the physics treatment is problem-dependent and if model results show differences from good experimental
3
K. Verghese, R P. Gardner / Designing specific-purpose Monte Carlo programs
benchmarks, it is prudent to consider adding a more detailed treatment of the collision physics . The physics of gamma-ray interactions are more straightforward than that of neutron interactions and are adequately treated in textbooks [7] . The major angular distribution functions that are encountered are Rayleigh scattering and Compton scattering with the classical angular distributions corrected by form factors at very low energies. The most common secondary particles emitted are electrons which, because of their relatively short range, need not be tracked in most cases .
2.3. Interaction-physics data The major source for neutron-physics data in the USA is the Evaluated Nuclear Data Files/Part B (ENDF/B) compilation at Brookhaven National Laboratory, and photon-physics data are available from the National Institute of Standards and Technology . These data banks are invaluable sources of pointwise data . Generally, the data sets are very large, and particular attention must be paid to reducing the data storage and the accessing times since the data sets might have to be accessed millions of times during the course of a Monte Carlo run . One way of compressing the data sets is to develop empirical functions that fit the pointwise data. For photon cross sections, this is feasible to do since the data generally behave smoothly . In our work, we use cubic-spline-fitted Rayleigh, Compton, photoelectric and pair-production cross sections. The fitting algorithms for continuous and piecewise continuous functions are discussed by Yacout et al . [8] . Also, for lowenergy photons, we have developed [9] and used two-dimensional cubic splines to fit cumulative probability density functions of coherent and incoherent scattering . The fits are accurate to within the error in the tabulated pointwise data and the speed-up in data access is of the order of 80 compared with binary interpolation of pointwise data . Neutron cross sections are too bumpy for fitting and are treated as pointwise binary files . Some access time can be saved by preprocessing the cross sections such that linear-linear interpolation can be performed between the data points (rather than log-linear or log-log interpolation) . Neutron cross-section files are generally the limiting aspect of computer memory requirements in applications that involve a large number of elements . The other large neutron data sets are capture gamma-ray energies and yields, nuclear-energy levels and secondary gamma-ray yields for inelastic scattering, anisotropicscattering angular distributions as a function of neutron energy, and thermal-neutron scattering data for the S(a,,8) treatment.
2.4. Detector response For all radiation measurement applications, the final tally in a Monte Carlo calculation is the score in the detector(s) . This score in many cases may be multidimensional . If only the counting rates are required, the mean score of particles interacting in the detector is multiplied by the source emission rate. However, if the counting rate as a function of energy deposition in the detector is required (as in energy spectroscopy), the entering particle must be followed through its interactions inside and outside the detector because of the potential to leave the detector and at a later time return to the detector . When applicable, a more efficient scheme is to devise and store the spectral response function of the detector as a probability-density function and to convolve the weight distribution of the particles interacting in the detector as a function of their incident energy with the response function . The spectral response function R(E, E') dE is the probability that radiation entering the detector with energy E' will deposit energy between E and E + d E in the detector . For detectors with low aspect ratio or those in applications where the radiation enters the detector only from one side, the response function can be deduced as a semi-empirical analytical function from experimental spectra for known source energies . Examples are provided in refs . [10] and [11] for a Ge detector and in refs . [12-14] for a Si(Li) detector . Jin et al . [15] provide an example of the use of the response function in conjunction with a Monte Carlo program to generate the pulseheight spectrum from a coal elemental analyzer and Yacout et al . [16] used the response function of a Si(Li) detector to calculate the spectral response of an X-ray fluorescence analyzer . In both cases, the simulated spectra matched up very well with spectra from experimental standards . An improved and more general Si(Li) response function for X-rays covering the range 5-60 keV is discussed by He et al . [17] .
2.5 . Variance reduction methods Many of the variance reduction methods that are useful in specific-purpose codes are discussed by Carter and Cashwell [5] . They generally fall within the category of biased sampling methods where instead of sampling from the actual probability-density function one samples from a fictitious density function which is known to enhance the chance of success (detection) for that particle . The bias introduced by the fictitious function is corrected by using a weight factor to alter the weight of the particle . In many cases the biasing scheme is physically obvious . The tricky problems are those of deep penetration where the particles have to get through a large number of mean free paths in order to reach the t . DETECTORS
K. Verghese, R. P. Gardner / Designing specific-purpose Monte Carlo programs
4
detector and narrowly collimated
systems where the
particles can be detected only by travelling in the "right" direction with the proper energy and weight . Unfor-
tunately many of the radiation measurement application problems fall in these categories. Examples are
nuclear well-logging systems [18] (such as the neutron porosity tool and the lithodensity tool) and photon tomography systems. The problem is that in order to
enhance the chance of success one must preferentially
sample those trajectories in the phase space that are
splitting weight . The method is a generalization of the DXTRAN scheme found in the MCNP code [2] manual and the details of implementation for photons are dis-
cussed by Prettyman et al . [23] .
The last of the powerful variance reduction schemes
that must be discussed is the concept of correlated
sampling . Correlated sampling is a scheme for studying
the sensitivity of the Monte Carlo result to small per-
turbations in the physical parameters of the problem. For example, it may be used to determine the change in
that one can a priori identify such trajectories . That is
the system response caused by small changes in the particle cross sections, density, elemental composition,
calculations or simulations. Assuming that such path-
accomplished without rerunning the Monte Carlo code
techniques to reduce the variance in deep-penetration
correlated sampling is not new, it has been implemented differently by different investigators . An efficient correlated sampling scheme that the authors have imple-
likely to contribute to the detector score and it assumes often not a simple task and must be deduced from prior
ways can be identified, we have successfully used two problems : one is path-length stretching or exponential transform [5] and the other is direction biasing [19] . In
etc., in specified regions of the phase space. This can be
for all of the perturbed samples. While the concept of
stretched using a fictitious cross section if the particle is
mented [24] consists of doing the Monte Carlo simulation in the usual way for the reference sample . Simulta-
shortened if the particle is travelling in the opposite
the perturbed samples are also calculated assuming that
path-length
travelling in
stretching,
the
distance
to
collision
is
the right direction and the distance is
direction. The net effect is that more particles are able
to get to the target (detector) . In direction biasing, the particle direction is biased towards the target following
neously during this simulation, the responses for all of
the same particle path is valid in both the reference sample and in the perturbed samples. With this assump-
a scattering interaction . A combination of the two
tion which is valid for small perturbations the only additional calculation needed for the perturbed sample
deep-penetration problems . Implementation of pathlength stretching is straightforward and requires one
weight for the reference sample is to be calculated . Thus the amount of extra calculation for a large set of per-
methods can yield significant reduction in variance in
input parameter while the implementation of direction biasing for the common neutron and photon scattering
response is the modified weight at each event where the
turbed samples (say, 50 samples) only adds a few percent to the total computer time but the results for the
kernels required considerable development work, as dis-
perturbed samples are highly correlated to the result for
Another variance reduction method which we have used in several of our simulations is an estimation
lated sampling provides an excellent way to calibrate
cussed by Gardner et al [20] .
method known as statistical estimation [21] . It is a way of estimating the score to the detector if a particle were to scatter from any point of interest towards the detector, reach the detector without further collisions along
the reference sample.
For nuclear gauge and analyzer applications, corre-
the system without having to resort to a large set of
well-characterized experimental
standard
samples. A
very small number of experiments might be necessary in
order to tie down the response of the reference sample,
the way and interact in the detector . Statistical estimation for amsotropic scattering can be a time-consuming
but the additional data points for calibration can be generated using the responses for perturbed samples.
mation .
simulation is
dure for variance reduction is possible. This scheme, which is called a combined splitting-biasing game, pre-
are discussed by Yuan et al . [26] .
calculation and Gardner et al . [21] and Mickael et al. [22] have discussed efficient ways to do statistical estiFor measurement systems with narrowly collimated geometry as in tomographic systems, a different proce-
The use of correlated sampling for density- and/or elemental-composition changes in nuclear-logging tool
discussed by Gardner et al . [25] . The
details of implementation for a nuclear coal analyzer to calculate the response as a function of coal composition Generally, the effect of changes m geometric dimen-
determines the appropriate types of collision paths that could potentially make a contribution to the detector
sions of the measurement system cannot be investigated
carrying a weight equal to the probability of achieving the desired path for that type . The nice feature of this
3. Specific-purpose Monte Carlo applications
score, and the original particle is split into multiple particles (one for each type) with each split particle
scheme is that if the split particle is sampled from an arbitrary distribution, one does not have to know the
by this method .
Three areas of applications of specific-purpose Monte Carlo simulation to nucleonic measurement systems are
K. Verghese, R . P. Gardner / Designing specific-purpose Monte Carlo programs
discussed. The versatility and accuracy of Monte Carlo simulations rather than the mechanics of implementation will be the focus of the discussion since the details of implementation can be found in the papers and student dissertations that are referenced . 3.1 . Nuclear analyzers for elemental composition We have simulated the response of two types of elemental analyzers: a neutron-capture prompt gammaray analyzer (NCPGA) for on-line composition measurement and a radioisotope source-excited X-ray fluorescence analyzer (XRF). 3.1 .1 . Neutron-capture prompt gamma-ray analyzer A NCPGA is used commercially for coal analysis primarily for the purpose of ash and sulfur measurements. The system that we have studied is similar in hardware to the CONAC system [27] developed by SAIC under the sponsorship of the Electric Power Re252Cf search Institute (EPRI) . It uses a source of fission neutrons and a Ge detector in transmission geometry . The neutrons are moderated and allowed to enter a 91 cm x 91 cm x 25 cm thick box of coal . The capture gamma-ray spectrum resulting from the absorption of
5
thermal neutrons in the coal elements is analyzed to determine the elemental concentrations . Fig. 1 schematically shows the simulated system . The measurement method is complicated by the fact that the response for any particular element is, in principle, nonlinearly dependent on the amounts of all of the other elements . As a result, accurate calibration of such a device as CONAC requires a very large set of accurately characterized, large-volume experimental standards in order to determine the interelement effects. Monte Carlo simulation combined with correlated sampling can be a powerful tool to determine the interelement effects and thereby simplify the calibration procedure. The details of the simulation procedure are given by Jin et al. [28] . Fig. 2 shows a recent comparison of the Monte Carlo simulated spectrum to the experimental spectrum for the same coal composition . Considering that the Ge detector response function that was used in the simulation is only valid between 1 .5 and 6 MeV, the uncertainties in the coal composition, and that the counting rates range over six decades, the agreement is quite good and adequate for calibration applications . It is important to note that there is no normalization factor used to normalize the two spectra.
Fig. 1. Schematic configuration of a laboratory simulation system for the elemental analysis of coal : (1) polypropylene, (2) paraffin work, (3) paraffin and lithium carbonate, (4) lead, (5) polyethylene, (6) coal, (7) paraffin wax and lithium carbonate, (8) Ge detector, ( " ) neutron source . All dimensions are in cm . 1. DETECTORS
K. herghese, R. P. Gardner / Designing specific-purpose Monte Carlo program ,
10
6
10
5
10
Uô
3
10
2
10
6 E
7
I
I
I
8
9
10
(MeV)
Fig. 2. Comparison of the measured gamma-ray pulse-height spectrum of a coal sample to the spectrum generated by Monte Carlo simulation for the same coal composition .
3.1 .2 . X-Ray fluorescence analyzer XRF analyzer systems are commonly used in analytical laboratories to do elemental analysis with a sensitivity in the range of parts per million for all but the very
light elements . Although
the spectral-response
simulation techniques are equally well applicable to
radioisotope-source-excited systems as well as to X-ray tube-excited analyzers, we have chosen to demonstrate
the versatility of Monte Carlo simulation on a sourceexcited system that is available for experimental validation in our laboratory .
Fig. 3 shows an XRF spectrum simulated by our
specific-purpose Monte Carlo code for austenitic stainless steel containing six elements and excited by a 109Cd
source and compared with the experimental spectrum .
Once again there is no normalization performed between the simulated and experimental spectra. The agreement is excellent in the fluorescent region . The
disagreement in the source backscatter region is because the Si(Li) detector response function that was used was valid only at lower energies than the backscattered source energies. The new response function by He et al . [17] which was referred to in section 2.3 is valid up to 60 keV and should perform better in the comparison.
So far, what we have demonstrated is that Monte Carlo simulation using fundamental data can be used to predict the spectral response of analyzer systems where the response is nonlinearly dependent on the elemental concentrations . This is referred to as the "forward"
problem m contrast to the "inverse" problem which is the problem of obtaining elemental concentrations from
measured spectral response of the system . The classical way to decompose a composite spectrum is by using the library least-squares (LLS) method, m which one re-
quires the shapes of spectra from each of the elements that are present in the composite. Unfortunately, in a
nonlinear system these shapes are, to say the least, difficult to obtain experimentally since, in principle, the
intensity and shape of the spectrum from each element will depend on the amounts of the other materials that are present in the sample matrix . However, they are easily obtained by Monte Carlo simulation by simply tallying the contribution to the spectral response from
each element separately . Thus a combination of the Monte Carlo and LLS methods may yield a way of solving the inverse problem. We are starting to apply I
o'
5
7
9
11l
I 13 15 17 ENERGY 1KEV1
19
21
23
25
27
Fig. 3 . Comparison of the measured and simulated pulse-height spectra of the backscattered 109Cd photons and the fluorescent X-rays from a 304 stainless steel sample containing Si, Cr, Mn, Fe, Ni and Mo (- Monte Carlo simulation, - - - experimental).
the
Monte
Carlo-Library
least-squares
(MCLLS)
method to both the coal analysis problem and to the
XRF analysis problem. The MCLLS procedure consists of generating the library spectra from Monte Carlo simulation based on an initial estimate of the sample composition, perform the LLS analysis using these
K Verghese, R.P Gardner / Designing specific-purpose Monte Carloprograms
library spectra to calculate the elemental composition of the sample and if the new composition is quite different from the assumed composition, alter the assumed composition and repeat the procedure . The preliminary results of MCLLS for the XRF problem appears promising [29] and the work on applying it to the coal analysis problem is continuing. If successful, the MCLLS method offers the potential for eliminating expensive and tedious experimental calibrations. 3 .2 . Nuclear-tool responses for oil-well logging Nuclear-measurement methods are commonly used in logging bore holes for petroleum exploration and production . Such measurements routinely provide information on porosity of the formation, salinity and percent saturation of pores, formation density, effective atomic number, etc ., and are used in conjunction with other measurements to deduce the necessary geophysical data . The tools used for these measurements consist of a radiation source (neutrons or gamma-rays) and a set of two or more detectors . The detector responses are analyzed to determine the geophysical quantity of interest as a function of depth down the bore hole . The tools are generally run along the side of the bore hole . There are three nuclear tools for which Monte Carlo codes are being developed . These are the dual-spaced neutron-porosity tool, the pulsed-neutron lifetime tool, and the lithodensity gamma-ray tool . 3.2.1 . Dual-spaced neutron-porosity tool The dual-spaced neutron-porosity tool is a 7 .6 cm diameter tool that contains a few-Ci Am-Be neutron source and two 3 He neutron detectors mounted axially at two different distances from the source . In addition, such tools contain various shield materials placed strategically to enhance the porosity sensitivity . The ratio of the counting rate on the near-spaced detector to that on the far-spaced detector correlates well with the porosity of the formation and it is nearly independent of the bore hole conditions . For Monte Carlo simulation, this system forms a deep-penetration problem since the far-spaced detector is placed generally at around 50 cm from the source and hence requires extensive use of variance reduction schemes. General-purpose codes such as MCNP [2] do this by splitting up the geometric regions of the problem into a grid and neutrons which enter a grid region are either split into multiple neutrons with reduced weight or killed by Russian Roulette based on how important that region is to the final score. The set of importance values for the grid must be provided by the user . For further variance reduction, other schemes such as weight windows are used to avoid excessive weight fluctuations for the neutrons which are detected . In contrast we have preferred the combined use of path length stretching,
7
direction biasing, and statistical estimation in our code, McDNL [30] . This makes McDNL considerably easier to use than MCNP for this problem . Correlated sampling is an integral part of this code and has proved invaluable for calculating environmental corrections in porosity logging [31] . The McDNL code has been in use in the petroleum industry now for about two years . Although the previously reported computing times [30] are comparable with that of MCNP, recent improvements make it significantly faster than MCNP . A major attractive feature of McDNL is its user-friendliness. 3.2 .2 . Pulsed-neutron lifetime tool The objective of the pulsed-neutron lifetime measurement is to determine the macroscopic absorption cross section of the formation . The tool is a 4 .3 cm diameter tube which houses a (13,T) neutron generator and two gamma-ray detectors in addition to shielding configurations . A 14-MeV neutron burst lasting for several microseconds is emitted from the generator . The neutrons moderate and are absorbed generating gamma-rays from inelastic neutron scattering and radiative capture. These gamma-rays are detected as a function of time measured from the neutron burst . This temporal spectrum is analyzed to yield the absorption cross section of the formation. Certain pulsed-neutron tools also can record the gamma-ray pulse-height spectra from which elemental analysis for the formation can be performed . For Monte Carlo simulation, this is a combined neutron/photon problem requiring detailed physics treatment because of the high neutron energy and the importance of inelastic neutron scattering. Also, the fact that responses in up to 100 time channels must be tallied with good statistical accuracy complicates the simulation somewhat . However, because of gamma-ray detection rather than thermal neutrons, path-length stretching and direction biasing have not been necessary as variance reduction methods. Statistical estimation for gamma-rays is the major variance reduction method . Also, as in the other examples, correlated sampling provides a useful addition to study environmental effects on log response . McPNL, our Monte Carlo code for the pulsed-neutron logging tool has been benchmarked against a set of laboratory test pit data and the simulated results agreed exceptionally well with the experimental results [32,33] . 3.2 .3 . Gamma-ray hthodenstty tool This tool is a gamma-ray backscatter tool for the determination of formation density and to provide a measure of the lithology . It consists of a gamma-ray source (typically 135 Cs) collimated to emit directly into a cone within the formation and contains two collimated gamma-ray detectors . The pulse-height spectrum of I. DETECTORS
8
K. Verghese, R P Gardner / Designing specific-purpose Monte Carloprograms
scattered photons is recorded and used to differentiate between different formation lithologies. The complex part of simulating this tool is the low yield of counting rate per source photon (of the order of 10 -s ) resulting from the narrow collimation. A code called McLDL is currently under development and we expect to have the first version ready by the end of the year. 3.3 . Tomographic scanners for industrial tomography
Industrial tomographic systems are becoming widespread and diverse in design from the tomographic systems used in radiology. These systems are complicated and expensive to design and build. The complicating features are stringent collimation, multiple detectors and generally complex geometry. Simulation will play a decisive role in their design development and optimization. Accurate simulated data can be used also to develop image reconstruction methods. The power of Monte Carlo for tomographic simulation has been long recognized and used [34] to some extent for transmission-tomography systems. We are investigating the design of a system that would use both transmitted and scattered photon detector signals to image the atomic number as well as the density distributions in industrial samples. Monte Carlo simulation verified by relatively simple experiments in the laboratory is the approach taken. A code called MCPT [23] has been written for this purpose and is being benchmarked using scattered photon pulse-height spectra from inhomogeneous phantoms . The code incorporates the geometric modelling package, HERMETOR, and a general combined splitting-biasing scheme to send photons to detectors without scattering and after single, double and higher-order scatters in the sample as well as in the other important zones of the scanner. The simulated responses will be used to develop reconstruction methodology from the combined transmission and scattered pulse-height spectra. 4. Discussion and conclusions We have reviewed the past seven years of our work on the development and uses of specific-purpose Monte Carlo programs for a variety of radiation-measurement apphcations . All of our programs are developed on the relatively small Vax system of computers including the MicroVax 11 and DecStation 3100, but they are written in standard Fortran and are portable to more powerful systems such as Cray systems. Accurate CPU times for the various applications and the relative standard deviations of the simulated responses are given in the individual references that are cited, but generally the run times are of the order of one to a few hours on the DecStation
which runs approximately a factor of 15-18 times faster than the MicroVax II . In the early days, each of our specific-purpose codes was written independently of each other, but over the years we have generated a large set of subroutines that are common to many applications . The methodologies used in these subroutines are reviewed in section 2 of this paper and for details the reader is again referred to the set of references that are cited in that section. For a completely new application, typically the time it takes to assemble a working code is no more than the time to obtain a general-purpose code, install it and learn to run it in the optimum manner . But the major advantages of the specific-purpose codes are their ease of use and the fact that by writing the code one learns a great deal more about that application than if one just learns to set up the input to a general-purpose code. We have chosen to illustrate the versatility of specific-purpose Monte Carlo by discussing a wide range of fairly complex problems in radiation applications . Many of the applications have necessitated the development of new or improved implementation methods for geometric modelling, data fitting, detector response functions and variance reduction. Programs for some of the applications are still under development and improved versions of the others such as McDNL and McPNL are also being developed. Acknowledgements The authors acknowledge the interest, encouragement and financial support provided by the Amoco Production Company, Mobil Research and Development Corporation, BP Research, Arco Oil & Gas Company, Exxon Production Research Company, Atlas Wireline Services, Teleco Oilfield Services Inc ., and Conoco Inc., under an Associates Program supporting Nuclear Oil Well Logging Research and BP America for support under a separate contract for Industrial Tomography research . References [I] E.A . Straker, W H. Scott, Jr, and N.R. Byrn, Radiation Shielding Information Center Report CCC-127D (Oak Ridge National Laboratory, 1972) (available from National Technical Information Center, Springfield, VA, USA) . [2] Los Alamos Monte Carlo Group, Los Alamos National Laboratory Report LA-7396-M(rev .) (Los Alamos National Laboratory, April 1981). [3] R.P . Gardner and K. Verghese, Isotopenpraxis 19 (1983) 329. [4] W.E. Selph and C.W . Garrett, m. Reactor Shielding for Nuclear Engineers, ed . N.M . Schaeffer (U .S . Atomic Energy Commission, TID-25951, 1973) p. 207.
K. Verghese, R. P. Gardner / Designing specific-purpose Monte Carlo programs [5] L.L . Carter and E.D . Cashwell, Particle-Transport Simulation with the Monte Carlo Method, ERDA Critical Review Series TID-26607 (Technical Information Center, Springfield, VA, 1975). [6] R.P. Gardner, T.X . Prettyman, M. Mickael and K. Verghese, Trans. Am . Nucl . Soc. 57 (1988) 115. [7] R.D . Evans, The Atomic Nucleus (McGraw-Hill, New
York, 1955). [8] A.M . Yacout, R.P . Gardner and K. Verghese, Nucl . Instr. and Meth. 220 (1984) 461. [9] A.M. Yacout, R.P . Gardner and K. Verghese, X-Ray Spear. 15 (1986) 256. [10] Y. An, R.P . Gardner and K. Verghese, Nucl. Instr. and Meth . A242 (1986) 416. [11] M.C . Lee, K. Verghese and R.P . Gardner, Nucl . Instr. and Meth . A262 (1987) . [12] R.P. Gardner and J.M. Doster, Nucl . Instr. and Meth. 198 (1982) 381. [13] A.M. Yacout, R.P. Gardner and K. Verghese, Nucl. Instr. and Meth . A249 (1986) 121. [14] R.P. Gardner, A.M . Yacout, J. Zhang and K. Verghese, Nucl. Instr. and Meth . A242 (1986) 399. [15] Y. Jin, R.P . Gardner and K. Verghese, Nucl. Geophys. 1 (1987) 167. [16] A.M . Yacout, R.P . Gardner and K. Verghese, Adv. X-Ray Anal. 90 (1987) 121. [17] T. He, R.P. Gardner and K. Verghese, these Proceedings (7th Symp. on X- and Gamma-Ray Sources and Applications, Ann Arbor, MI, 1990), Nucl. Instr. and Meth . A299 (1990) 354. [18] D.V . Ellis, Well Logging for Earth Scientists (Elsevier Science Publ. Co., Amsterdam, 1987) chs. 9, 12. [19] J.S . Hendricks and L.L . Carter, Nucl . Sci. Eng. 89 (1985) 118. [20] R.P . Gardner, M. Mickael and K. Verghese, Nucl . Sci. Eng. 98 (1988) 51 .
9
[21] R.P. Gardner, H.K. Choi, M. Mickael and K. Verghese, Nucl. Sci. Eng. 95 (1987) 245 [22] M. Mickael, R.P. Gardner and K. Verghese, Nucl . Sci. Eng. 99 (1988) 251 . [23] T.H . Prettyman, R.P . Gardner and K. Verghese, these Proceedings (7th Symp . on X- and Gamma-Ray Sources and Applications, Ann Arbor, MI, 1990), Nucl . Instr. and Meth . A299 (1990) 516. [24] T.C . Clark, R.P . Gardner and K. Verghese, Nucl . Instr. and Meth . 193 (1982) 365. [25] R.P . Gardner, M.W . Mickael, C.W. Towsley, C.M . Shyu and K. Verghese, IEEE Trans . Nucl . Sci., in press. [26] Y.L . Yuan, R.P. Gardner and K. Verghese, Nucl. Technol.
77 (1987) 97 . [27] G. Reynolds, H. Bozorgmanesh, E. Ehas, T. Gozani, T. Maung and V. Orphan, EPRI Report EPRI-989, vol. 1 (1979) . [28] Y. Jm, unpublished Ph .D . Dissertation, North Carolina State University (1986) (available from University Microfilms, Ann Arbor, MI). [29] K. Verghese, M. Mickael, T. He and R.P . Gardner, Adv. X-Ray Anal . 31 (1988) 461. [30] M. Mickael, R.P . Gardner and K. Verghese, Proc . SPWLA 29th Ann. Logging Symp ., Paper MM (Society of Professional Well Log Analysts, 1988). [31] M. Mickael, H.K . Choi, R.P . Gardner and K. Verghese, Nucl . Geophys. 3 (1989) 339. [32] H.K . Choi, K. Verghese and R.P . Gardner, Nucl. Geophys . 1 (1987) 71 . [33] H.K . Choi, K. Verghese and R.P . Gardner, Nucl. Geophys . 3 (1989) 167. [34] G.T . Herman and S. Rowland, Technical Report no. 130 (State University of New York, Buffalo, 1978).
I. DETECTORS
1
Section I. Detectors
A generalized approach to designing specific-purpose Monte Carlo programs : a powerful simulation tool for radiation applications Kuruvilla Verghese and Robin P. Gardner Center for Engineering Applications of Radioisotopes, North Carolina State University, Raleigh, NC 27695-7909, USA
Well-designed Monte Carlo simulation programs provide the most accurate means of predicting the response of radiation measurement systems. Such simulations can replace large sets of experiments for the purpose of design optimization, calibration and analysis of sensitivity to measurement interferences. The use of detector response functions can extend Monte Carlo simulation to systems involving spectral responses. Compared to the use of general-purpose Monte Carlo codes, our work on developing specific-purpose codes has shown various advantages : greater "user-friendliness", use of an optimum set of variance reduction methods, superior code performance in most cases, and better understanding of the basic physics of the problem. The specific-purpose codes are relatively easy to develop using modules that we have developed for geometry, cross sections, detector response functions, standard as well as new and/or improved variance reduction schemes, etc. The measurement applications for which we have designed and used specific-purpose Monte Carlo simulation codes in recent years include on-line capture gamma-ray analysis of coal ; X-ray fluorescence analyzers; nuclear tools for oil-well logging by neutron moderation, pulsed neutron lifetime analysis and by gamma-ray backscatter; and industrial tomograpluc scanners. 1. Introduction Systems for radioisotope measurement applications are commonly designed by converting a measurement principle into a prototype measurement system, carrying out design calculations based on simplified models in order to approximately optimize the design, verifying the adequacy of the prototype design using an extensive set of laboratory measurements, refining the prototype system for commercial application, and calibrating the commercial device using sets of well-characterized calibration standards. These devices range in complexity of design and of the measurement principle all the way from relatively simple nuclear gauges to radiation imaging systems such as tomography scanners . As the complexity of the measurement application increases, it becomes essential to increase the sophistication of the design calculation model . One goal of this paper is to convincingly show that Monte Carlo simulation models present the best option for design calculations and that they are often sufficiently accurate and versatile for sensitivity evaluations and even for the measurement system calibration . With the recent revolution in computer capabilities, Monte Carlo simulation has become feasible without requiring access to a large-scale computational facility . It is possible to develop Monte Carlo programs which can run on relatively small computers and can simulate in reasonable time very complex geometries using detailed treatment of radiation-interaction physics and
microscopic data. What is required, however, is to tailor the program for the specific type of application . General-purpose codes such as MORSE [11 and MCNP [2] cover too wide a range of applications to be efficient and easy to use for radioisotope measurement applications. Another goal of this publication is to encourage wider use of Monte Carlo simulation within the radiation measurement community through a discussion of examples of successful specific-purpose codes which the authors and their group have developed over the last several years. Most of our Monte Carlo studies have been on systems that use neutrons, gamma-rays, or their combinations. Charged-particle Monte Carlo methods and examples will not be discussed here . In 1983, the authors published a paper [3] reviewing the development of specific-purpose Monte Carlo codes for radiation gauges and analyzers. It contrasted our experiences with general-purpose neutron-photon codes against a few of our specific-purpose codes (neutron moisture gauge, gamma-ray density gauge, X-ray fluorescence analyzers and a neutron-capture gamma-ray analyzer) and concluded that the latter are cost-effective and accurate enough to provide a nonexperimental method for calibration . Since that time, our applications have been extended to include nuclear geophyics problems, computed tomography, and the inverse problem (determination of elemental amounts from pulse-height spectra) of nuclear analyzers . The scheme for developing Monte Carlo programs has been generalized through a modular approach, and several improved methods of
O16ß-9002/90/$03 .50 © 1990 - Elsevier Science Publishers B V. (North-Holland)
I. DETECTORS
2
K Verghese, R. P. Gardner / Designing specific-purpose Monte Carloprograms
variance reduction have been formulated and implemented for deep-penetration neutron and photon problems . Specific-purpose Monte Carlo programs for radiation measurement applications are optimally formulated in a modular manner using subroutines which are common to many of the systems. Basically, these modules cover the following aspects of Monte Carlo: Generalized geometry and tracking of particle paths through the various zones, particle-collision physics, microscopic data pertinent to the radiation interactions, response functions of detectors, and variance reduction methods. Our approach to the development of these modules is briefly discussed. It is hoped that this discussion will help to remove the mystery of Monte Carlo techniques which has been perpetuated through the use of large general-purpose codes. Because of space constraints, it will be assumed that the reader is familiar with the fundamental methods of Monte Carlo for radiation transport simulation as outlined m refs . [4] and [5].
2. Generalized approach to specific-purpose Monte Carlo programs A Monte Carlo program for radiation measurement applications basically must have the following features : (1) it must be flexible enough for the user to easily make changes in the geometry of the measurement system and to track paths of particles through the various zones ; (2) the physics of all of the relevant radiation interactions must be included and sampling schemes for the interaction probability density functions must be sampled efficiently ; the microscopic physics data, such as cross sections, secondary-particle yields, etc., must be stored and accessed by the program in the most efficient manner ; (4) the response of the radiation detector (such as efficiencies, pulse-height spectral response, etc.) must be included and used in the most efficient manner ; and the variance reduction methods must be optimum for the type of measurement system, and it must be easy for the user to specify the parameters, if any, of the variance reduction schemes . It is in the last two items that specific-purpose codes differ significantly from the general-purpose codes. Spectral responses of detectors are not easily incorporated into general-purpose codes. Also, the user of a general-purpose code is usually offered a menu of variance reduction methods to choose from with the recommended method being the concept of splitting and
Russian Roulette [5] based on user-supplied importance regions and associated importance values. This method is generally not the optimum for radiation measurement problems and requires considerable experience to set up . In contrast, specific-purpose codes can be hard-wired with the optimum set of variance reduction methods for the type of problem. We will now proceed to discuss our approach to each of the above features. 2.1 . General geometric modelling
Radiation measurement system geometry can get quite complex m many applications and it is essential to have a general geometry modelling and particle-tracking program. T.H . Prettyman of our group has developed such a program called HERMETOR [6] which is being used to some of our Monte Carlo codes. HERMETOR allows the user to describe heterogeneous systems made up of homogeneous zones whose boundaries can be described by second-order surfaces . It uses a combination of boundary representation to construct zones and constructive solid geometry logic to combine the zone surfaces to form the desired geometric system . The input file for HERMETOR is simple and a 2D plotting package is used to verify the accuracy of the input. HERMETOR has proven invaluable in modelling complex systems such as nuclear well-logging tools [6] and tomographic scanners [6] which contain multiple zones of shielding, collimators and detectors. 2.2 . Neutron and photon interaction physics
The probability distributions and the appropriate methods for sampling these distributions for neutronphoton Monte Carlo simulations are well-documented in the literature (for example, ref. [5]) . In most radiation measurement applications, the source can be assumed to be a point emitting isotropically and the source energy spectra are well known either as empirical for252Cf mulae (for example, the fission spectrum) or as pointwtse tabulations as in the case of most (a, n) neutron sources . The major interactions for neutrons are elastic- and inelastic scattering and neutron-capture reactions. Secondary neutron emission reactions such as fission are rarely encountered. The angular distribution of elastically scattered fast neutrons can be treated in many cases with sufficient accuracy using the isotropic scattering assumption in the center-of-mass system . Thermal-neutron scattering is often treated adequately by the free-gas model . However, one word of caution in this regard is that the degree of sophistication required in the physics treatment is problem-dependent and if model results show differences from good experimental
3
K. Verghese, R P. Gardner / Designing specific-purpose Monte Carlo programs
benchmarks, it is prudent to consider adding a more detailed treatment of the collision physics . The physics of gamma-ray interactions are more straightforward than that of neutron interactions and are adequately treated in textbooks [7] . The major angular distribution functions that are encountered are Rayleigh scattering and Compton scattering with the classical angular distributions corrected by form factors at very low energies. The most common secondary particles emitted are electrons which, because of their relatively short range, need not be tracked in most cases .
2.3. Interaction-physics data The major source for neutron-physics data in the USA is the Evaluated Nuclear Data Files/Part B (ENDF/B) compilation at Brookhaven National Laboratory, and photon-physics data are available from the National Institute of Standards and Technology . These data banks are invaluable sources of pointwise data . Generally, the data sets are very large, and particular attention must be paid to reducing the data storage and the accessing times since the data sets might have to be accessed millions of times during the course of a Monte Carlo run . One way of compressing the data sets is to develop empirical functions that fit the pointwise data. For photon cross sections, this is feasible to do since the data generally behave smoothly . In our work, we use cubic-spline-fitted Rayleigh, Compton, photoelectric and pair-production cross sections. The fitting algorithms for continuous and piecewise continuous functions are discussed by Yacout et al . [8] . Also, for lowenergy photons, we have developed [9] and used two-dimensional cubic splines to fit cumulative probability density functions of coherent and incoherent scattering . The fits are accurate to within the error in the tabulated pointwise data and the speed-up in data access is of the order of 80 compared with binary interpolation of pointwise data . Neutron cross sections are too bumpy for fitting and are treated as pointwise binary files . Some access time can be saved by preprocessing the cross sections such that linear-linear interpolation can be performed between the data points (rather than log-linear or log-log interpolation) . Neutron cross-section files are generally the limiting aspect of computer memory requirements in applications that involve a large number of elements . The other large neutron data sets are capture gamma-ray energies and yields, nuclear-energy levels and secondary gamma-ray yields for inelastic scattering, anisotropicscattering angular distributions as a function of neutron energy, and thermal-neutron scattering data for the S(a,,8) treatment.
2.4. Detector response For all radiation measurement applications, the final tally in a Monte Carlo calculation is the score in the detector(s) . This score in many cases may be multidimensional . If only the counting rates are required, the mean score of particles interacting in the detector is multiplied by the source emission rate. However, if the counting rate as a function of energy deposition in the detector is required (as in energy spectroscopy), the entering particle must be followed through its interactions inside and outside the detector because of the potential to leave the detector and at a later time return to the detector . When applicable, a more efficient scheme is to devise and store the spectral response function of the detector as a probability-density function and to convolve the weight distribution of the particles interacting in the detector as a function of their incident energy with the response function . The spectral response function R(E, E') dE is the probability that radiation entering the detector with energy E' will deposit energy between E and E + d E in the detector . For detectors with low aspect ratio or those in applications where the radiation enters the detector only from one side, the response function can be deduced as a semi-empirical analytical function from experimental spectra for known source energies . Examples are provided in refs . [10] and [11] for a Ge detector and in refs . [12-14] for a Si(Li) detector . Jin et al . [15] provide an example of the use of the response function in conjunction with a Monte Carlo program to generate the pulseheight spectrum from a coal elemental analyzer and Yacout et al . [16] used the response function of a Si(Li) detector to calculate the spectral response of an X-ray fluorescence analyzer . In both cases, the simulated spectra matched up very well with spectra from experimental standards . An improved and more general Si(Li) response function for X-rays covering the range 5-60 keV is discussed by He et al . [17] .
2.5 . Variance reduction methods Many of the variance reduction methods that are useful in specific-purpose codes are discussed by Carter and Cashwell [5] . They generally fall within the category of biased sampling methods where instead of sampling from the actual probability-density function one samples from a fictitious density function which is known to enhance the chance of success (detection) for that particle . The bias introduced by the fictitious function is corrected by using a weight factor to alter the weight of the particle . In many cases the biasing scheme is physically obvious . The tricky problems are those of deep penetration where the particles have to get through a large number of mean free paths in order to reach the t . DETECTORS
K. Verghese, R. P. Gardner / Designing specific-purpose Monte Carlo programs
4
detector and narrowly collimated
systems where the
particles can be detected only by travelling in the "right" direction with the proper energy and weight . Unfor-
tunately many of the radiation measurement application problems fall in these categories. Examples are
nuclear well-logging systems [18] (such as the neutron porosity tool and the lithodensity tool) and photon tomography systems. The problem is that in order to
enhance the chance of success one must preferentially
sample those trajectories in the phase space that are
splitting weight . The method is a generalization of the DXTRAN scheme found in the MCNP code [2] manual and the details of implementation for photons are dis-
cussed by Prettyman et al . [23] .
The last of the powerful variance reduction schemes
that must be discussed is the concept of correlated
sampling . Correlated sampling is a scheme for studying
the sensitivity of the Monte Carlo result to small per-
turbations in the physical parameters of the problem. For example, it may be used to determine the change in
that one can a priori identify such trajectories . That is
the system response caused by small changes in the particle cross sections, density, elemental composition,
calculations or simulations. Assuming that such path-
accomplished without rerunning the Monte Carlo code
techniques to reduce the variance in deep-penetration
correlated sampling is not new, it has been implemented differently by different investigators . An efficient correlated sampling scheme that the authors have imple-
likely to contribute to the detector score and it assumes often not a simple task and must be deduced from prior
ways can be identified, we have successfully used two problems : one is path-length stretching or exponential transform [5] and the other is direction biasing [19] . In
etc., in specified regions of the phase space. This can be
for all of the perturbed samples. While the concept of
stretched using a fictitious cross section if the particle is
mented [24] consists of doing the Monte Carlo simulation in the usual way for the reference sample . Simulta-
shortened if the particle is travelling in the opposite
the perturbed samples are also calculated assuming that
path-length
travelling in
stretching,
the
distance
to
collision
is
the right direction and the distance is
direction. The net effect is that more particles are able
to get to the target (detector) . In direction biasing, the particle direction is biased towards the target following
neously during this simulation, the responses for all of
the same particle path is valid in both the reference sample and in the perturbed samples. With this assump-
a scattering interaction . A combination of the two
tion which is valid for small perturbations the only additional calculation needed for the perturbed sample
deep-penetration problems . Implementation of pathlength stretching is straightforward and requires one
weight for the reference sample is to be calculated . Thus the amount of extra calculation for a large set of per-
methods can yield significant reduction in variance in
input parameter while the implementation of direction biasing for the common neutron and photon scattering
response is the modified weight at each event where the
turbed samples (say, 50 samples) only adds a few percent to the total computer time but the results for the
kernels required considerable development work, as dis-
perturbed samples are highly correlated to the result for
Another variance reduction method which we have used in several of our simulations is an estimation
lated sampling provides an excellent way to calibrate
cussed by Gardner et al [20] .
method known as statistical estimation [21] . It is a way of estimating the score to the detector if a particle were to scatter from any point of interest towards the detector, reach the detector without further collisions along
the reference sample.
For nuclear gauge and analyzer applications, corre-
the system without having to resort to a large set of
well-characterized experimental
standard
samples. A
very small number of experiments might be necessary in
order to tie down the response of the reference sample,
the way and interact in the detector . Statistical estimation for amsotropic scattering can be a time-consuming
but the additional data points for calibration can be generated using the responses for perturbed samples.
mation .
simulation is
dure for variance reduction is possible. This scheme, which is called a combined splitting-biasing game, pre-
are discussed by Yuan et al . [26] .
calculation and Gardner et al . [21] and Mickael et al. [22] have discussed efficient ways to do statistical estiFor measurement systems with narrowly collimated geometry as in tomographic systems, a different proce-
The use of correlated sampling for density- and/or elemental-composition changes in nuclear-logging tool
discussed by Gardner et al . [25] . The
details of implementation for a nuclear coal analyzer to calculate the response as a function of coal composition Generally, the effect of changes m geometric dimen-
determines the appropriate types of collision paths that could potentially make a contribution to the detector
sions of the measurement system cannot be investigated
carrying a weight equal to the probability of achieving the desired path for that type . The nice feature of this
3. Specific-purpose Monte Carlo applications
score, and the original particle is split into multiple particles (one for each type) with each split particle
scheme is that if the split particle is sampled from an arbitrary distribution, one does not have to know the
by this method .
Three areas of applications of specific-purpose Monte Carlo simulation to nucleonic measurement systems are
K. Verghese, R . P. Gardner / Designing specific-purpose Monte Carlo programs
discussed. The versatility and accuracy of Monte Carlo simulations rather than the mechanics of implementation will be the focus of the discussion since the details of implementation can be found in the papers and student dissertations that are referenced . 3.1 . Nuclear analyzers for elemental composition We have simulated the response of two types of elemental analyzers: a neutron-capture prompt gammaray analyzer (NCPGA) for on-line composition measurement and a radioisotope source-excited X-ray fluorescence analyzer (XRF). 3.1 .1 . Neutron-capture prompt gamma-ray analyzer A NCPGA is used commercially for coal analysis primarily for the purpose of ash and sulfur measurements. The system that we have studied is similar in hardware to the CONAC system [27] developed by SAIC under the sponsorship of the Electric Power Re252Cf search Institute (EPRI) . It uses a source of fission neutrons and a Ge detector in transmission geometry . The neutrons are moderated and allowed to enter a 91 cm x 91 cm x 25 cm thick box of coal . The capture gamma-ray spectrum resulting from the absorption of
5
thermal neutrons in the coal elements is analyzed to determine the elemental concentrations . Fig. 1 schematically shows the simulated system . The measurement method is complicated by the fact that the response for any particular element is, in principle, nonlinearly dependent on the amounts of all of the other elements . As a result, accurate calibration of such a device as CONAC requires a very large set of accurately characterized, large-volume experimental standards in order to determine the interelement effects. Monte Carlo simulation combined with correlated sampling can be a powerful tool to determine the interelement effects and thereby simplify the calibration procedure. The details of the simulation procedure are given by Jin et al. [28] . Fig. 2 shows a recent comparison of the Monte Carlo simulated spectrum to the experimental spectrum for the same coal composition . Considering that the Ge detector response function that was used in the simulation is only valid between 1 .5 and 6 MeV, the uncertainties in the coal composition, and that the counting rates range over six decades, the agreement is quite good and adequate for calibration applications . It is important to note that there is no normalization factor used to normalize the two spectra.
Fig. 1. Schematic configuration of a laboratory simulation system for the elemental analysis of coal : (1) polypropylene, (2) paraffin work, (3) paraffin and lithium carbonate, (4) lead, (5) polyethylene, (6) coal, (7) paraffin wax and lithium carbonate, (8) Ge detector, ( " ) neutron source . All dimensions are in cm . 1. DETECTORS
K. herghese, R. P. Gardner / Designing specific-purpose Monte Carlo program ,
10
6
10
5
10
Uô
3
10
2
10
6 E
7
I
I
I
8
9
10
(MeV)
Fig. 2. Comparison of the measured gamma-ray pulse-height spectrum of a coal sample to the spectrum generated by Monte Carlo simulation for the same coal composition .
3.1 .2 . X-Ray fluorescence analyzer XRF analyzer systems are commonly used in analytical laboratories to do elemental analysis with a sensitivity in the range of parts per million for all but the very
light elements . Although
the spectral-response
simulation techniques are equally well applicable to
radioisotope-source-excited systems as well as to X-ray tube-excited analyzers, we have chosen to demonstrate
the versatility of Monte Carlo simulation on a sourceexcited system that is available for experimental validation in our laboratory .
Fig. 3 shows an XRF spectrum simulated by our
specific-purpose Monte Carlo code for austenitic stainless steel containing six elements and excited by a 109Cd
source and compared with the experimental spectrum .
Once again there is no normalization performed between the simulated and experimental spectra. The agreement is excellent in the fluorescent region . The
disagreement in the source backscatter region is because the Si(Li) detector response function that was used was valid only at lower energies than the backscattered source energies. The new response function by He et al . [17] which was referred to in section 2.3 is valid up to 60 keV and should perform better in the comparison.
So far, what we have demonstrated is that Monte Carlo simulation using fundamental data can be used to predict the spectral response of analyzer systems where the response is nonlinearly dependent on the elemental concentrations . This is referred to as the "forward"
problem m contrast to the "inverse" problem which is the problem of obtaining elemental concentrations from
measured spectral response of the system . The classical way to decompose a composite spectrum is by using the library least-squares (LLS) method, m which one re-
quires the shapes of spectra from each of the elements that are present in the composite. Unfortunately, in a
nonlinear system these shapes are, to say the least, difficult to obtain experimentally since, in principle, the
intensity and shape of the spectrum from each element will depend on the amounts of the other materials that are present in the sample matrix . However, they are easily obtained by Monte Carlo simulation by simply tallying the contribution to the spectral response from
each element separately . Thus a combination of the Monte Carlo and LLS methods may yield a way of solving the inverse problem. We are starting to apply I
o'
5
7
9
11l
I 13 15 17 ENERGY 1KEV1
19
21
23
25
27
Fig. 3 . Comparison of the measured and simulated pulse-height spectra of the backscattered 109Cd photons and the fluorescent X-rays from a 304 stainless steel sample containing Si, Cr, Mn, Fe, Ni and Mo (- Monte Carlo simulation, - - - experimental).
the
Monte
Carlo-Library
least-squares
(MCLLS)
method to both the coal analysis problem and to the
XRF analysis problem. The MCLLS procedure consists of generating the library spectra from Monte Carlo simulation based on an initial estimate of the sample composition, perform the LLS analysis using these
K Verghese, R.P Gardner / Designing specific-purpose Monte Carloprograms
library spectra to calculate the elemental composition of the sample and if the new composition is quite different from the assumed composition, alter the assumed composition and repeat the procedure . The preliminary results of MCLLS for the XRF problem appears promising [29] and the work on applying it to the coal analysis problem is continuing. If successful, the MCLLS method offers the potential for eliminating expensive and tedious experimental calibrations. 3 .2 . Nuclear-tool responses for oil-well logging Nuclear-measurement methods are commonly used in logging bore holes for petroleum exploration and production . Such measurements routinely provide information on porosity of the formation, salinity and percent saturation of pores, formation density, effective atomic number, etc ., and are used in conjunction with other measurements to deduce the necessary geophysical data . The tools used for these measurements consist of a radiation source (neutrons or gamma-rays) and a set of two or more detectors . The detector responses are analyzed to determine the geophysical quantity of interest as a function of depth down the bore hole . The tools are generally run along the side of the bore hole . There are three nuclear tools for which Monte Carlo codes are being developed . These are the dual-spaced neutron-porosity tool, the pulsed-neutron lifetime tool, and the lithodensity gamma-ray tool . 3.2.1 . Dual-spaced neutron-porosity tool The dual-spaced neutron-porosity tool is a 7 .6 cm diameter tool that contains a few-Ci Am-Be neutron source and two 3 He neutron detectors mounted axially at two different distances from the source . In addition, such tools contain various shield materials placed strategically to enhance the porosity sensitivity . The ratio of the counting rate on the near-spaced detector to that on the far-spaced detector correlates well with the porosity of the formation and it is nearly independent of the bore hole conditions . For Monte Carlo simulation, this system forms a deep-penetration problem since the far-spaced detector is placed generally at around 50 cm from the source and hence requires extensive use of variance reduction schemes. General-purpose codes such as MCNP [2] do this by splitting up the geometric regions of the problem into a grid and neutrons which enter a grid region are either split into multiple neutrons with reduced weight or killed by Russian Roulette based on how important that region is to the final score. The set of importance values for the grid must be provided by the user . For further variance reduction, other schemes such as weight windows are used to avoid excessive weight fluctuations for the neutrons which are detected . In contrast we have preferred the combined use of path length stretching,
7
direction biasing, and statistical estimation in our code, McDNL [30] . This makes McDNL considerably easier to use than MCNP for this problem . Correlated sampling is an integral part of this code and has proved invaluable for calculating environmental corrections in porosity logging [31] . The McDNL code has been in use in the petroleum industry now for about two years . Although the previously reported computing times [30] are comparable with that of MCNP, recent improvements make it significantly faster than MCNP . A major attractive feature of McDNL is its user-friendliness. 3.2 .2 . Pulsed-neutron lifetime tool The objective of the pulsed-neutron lifetime measurement is to determine the macroscopic absorption cross section of the formation . The tool is a 4 .3 cm diameter tube which houses a (13,T) neutron generator and two gamma-ray detectors in addition to shielding configurations . A 14-MeV neutron burst lasting for several microseconds is emitted from the generator . The neutrons moderate and are absorbed generating gamma-rays from inelastic neutron scattering and radiative capture. These gamma-rays are detected as a function of time measured from the neutron burst . This temporal spectrum is analyzed to yield the absorption cross section of the formation. Certain pulsed-neutron tools also can record the gamma-ray pulse-height spectra from which elemental analysis for the formation can be performed . For Monte Carlo simulation, this is a combined neutron/photon problem requiring detailed physics treatment because of the high neutron energy and the importance of inelastic neutron scattering. Also, the fact that responses in up to 100 time channels must be tallied with good statistical accuracy complicates the simulation somewhat . However, because of gamma-ray detection rather than thermal neutrons, path-length stretching and direction biasing have not been necessary as variance reduction methods. Statistical estimation for gamma-rays is the major variance reduction method . Also, as in the other examples, correlated sampling provides a useful addition to study environmental effects on log response . McPNL, our Monte Carlo code for the pulsed-neutron logging tool has been benchmarked against a set of laboratory test pit data and the simulated results agreed exceptionally well with the experimental results [32,33] . 3.2 .3 . Gamma-ray hthodenstty tool This tool is a gamma-ray backscatter tool for the determination of formation density and to provide a measure of the lithology . It consists of a gamma-ray source (typically 135 Cs) collimated to emit directly into a cone within the formation and contains two collimated gamma-ray detectors . The pulse-height spectrum of I. DETECTORS
8
K. Verghese, R P Gardner / Designing specific-purpose Monte Carloprograms
scattered photons is recorded and used to differentiate between different formation lithologies. The complex part of simulating this tool is the low yield of counting rate per source photon (of the order of 10 -s ) resulting from the narrow collimation. A code called McLDL is currently under development and we expect to have the first version ready by the end of the year. 3.3 . Tomographic scanners for industrial tomography
Industrial tomographic systems are becoming widespread and diverse in design from the tomographic systems used in radiology. These systems are complicated and expensive to design and build. The complicating features are stringent collimation, multiple detectors and generally complex geometry. Simulation will play a decisive role in their design development and optimization. Accurate simulated data can be used also to develop image reconstruction methods. The power of Monte Carlo for tomographic simulation has been long recognized and used [34] to some extent for transmission-tomography systems. We are investigating the design of a system that would use both transmitted and scattered photon detector signals to image the atomic number as well as the density distributions in industrial samples. Monte Carlo simulation verified by relatively simple experiments in the laboratory is the approach taken. A code called MCPT [23] has been written for this purpose and is being benchmarked using scattered photon pulse-height spectra from inhomogeneous phantoms . The code incorporates the geometric modelling package, HERMETOR, and a general combined splitting-biasing scheme to send photons to detectors without scattering and after single, double and higher-order scatters in the sample as well as in the other important zones of the scanner. The simulated responses will be used to develop reconstruction methodology from the combined transmission and scattered pulse-height spectra. 4. Discussion and conclusions We have reviewed the past seven years of our work on the development and uses of specific-purpose Monte Carlo programs for a variety of radiation-measurement apphcations . All of our programs are developed on the relatively small Vax system of computers including the MicroVax 11 and DecStation 3100, but they are written in standard Fortran and are portable to more powerful systems such as Cray systems. Accurate CPU times for the various applications and the relative standard deviations of the simulated responses are given in the individual references that are cited, but generally the run times are of the order of one to a few hours on the DecStation
which runs approximately a factor of 15-18 times faster than the MicroVax II . In the early days, each of our specific-purpose codes was written independently of each other, but over the years we have generated a large set of subroutines that are common to many applications . The methodologies used in these subroutines are reviewed in section 2 of this paper and for details the reader is again referred to the set of references that are cited in that section. For a completely new application, typically the time it takes to assemble a working code is no more than the time to obtain a general-purpose code, install it and learn to run it in the optimum manner . But the major advantages of the specific-purpose codes are their ease of use and the fact that by writing the code one learns a great deal more about that application than if one just learns to set up the input to a general-purpose code. We have chosen to illustrate the versatility of specific-purpose Monte Carlo by discussing a wide range of fairly complex problems in radiation applications . Many of the applications have necessitated the development of new or improved implementation methods for geometric modelling, data fitting, detector response functions and variance reduction. Programs for some of the applications are still under development and improved versions of the others such as McDNL and McPNL are also being developed. Acknowledgements The authors acknowledge the interest, encouragement and financial support provided by the Amoco Production Company, Mobil Research and Development Corporation, BP Research, Arco Oil & Gas Company, Exxon Production Research Company, Atlas Wireline Services, Teleco Oilfield Services Inc ., and Conoco Inc., under an Associates Program supporting Nuclear Oil Well Logging Research and BP America for support under a separate contract for Industrial Tomography research . References [I] E.A . Straker, W H. Scott, Jr, and N.R. Byrn, Radiation Shielding Information Center Report CCC-127D (Oak Ridge National Laboratory, 1972) (available from National Technical Information Center, Springfield, VA, USA) . [2] Los Alamos Monte Carlo Group, Los Alamos National Laboratory Report LA-7396-M(rev .) (Los Alamos National Laboratory, April 1981). [3] R.P . Gardner and K. Verghese, Isotopenpraxis 19 (1983) 329. [4] W.E. Selph and C.W . Garrett, m. Reactor Shielding for Nuclear Engineers, ed . N.M . Schaeffer (U .S . Atomic Energy Commission, TID-25951, 1973) p. 207.
K. Verghese, R. P. Gardner / Designing specific-purpose Monte Carlo programs [5] L.L . Carter and E.D . Cashwell, Particle-Transport Simulation with the Monte Carlo Method, ERDA Critical Review Series TID-26607 (Technical Information Center, Springfield, VA, 1975). [6] R.P. Gardner, T.X . Prettyman, M. Mickael and K. Verghese, Trans. Am . Nucl . Soc. 57 (1988) 115. [7] R.D . Evans, The Atomic Nucleus (McGraw-Hill, New
York, 1955). [8] A.M . Yacout, R.P . Gardner and K. Verghese, Nucl . Instr. and Meth. 220 (1984) 461. [9] A.M. Yacout, R.P . Gardner and K. Verghese, X-Ray Spear. 15 (1986) 256. [10] Y. An, R.P . Gardner and K. Verghese, Nucl. Instr. and Meth . A242 (1986) 416. [11] M.C . Lee, K. Verghese and R.P . Gardner, Nucl . Instr. and Meth . A262 (1987) . [12] R.P. Gardner and J.M. Doster, Nucl . Instr. and Meth. 198 (1982) 381. [13] A.M. Yacout, R.P. Gardner and K. Verghese, Nucl. Instr. and Meth . A249 (1986) 121. [14] R.P. Gardner, A.M . Yacout, J. Zhang and K. Verghese, Nucl. Instr. and Meth . A242 (1986) 399. [15] Y. Jin, R.P . Gardner and K. Verghese, Nucl. Geophys. 1 (1987) 167. [16] A.M . Yacout, R.P . Gardner and K. Verghese, Adv. X-Ray Anal. 90 (1987) 121. [17] T. He, R.P. Gardner and K. Verghese, these Proceedings (7th Symp. on X- and Gamma-Ray Sources and Applications, Ann Arbor, MI, 1990), Nucl. Instr. and Meth . A299 (1990) 354. [18] D.V . Ellis, Well Logging for Earth Scientists (Elsevier Science Publ. Co., Amsterdam, 1987) chs. 9, 12. [19] J.S . Hendricks and L.L . Carter, Nucl . Sci. Eng. 89 (1985) 118. [20] R.P . Gardner, M. Mickael and K. Verghese, Nucl . Sci. Eng. 98 (1988) 51 .
9
[21] R.P. Gardner, H.K. Choi, M. Mickael and K. Verghese, Nucl. Sci. Eng. 95 (1987) 245 [22] M. Mickael, R.P. Gardner and K. Verghese, Nucl . Sci. Eng. 99 (1988) 251 . [23] T.H . Prettyman, R.P . Gardner and K. Verghese, these Proceedings (7th Symp . on X- and Gamma-Ray Sources and Applications, Ann Arbor, MI, 1990), Nucl . Instr. and Meth . A299 (1990) 516. [24] T.C . Clark, R.P . Gardner and K. Verghese, Nucl . Instr. and Meth . 193 (1982) 365. [25] R.P . Gardner, M.W . Mickael, C.W. Towsley, C.M . Shyu and K. Verghese, IEEE Trans . Nucl . Sci., in press. [26] Y.L . Yuan, R.P. Gardner and K. Verghese, Nucl. Technol.
77 (1987) 97 . [27] G. Reynolds, H. Bozorgmanesh, E. Ehas, T. Gozani, T. Maung and V. Orphan, EPRI Report EPRI-989, vol. 1 (1979) . [28] Y. Jm, unpublished Ph .D . Dissertation, North Carolina State University (1986) (available from University Microfilms, Ann Arbor, MI). [29] K. Verghese, M. Mickael, T. He and R.P . Gardner, Adv. X-Ray Anal . 31 (1988) 461. [30] M. Mickael, R.P . Gardner and K. Verghese, Proc . SPWLA 29th Ann. Logging Symp ., Paper MM (Society of Professional Well Log Analysts, 1988). [31] M. Mickael, H.K . Choi, R.P . Gardner and K. Verghese, Nucl . Geophys. 3 (1989) 339. [32] H.K . Choi, K. Verghese and R.P . Gardner, Nucl. Geophys . 1 (1987) 71 . [33] H.K . Choi, K. Verghese and R.P . Gardner, Nucl. Geophys . 3 (1989) 167. [34] G.T . Herman and S. Rowland, Technical Report no. 130 (State University of New York, Buffalo, 1978).
I. DETECTORS