- Email: [email protected]
Algorithm for computer processing and analysis of gamma ray spectra for large-scale laboratory operations
NUCLEAR
INSTRUMENTS
AND
METHODS
I37
(t976)
I73-I78;
©
NORTH-HOLLAND
PUBLISHING
CO.
A L G O R I T H M FOR C O M P U T E R P R O C E S S I N G A N D ANALYSIS OF G A M M A RAY SPECTRA FOR LARGE-SCALE L A B O R A T O R Y O P E R A T I O N S B R A D L E Y A. R O S C O E and A. K E I T H F U R R
Nuclear Science and Engineering, Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061, U.S.A. Received 11 May 1976 For large scale neutron activation analysis laboratories, computer analysis o f data is necessary. This work discusses major features o f the analytical program which has been developed and is currently being used at Virginia Polytechnic Institute and State University. Specific attention is drawn to fitting function choice, controlled fitting techniques, and methods o f error analysis. In order to test the program, samples o f standard reference material from the National Bureau o f Standards were obtained and analyzed. The results o f this analysis are given and demonstrate that the code is effectively doing its intended task.
1. Introduction With the acceptance of neutron activation analysis (NAA) as an accurate method of elemental analysis, several laboratories, such as the one at Virginia Polytechnic Institute and State University (VPI & SU), have made their services widely available. Analysis has been carried out on automobile emissions, organic materials, soils, industrial effluents, and a multitude of other types of samples in efforts to establish the environmental impact of today's society on the ecology. Another important application is determining trace element concentrations in an attempt to find correlations between these trace elements and the sample's general characteristics, i.e. leanness of beef, susceptibility of teeth to decay, etc. As a semi-commercial service, the N A A Laboratory at VPI & SU analyzed over 3000 samples in 1975. This service was extended primarily to on-campus groups but was also available to other institutions on a timeavailable basis. For each sample analyzed, there were at least three gamma spectra associated with it such that over 9000 individual g a m m a spectra were analyzed in that year. For this reason, the need of an effective computer code to analyze g a m m a spectra can be realized. There have been several programs for g a m m a spectra analysis developed over the years1-3); however, when development of our current program started at VPI & SU in 1972, only relatively unsophisticated ones were readily available. The program currently being used at VPI & SU has developed into a very powerful and inexpensive analytical tool for spectra analysis. Typically, the cost of a single element identification and quantitative determination is on the order of 2-3 cents for computer time. It runs on an IBM 370,
mod 158 computer and requires approximately 120k for compilation and execution. It is written in modular form so that, if necessary, it could be adapted to run on a much smaller system. 2. The fitting function The selection of a proper fitting function is of primary importance for data evaluation. In the selection of this function, one must consider all the limiting factors in the entire counting system from detector to memory. The first consideration is to the energy range to be observed and the number of data channels available with which to store data. In fitting a function to data using a least-squares method, there must be at least one more data point than parameter in the function. A fine point of least-squares fitting techniques that should be noted is that in order to obtain a good quality fit, the number of points in the actual peak, not including the background portion of the data, should be greater than the number of parameters in the fitting equation. For the counting system currently used at VPI & SU, an energy range from 60 keV to 3200 keV is observed using 4096 channels. This configuration yields about 0.8 keV per channel for acquired data. With peak width (full width at half maximum) of about 2.0 keV and a total width at 0.1 maximum of about 4 keV, an average peak will only have about 7 or 8 data points clearly distinguishable from background. Also to be considered in choosing the fitting function is the actual shape of the photopeak. Ideally, a photopeak should take the form of a Gaussian. However, due to imperfections in detectors and amplifiers, the Gaussian is always distorted. The most prominent cause of the distortion is produced by the amplifier
174
B R A D L E Y A. R O S C O E A N D A. K E I T H F U R R
giving the Gaussian photopeak an exponential tail on the low energy end due to inadequate base-line restoration at elevated counting rates. In order to reduce this effect, it is necessary to use an amplifier with superior baseline restoration such as the O R T E C model 472 used at VPI & SU. Possible fitting functions suggested by other researchers have included Gaussians plus: 1) an exponential tail on the low end4), 2) exponential tails on both endsS), or 3) exponential tails on both ends distorted by a power series6). In all the above functions, a quadratic was included in the function to account for background. It should be noted that in the formulation of these functional shapes, the experimenters were using data with at least 20 points and frequently more in the photopeak. In the analytical program discussed here, the exponential tail is ignored since an amplifier with good baseline restoration is used, moderate count rates are employed, and 8 data points will not show the exponential tail structure prominently. Thus the fitting function used is as simple as possible, i.e.,
established. The fitting function used at VP1 & SU is obviously non-linear, thus a method of non-linear leastsquares fitting is used. For this purpose, the "algorithm of Marquardt"v), which uses an analytical type gradient search, is used with extremely good results (see figs. 1 and 2). A typical peak will require only 3 iterations to obtain a satisfactory fit but up to 8 are allowed in attempting to fit the data. Since more than one solution is possible with a non-linear least-squares fit, it is necessary to control the parameters during the fit so that the best solution is obtained; thus the lowest chi-square is obtained instead of a minimum in the chi-square hypersurfaceT). This is achieved by making good first guesses for parameters and by preventing them from wandering from established limits. Initial guesses for parameters must be picked as closely as possible to the nominally correct solution. The peak location and width are the most important parameters that need to have good first guesses. In the 8
8
Y(1) = A(1) exp -
\
A(5)
] _]
+ A(4) + A(2) X(I) + A(6) X(I) 2 ,
~
BackgrOund+ Bromine Peak
|
(1)
where Y(I) is the number of counts in the X(1) channel; A(2), A(4), and A(6) represent the quadratic background; and A(1), A(3), and A(5) are the height, location, and width of the peak. In the event of a doublet being present (two peaks so close together that they overlap), an additional term is added to the functional equation to compensate for the interference,
Z
.
o'
.
xooLoo
.
,
I
.
.
Loos.~o
.
:u:..
i rol.o,°o P: ,k r
Lo~oo
I
xozx oo
I
Lo~7. ~
r
z~z. oo
NORPIRLIZED CHRNNF--L LOCRT[ON
~.m
J
l . o ~ oo
Fig. 1. Least-squares fit to a bromine peak at 698 keV using 6 variables [eq. (1)] and 29 data points. 8
A(7) exp
- 2
A(5)
'
(2) 8.
where A(7) and A(8) are the height and differential location of the interference peak. In order to minimize the number of free parameters, it is assumed that the width of the two overlapping peaks will be identical. It has also been found that having A(3) and A(8) coupled together for the location of the second peak will aid in the rate of convergence of the routine.
8 Background/
o~ 3. Method of curve fitting
In photopeak analysis, it is necessary to fit the acquired data to the desired fitting function so that the location and area of the peak can be accurately
~aoa
L'o~_00 ~ 0
NORHRLIZE0
~00
~co
~'a~o0
CHFiNNIZL L O C R T I O N
LIOLZ.EIO
Fig. 2. Least-squares fit to an iron peak at 1099.2 keV and an unknown interference using 8 variables [eqs. (1) and (2)] and l0 data points.
ALGORITHM
175
FOR C O M P U T E R P R O C E S S I N G
analysis program at VPI & SU, the program first calculates a channel vs energy and width vs energy calibration for each spectrum so that a good first guess can be made for the location and width of the peak. For the channel vs energy calibration, two known peaks are located such that a linear correlation may be made. One of the peaks is from 241Am at 59.543 keV which is present in every gamma spectrum from a source mounted near the detector. Due to the numerous possible interferences at very low energies from many of the rare earths, it has been decided not to attempt to use data below approximately 70 keV. Hence, the presence of the 2 4 1 A m peak does not adversely affect any data of interest. For a second calibration point, one of the following peaks is located: :4Na at 2754.1 keV, :SA1 at 1778.8 keV, or ~°K at 1460.75 keV, the choice being dependent upon the decay interval. In almost all cases, one of the above peaks will be present, such that the linear calibration can be made. Minor deviations from a linear fit can be accommodated by slight modifications of the library photopeak energies if necessary. The width vs energy calibration is made from the determined width of either the Na, AI, or K peaks. It has been experimentally verified that a linear width vs energy curve will vary only in the intercept of the linear equation and not in the slope term (see fig. 3). Thus, since the slope of this calibration is known, a new intercept can easily be found for the calibration from the one data point. Between iterations, the width, height, and location parameters are checked to insure that they have not wandered from pre-established limits. If one has deviated too much, it is reset to a more reasonable
value and the interation continued such that the desired solution may be obtained. When looking for a specific photopeak, the program will fit the data to a gaussian plus a quadratic background. If in the final fit, the width is too wide (see fig. 4), the data will be re-fitted using the equation for a doublet, thus accounting for interference peaks. The method of controlled fit to the non-linear equation has been proven effective in the analysis of complex spectra. It is quick, accurate, and needed in the computer analysis of data for large scale analysis operations. 4. Error estimates
As in any experimental procedure, it is desirable to have a method of error analysis so that a confidence level in results can be realized. In the analysis of gamma spectra, two different types of error are considered, statistical error and fitting error, both of which can be calculated. By examining the value of these two quantities, the interpreter of the data can achieve a level of confidence in the final data. The statistical error is obtained from basic counting theory which states that the standard deviation of the number of counts is equal to the square root of the number of counts. Since peak areas are calculated over several channels, a slight variation is used such that the statistical error, equivalent to one standard deviation, is calculated from S.E. = ~/[A(peak) + A(bkg)] A (peak)
x
loo Vo,
(3)
where A(peak) and A(bkg) are the areas of the
D ipl"~
Oe~dTm i e~42%~ UpperLimit
(I wtq
0-_ i -
~
^
•
~ . _ ~DeadTi. . . .
ENERGY (XIDw Jl
Fig. 3. Peak width times. Data taken represented by ' © ' peak width is the channels.
vs energy (keV) for different system dead at 42% dead time and 2% dead time are and ' × ', respectively. On the ordinate, the Gaussian width o f the peak measured in
5.oo
"~.m
n~rn
ENE'OO FIGY C~ID,, @
Fig. 4. G r a p h of peak width vs energy (keV). The lower and upper limits are used as trips by the computer p r o g r a m to aid in the peak fitting process. If during the fit the new calculated width of the peak falls outside the limits, it is reset to the calibration line and fitting continued with this new value.
176
BRADLEY
A. R O S C O E
AND
A. K E I T H
FURR
TABLE 1A Orchard leaves.
Element
Ca K Mg Fe Mn
Na Zn As Cu Rb Hg Br CI Cr Co S
Gamma energy (keV)
NBS a cone. (ppm)
3084.4 1524.7 1014.4 1099.2 1291.5 846.6 1811.2 2112.6 1368.5 2754.1 438.9 1115.5 559.1 657.1 511.0 1039.0 1077.0 77.6 698.3 776.5 1642.4 2167.5 320.1 1173.2 1332.5 3102.4
20900+300 14700±300 6200±200 300 ± 20
Cone. (ppm)
91 ±
4
82±
6
25±
3
14±
2
12±
1
12± 1 0.155±0.015 (10) (700) (2.3) (0.2) (2300)
21600 15390 6174 279 274 93.4 91.2 95.2 98.6 91.2 17.6(N) 24.2 14.9 13.8 12.9 9.1 12.0 0.188 8.6 8.6 695 732 2.17 O. 18 0.21 2744 (N)
Sample 1 Count error (%)
Fit error (%)
1.3 0.7 1.6 7.6 10.4 0.5 1.5 2.2 0.8 1.1 52.4 5.8 0.6 0.1 3.8 22.6 7.8 9.7 1.5 0.8 2.4 2.0 I 1.1 18.4 16.7 58.8
67.4 2.3 8.7 184.4 40.9 0.9 5.0 0.3 1.0 14.5 3.3 169.6 0.I 2.6 5.1 34.6 130.4 0.3 2.0 1.9 6.8 65.6 18.2 65.1 217.4 69.0
Cone. (ppm) 21200 14540 6436 294 369 92.4 89.6 91.3 65.1 68.2 23.8 24.2 14.4 13.4 16.8 (W) 7.6 (N) 12.0 0.175 (N) 8.2 8.4 709 688 2.9 0.22 0.26 4818 (N)
Sample 2 Count error (%)
Fit error (%)
1.4 0.7 1.7 6.7 7.9 0.5 1.6 2.4 0.9 1.2 47.6 5.7 0.5 2.2 3.2 25.5 7.3 8.5 1.5 0.7 2.5 2.2 8.6 14.4 14.4 34.7
26.5 ! .3 10.1 44.3 37.7 0.5 21.0 32.2 0.9 10.2 128.6 75.3 0.1 3.3 0.8 3.8 221.9 0.2 0.9 2.1 5.8 31.4 7.0 0.9 172.6 192.9
a "values in parentheses are uncertl.eU. Gauss±an peak and quadratic background over the i n t e r v a l o f t h e p e a k ' s full w i d t h a t h a l f m a x i m u m . The method of determining the error due to the least-squares fitting process has only been developed for l i n e a r cases. F o r n o n - l i n e a r t e c h n i q u e s , it is a s s u m e d t h a t t h e e r r o r c a n be c a l c u l a t e d in a s i m i l a r m a n n e r . I n l i n e a r l e a s t - s q u a r e s fitting, a s q u a r e s y m m e t r i c m a t r i x is g e n e r a t e d t h a t c o n t a i n s t h e w e i g h t e d d e r i v a t i v e s o f all t h e p a r a m e t e r s in t h e function. Taking the inverse of this matrix forms an e r r o r m a t r i x w h i c h c o n t a i n s t h e v a r i a n c e s a n d cov a r i a n c e s o f all t h e f u n c t i o n a l p a r a m e t e r s . T h e determinant of this matrix represents the overall v a r i a n c e o f t h e fit. By d i r e c t a n a l o g y , t h e v a r i a n c e o f t h e fit f o r t h e n o n - l i n e a r case is t h e s a m e d e t e r m i n a n t . T h u s , since t h e s t a n d a r d d e v i a t i o n is t h e s q u a r e r o o t o f t h e v a r i a n c e , t h e fit e r r o r is c a l c u l a t e d f r o m : F.E. =
fNSD " "
] '/N
IA(od
x 100 °/ /o,
(4)
where S.D. = s t a n d a r d d e v i a t i o n o f t h e fit, N = number of variable parameters, and A (i) = p a r a m e t e r n u m b e r i. It s h o u l d b e n o t e d t h a t t h i s m e t h o d o f c a l c u l a t i n g e r r o r is a q u a l i t a t i v e m e a s u r e o f t h e g o o d n e s s o f fit a n d is not absolute. T h e s t a t i s t i c a l e r r o r a n d t h e fitting e r r o r p r o v i d e useful i n f o r m a t i o n t o t h e e x p e r i m e n t e r a n d a l l o w h i m t o e s t a b l i s h l i m i t s b e l o w w h i c h h e will c o n s i d e r d a t a t o b e a c c e p t a b l e . H o w e v e r , in s o m e i n s t a n c e s , r e s e a r c h w o r k e r s a r e i n t e r e s t e d in w h e t h e r t h e r e is a n y p o s s i b i l i t y a t all o f a n e l e m e n t b e i n g p r e s e n t in t h e s a m p l e . T h e r e fore by providing the raw values of the two error e s t i m a t e s , t h e c h o i c e o f w h e t h e r t o use d a t a is t r a n s fered to the researcher and has not been predetermined by the program.
5. Calculation of elemental concentration T o m i n i m i z e e r r o r s in c a l c u l a t i n g e l e m e n t a l c o n c e n t r a t i o n , s t a n d a r d s h a v e b e e n a n a l y z e d o n all
ALGORITHM
177
FOR C O M P U T E R P R O C E S S I N G
TABLE 1B Orchard leaves.
Element
G a m m a energy (keV)
Dy La
Mo Sb
Sm Ba Cr Sc Cs W Au
Conc. (ppm)
94.6 108.2 328.7 487.0 815.8 1596.2 141.0 564.0 645.8 1691.0 69.6 103.1 216.0 496.3 320.1 889.0 1120.5 795.7 479.4 685.7 411.7
0.14 0.14 0.84 0.91 0.92 0.91 0.19 (N) 2.52 4.6 3.5 0.64 0.48 26. 38. 1.25 0.038 0.035 0.087 0.139 0.097 0.00358
Sample 1 Count error (%)
Fit error (%)
31.9 6.0 2.8 1.8 3.8 2.0 35.9 0.6 12.5 5.0 3.5 0.8 21.0 14.9 14.1 5.5 6.2 32.9 42.8 38.9 12.1
Conc. (ppm)
3.1 0.9 0.7 3.7 2.9 4.4 0.0 0.0 36.1 205.0 0.3 0.1 6.3 24.1 5.7 35.1 223.2 22.4 3.7 19.1 1.8
Sample 2 Count error
0.179 0.156 0.88 0.92 0.83 0.91 0.202 (N) 3.2 4.0 (N) 3.9 0.60 0.44 36. 39. 1.04 0.0038 0.0043 0.067 0.103 0.16 0.00398
(%)
Fit error
(%)
22.1 6.4 2.5 1.7 3.6 1.8 30.8 0.5 12.7 4.3 3.4 0.8 15.6 13.2 14.4 5.1 5.0 43.4 46.5 28.9 9.7
1.2 2.5 0.7 3.0 3.5 4.3 0.4 0.1 I.I 189.4 0.2 0.1 6.2 6.3 2.8 30.1 47.9 22.5 1.4 4.6 2.8
Sample 2 Count error
Fit error
(%)
(%)
1.0 0.2 0.2 8.5 10.2 0.6 1.4 27.1 2.6 5.8 1.3 5.4 8.4 16.5 34.3 0.8 0.8 11.8 10.7 7.1 5.4
4.8 0.2 0.3 211.9 129.7 0.1 7.5 1.5 22.5 38.8 1.3 14.5 8.5 23.4 37.2 8.5 1.8 319.7 187.6 0.9 0.2
TABLE 2 Bovine liver.
Element
K Na Fe Cu Zn Rb Mn
Se Ca CI Co Mg Mo
Gamma energy (keV)
1524.7 1368.5 2754.1 1099.2 1291.5 511.0 1039.0 438.9 1115.5 1077.0 846.6 1811.2 2112.6 161.6 3085.0 1642.4 2167.5 1173.2 1332.5 1014.4 141.0
NBS a conc. (ppm)
9700 2430
4-600 4- 136
270
4- 20
193
4- 10
130
4- 10
18.3410.3 4-
I I
1.1 + 0.1 (123) (2600) (0.18) (605) (3.2)
a Values in parentheses are uncertified.
Conc. (ppm)
Sample 1 Count error (%)
10320 2003 2027 234 315 275 (W) 192 89 103 19.3 11.2 11.0 8.9 1.8 74 (N) 2891 2806 0.39 0.38 795 2.6
1.1 0.8 0.3 9.8 10.4 0.6 1.5 24.7 3.0 6.5 1.4 6.2 11.8 26.7 66.2 0.9 0.9 11.9 11.8 6.2 6.0
Fit error (%)
3.8 0.2 0.5 462.2 466.5 0.8 11.0 2.6 250.4 118.8 1.4 6.3 0.5 4.4 34.0 5.8 10.6 118.1 225.2 5.9 0.1
Conc. (ppm)
10080 1920 1955 224 246 264 (W) 181 67 107 18.2 10.3 10.6 11.7 1.6 125 2640 2666 0.30 0.35 573 (N) 2.4
178
BRADLEY A. ROSCOE AND A. KEITH FURR
elements that are sought by the VPI & SU laboratory. In this manner, inaccuracies in branching ratios, cross section, neutron spectrum, and natural abundances have no effect on final data evaluation. Elemental concentrations are then calculated by standard methods from the peak area and activation times. In order to minimize computer time and eliminate confusing answers appearing in the output, numerous flags have been incorporated in the program so that the program only looks for individual elements under conditions predetermined by the activation, decay, and counting times. Although certain rough rules have been incorporated in setting these flags, empirical factors, such as the presence of interferences, have proven more important and hence the flags have been determined laboriously, element by element.
6. Program testing In order to test the accuracy of VPI & S U ' s analysis program, samples of standard reference material, such as bovine liver and orchard leaves, were obtained from the National Bureau o f Standards. In the data presented here, two samples of each material were prepared and analyzed using the same procedures that any other samples undergoing analysis would receive. In this manner, the overall accuracy of the laboratory was tested. The results of the analysis of the reference materials are given in tables 1A, 1B, and 2. Tables 1A and 2 give comparison of certified and uncertified NBS concentrations while table 1B gives elemental concentrations contained in the sample that were found at VPI & SU and not reported by the NBS. In the tables, an (N) or (W) represents a narrow or wide peak as determined by the c o m p u t e r ' s width vs energy calibration. As should be expected, the concentration for a narrow peak tends to be low while the value for a wide peak tends to be high. It should also be noted
again that the fit error is only a relative measure and not absolute. As can be seen from the data, the results obtained correlate very well with NBS certified and uncertified values except for the cobalt in bovine liver and sulfur in orchard leaves. Since N A A is not very sensitive to sulfur, the result is not surprising; however, for sulfur, an order o f magnitude answer is reasonable. The cobalt in the bovine liver disagrees by about a factor of two. We believe that the uncertified NBS value is off by this factor due to our extreme sensitivity to cobalt and the fact that the cobalt data for the orchard leaves compared reasonably well. Another factor providing confidence in our cobalt value is that cobalt is our basic standard for determining the thermal neutron flux in the VPI & SU reactor. F r o m the results obtained from the NBS reference material, it can be seen that the analysis program employed at VPI & SU is doing its intended task very effectively. It has vastly reduced data analysis time, decreased errors and increased the a m o u n t of data that can be processed at the laboratory.
References 1) T. Meyers and M. Denies, Longterm and Peakscan: neutron activation analysis computer programs, Michigan: Museum of Anthropology (1972). 2) Radionuclide program: Peakfind and Isoquan, Ortec Incorporated (1975). 3) T. Harper and T. Inouye, GAMANL, a computer program applying fourier transforms to the analysis of gamma spectral data, Massachusetts Institute of Technology, MIT-3944-2 (1968). 4) L. Varnell and J. Trischuk, Nucl. Instr. and Meth. 76 (1969) 109. 5) j. Kern, Nucl. Instr. and Meth. 76 (1970) 233. 6) j.T. Routti and S.G. Prussin, Nucl. Instr. and Meth. 72 (1969) 125. 7) p. R. Bevington, Data reduction and error analysis f o r the physical sciences (McGraw-Hill, New York, 1969).
INSTRUMENTS
AND
METHODS
I37
(t976)
I73-I78;
©
NORTH-HOLLAND
PUBLISHING
CO.
A L G O R I T H M FOR C O M P U T E R P R O C E S S I N G A N D ANALYSIS OF G A M M A RAY SPECTRA FOR LARGE-SCALE L A B O R A T O R Y O P E R A T I O N S B R A D L E Y A. R O S C O E and A. K E I T H F U R R
Nuclear Science and Engineering, Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061, U.S.A. Received 11 May 1976 For large scale neutron activation analysis laboratories, computer analysis o f data is necessary. This work discusses major features o f the analytical program which has been developed and is currently being used at Virginia Polytechnic Institute and State University. Specific attention is drawn to fitting function choice, controlled fitting techniques, and methods o f error analysis. In order to test the program, samples o f standard reference material from the National Bureau o f Standards were obtained and analyzed. The results o f this analysis are given and demonstrate that the code is effectively doing its intended task.
1. Introduction With the acceptance of neutron activation analysis (NAA) as an accurate method of elemental analysis, several laboratories, such as the one at Virginia Polytechnic Institute and State University (VPI & SU), have made their services widely available. Analysis has been carried out on automobile emissions, organic materials, soils, industrial effluents, and a multitude of other types of samples in efforts to establish the environmental impact of today's society on the ecology. Another important application is determining trace element concentrations in an attempt to find correlations between these trace elements and the sample's general characteristics, i.e. leanness of beef, susceptibility of teeth to decay, etc. As a semi-commercial service, the N A A Laboratory at VPI & SU analyzed over 3000 samples in 1975. This service was extended primarily to on-campus groups but was also available to other institutions on a timeavailable basis. For each sample analyzed, there were at least three gamma spectra associated with it such that over 9000 individual g a m m a spectra were analyzed in that year. For this reason, the need of an effective computer code to analyze g a m m a spectra can be realized. There have been several programs for g a m m a spectra analysis developed over the years1-3); however, when development of our current program started at VPI & SU in 1972, only relatively unsophisticated ones were readily available. The program currently being used at VPI & SU has developed into a very powerful and inexpensive analytical tool for spectra analysis. Typically, the cost of a single element identification and quantitative determination is on the order of 2-3 cents for computer time. It runs on an IBM 370,
mod 158 computer and requires approximately 120k for compilation and execution. It is written in modular form so that, if necessary, it could be adapted to run on a much smaller system. 2. The fitting function The selection of a proper fitting function is of primary importance for data evaluation. In the selection of this function, one must consider all the limiting factors in the entire counting system from detector to memory. The first consideration is to the energy range to be observed and the number of data channels available with which to store data. In fitting a function to data using a least-squares method, there must be at least one more data point than parameter in the function. A fine point of least-squares fitting techniques that should be noted is that in order to obtain a good quality fit, the number of points in the actual peak, not including the background portion of the data, should be greater than the number of parameters in the fitting equation. For the counting system currently used at VPI & SU, an energy range from 60 keV to 3200 keV is observed using 4096 channels. This configuration yields about 0.8 keV per channel for acquired data. With peak width (full width at half maximum) of about 2.0 keV and a total width at 0.1 maximum of about 4 keV, an average peak will only have about 7 or 8 data points clearly distinguishable from background. Also to be considered in choosing the fitting function is the actual shape of the photopeak. Ideally, a photopeak should take the form of a Gaussian. However, due to imperfections in detectors and amplifiers, the Gaussian is always distorted. The most prominent cause of the distortion is produced by the amplifier
174
B R A D L E Y A. R O S C O E A N D A. K E I T H F U R R
giving the Gaussian photopeak an exponential tail on the low energy end due to inadequate base-line restoration at elevated counting rates. In order to reduce this effect, it is necessary to use an amplifier with superior baseline restoration such as the O R T E C model 472 used at VPI & SU. Possible fitting functions suggested by other researchers have included Gaussians plus: 1) an exponential tail on the low end4), 2) exponential tails on both endsS), or 3) exponential tails on both ends distorted by a power series6). In all the above functions, a quadratic was included in the function to account for background. It should be noted that in the formulation of these functional shapes, the experimenters were using data with at least 20 points and frequently more in the photopeak. In the analytical program discussed here, the exponential tail is ignored since an amplifier with good baseline restoration is used, moderate count rates are employed, and 8 data points will not show the exponential tail structure prominently. Thus the fitting function used is as simple as possible, i.e.,
established. The fitting function used at VP1 & SU is obviously non-linear, thus a method of non-linear leastsquares fitting is used. For this purpose, the "algorithm of Marquardt"v), which uses an analytical type gradient search, is used with extremely good results (see figs. 1 and 2). A typical peak will require only 3 iterations to obtain a satisfactory fit but up to 8 are allowed in attempting to fit the data. Since more than one solution is possible with a non-linear least-squares fit, it is necessary to control the parameters during the fit so that the best solution is obtained; thus the lowest chi-square is obtained instead of a minimum in the chi-square hypersurfaceT). This is achieved by making good first guesses for parameters and by preventing them from wandering from established limits. Initial guesses for parameters must be picked as closely as possible to the nominally correct solution. The peak location and width are the most important parameters that need to have good first guesses. In the 8
8
Y(1) = A(1) exp -
\
A(5)
] _]
+ A(4) + A(2) X(I) + A(6) X(I) 2 ,
~
BackgrOund+ Bromine Peak
|
(1)
where Y(I) is the number of counts in the X(1) channel; A(2), A(4), and A(6) represent the quadratic background; and A(1), A(3), and A(5) are the height, location, and width of the peak. In the event of a doublet being present (two peaks so close together that they overlap), an additional term is added to the functional equation to compensate for the interference,
Z
.
o'
.
xooLoo
.
,
I
.
.
Loos.~o
.
:u:..
i rol.o,°o P: ,k r
Lo~oo
I
xozx oo
I
Lo~7. ~
r
z~z. oo
NORPIRLIZED CHRNNF--L LOCRT[ON
~.m
J
l . o ~ oo
Fig. 1. Least-squares fit to a bromine peak at 698 keV using 6 variables [eq. (1)] and 29 data points. 8
A(7) exp
- 2
A(5)
'
(2) 8.
where A(7) and A(8) are the height and differential location of the interference peak. In order to minimize the number of free parameters, it is assumed that the width of the two overlapping peaks will be identical. It has also been found that having A(3) and A(8) coupled together for the location of the second peak will aid in the rate of convergence of the routine.
8 Background/
o~ 3. Method of curve fitting
In photopeak analysis, it is necessary to fit the acquired data to the desired fitting function so that the location and area of the peak can be accurately
~aoa
L'o~_00 ~ 0
NORHRLIZE0
~00
~co
~'a~o0
CHFiNNIZL L O C R T I O N
LIOLZ.EIO
Fig. 2. Least-squares fit to an iron peak at 1099.2 keV and an unknown interference using 8 variables [eqs. (1) and (2)] and l0 data points.
ALGORITHM
175
FOR C O M P U T E R P R O C E S S I N G
analysis program at VPI & SU, the program first calculates a channel vs energy and width vs energy calibration for each spectrum so that a good first guess can be made for the location and width of the peak. For the channel vs energy calibration, two known peaks are located such that a linear correlation may be made. One of the peaks is from 241Am at 59.543 keV which is present in every gamma spectrum from a source mounted near the detector. Due to the numerous possible interferences at very low energies from many of the rare earths, it has been decided not to attempt to use data below approximately 70 keV. Hence, the presence of the 2 4 1 A m peak does not adversely affect any data of interest. For a second calibration point, one of the following peaks is located: :4Na at 2754.1 keV, :SA1 at 1778.8 keV, or ~°K at 1460.75 keV, the choice being dependent upon the decay interval. In almost all cases, one of the above peaks will be present, such that the linear calibration can be made. Minor deviations from a linear fit can be accommodated by slight modifications of the library photopeak energies if necessary. The width vs energy calibration is made from the determined width of either the Na, AI, or K peaks. It has been experimentally verified that a linear width vs energy curve will vary only in the intercept of the linear equation and not in the slope term (see fig. 3). Thus, since the slope of this calibration is known, a new intercept can easily be found for the calibration from the one data point. Between iterations, the width, height, and location parameters are checked to insure that they have not wandered from pre-established limits. If one has deviated too much, it is reset to a more reasonable
value and the interation continued such that the desired solution may be obtained. When looking for a specific photopeak, the program will fit the data to a gaussian plus a quadratic background. If in the final fit, the width is too wide (see fig. 4), the data will be re-fitted using the equation for a doublet, thus accounting for interference peaks. The method of controlled fit to the non-linear equation has been proven effective in the analysis of complex spectra. It is quick, accurate, and needed in the computer analysis of data for large scale analysis operations. 4. Error estimates
As in any experimental procedure, it is desirable to have a method of error analysis so that a confidence level in results can be realized. In the analysis of gamma spectra, two different types of error are considered, statistical error and fitting error, both of which can be calculated. By examining the value of these two quantities, the interpreter of the data can achieve a level of confidence in the final data. The statistical error is obtained from basic counting theory which states that the standard deviation of the number of counts is equal to the square root of the number of counts. Since peak areas are calculated over several channels, a slight variation is used such that the statistical error, equivalent to one standard deviation, is calculated from S.E. = ~/[A(peak) + A(bkg)] A (peak)
x
loo Vo,
(3)
where A(peak) and A(bkg) are the areas of the
D ipl"~
Oe~dTm i e~42%~ UpperLimit
(I wtq
0-_ i -
~
^
•
~ . _ ~DeadTi. . . .
ENERGY (XIDw Jl
Fig. 3. Peak width times. Data taken represented by ' © ' peak width is the channels.
vs energy (keV) for different system dead at 42% dead time and 2% dead time are and ' × ', respectively. On the ordinate, the Gaussian width o f the peak measured in
5.oo
"~.m
n~rn
ENE'OO FIGY C~ID,, @
Fig. 4. G r a p h of peak width vs energy (keV). The lower and upper limits are used as trips by the computer p r o g r a m to aid in the peak fitting process. If during the fit the new calculated width of the peak falls outside the limits, it is reset to the calibration line and fitting continued with this new value.
176
BRADLEY
A. R O S C O E
AND
A. K E I T H
FURR
TABLE 1A Orchard leaves.
Element
Ca K Mg Fe Mn
Na Zn As Cu Rb Hg Br CI Cr Co S
Gamma energy (keV)
NBS a cone. (ppm)
3084.4 1524.7 1014.4 1099.2 1291.5 846.6 1811.2 2112.6 1368.5 2754.1 438.9 1115.5 559.1 657.1 511.0 1039.0 1077.0 77.6 698.3 776.5 1642.4 2167.5 320.1 1173.2 1332.5 3102.4
20900+300 14700±300 6200±200 300 ± 20
Cone. (ppm)
91 ±
4
82±
6
25±
3
14±
2
12±
1
12± 1 0.155±0.015 (10) (700) (2.3) (0.2) (2300)
21600 15390 6174 279 274 93.4 91.2 95.2 98.6 91.2 17.6(N) 24.2 14.9 13.8 12.9 9.1 12.0 0.188 8.6 8.6 695 732 2.17 O. 18 0.21 2744 (N)
Sample 1 Count error (%)
Fit error (%)
1.3 0.7 1.6 7.6 10.4 0.5 1.5 2.2 0.8 1.1 52.4 5.8 0.6 0.1 3.8 22.6 7.8 9.7 1.5 0.8 2.4 2.0 I 1.1 18.4 16.7 58.8
67.4 2.3 8.7 184.4 40.9 0.9 5.0 0.3 1.0 14.5 3.3 169.6 0.I 2.6 5.1 34.6 130.4 0.3 2.0 1.9 6.8 65.6 18.2 65.1 217.4 69.0
Cone. (ppm) 21200 14540 6436 294 369 92.4 89.6 91.3 65.1 68.2 23.8 24.2 14.4 13.4 16.8 (W) 7.6 (N) 12.0 0.175 (N) 8.2 8.4 709 688 2.9 0.22 0.26 4818 (N)
Sample 2 Count error (%)
Fit error (%)
1.4 0.7 1.7 6.7 7.9 0.5 1.6 2.4 0.9 1.2 47.6 5.7 0.5 2.2 3.2 25.5 7.3 8.5 1.5 0.7 2.5 2.2 8.6 14.4 14.4 34.7
26.5 ! .3 10.1 44.3 37.7 0.5 21.0 32.2 0.9 10.2 128.6 75.3 0.1 3.3 0.8 3.8 221.9 0.2 0.9 2.1 5.8 31.4 7.0 0.9 172.6 192.9
a "values in parentheses are uncertl.eU. Gauss±an peak and quadratic background over the i n t e r v a l o f t h e p e a k ' s full w i d t h a t h a l f m a x i m u m . The method of determining the error due to the least-squares fitting process has only been developed for l i n e a r cases. F o r n o n - l i n e a r t e c h n i q u e s , it is a s s u m e d t h a t t h e e r r o r c a n be c a l c u l a t e d in a s i m i l a r m a n n e r . I n l i n e a r l e a s t - s q u a r e s fitting, a s q u a r e s y m m e t r i c m a t r i x is g e n e r a t e d t h a t c o n t a i n s t h e w e i g h t e d d e r i v a t i v e s o f all t h e p a r a m e t e r s in t h e function. Taking the inverse of this matrix forms an e r r o r m a t r i x w h i c h c o n t a i n s t h e v a r i a n c e s a n d cov a r i a n c e s o f all t h e f u n c t i o n a l p a r a m e t e r s . T h e determinant of this matrix represents the overall v a r i a n c e o f t h e fit. By d i r e c t a n a l o g y , t h e v a r i a n c e o f t h e fit f o r t h e n o n - l i n e a r case is t h e s a m e d e t e r m i n a n t . T h u s , since t h e s t a n d a r d d e v i a t i o n is t h e s q u a r e r o o t o f t h e v a r i a n c e , t h e fit e r r o r is c a l c u l a t e d f r o m : F.E. =
fNSD " "
] '/N
IA(od
x 100 °/ /o,
(4)
where S.D. = s t a n d a r d d e v i a t i o n o f t h e fit, N = number of variable parameters, and A (i) = p a r a m e t e r n u m b e r i. It s h o u l d b e n o t e d t h a t t h i s m e t h o d o f c a l c u l a t i n g e r r o r is a q u a l i t a t i v e m e a s u r e o f t h e g o o d n e s s o f fit a n d is not absolute. T h e s t a t i s t i c a l e r r o r a n d t h e fitting e r r o r p r o v i d e useful i n f o r m a t i o n t o t h e e x p e r i m e n t e r a n d a l l o w h i m t o e s t a b l i s h l i m i t s b e l o w w h i c h h e will c o n s i d e r d a t a t o b e a c c e p t a b l e . H o w e v e r , in s o m e i n s t a n c e s , r e s e a r c h w o r k e r s a r e i n t e r e s t e d in w h e t h e r t h e r e is a n y p o s s i b i l i t y a t all o f a n e l e m e n t b e i n g p r e s e n t in t h e s a m p l e . T h e r e fore by providing the raw values of the two error e s t i m a t e s , t h e c h o i c e o f w h e t h e r t o use d a t a is t r a n s fered to the researcher and has not been predetermined by the program.
5. Calculation of elemental concentration T o m i n i m i z e e r r o r s in c a l c u l a t i n g e l e m e n t a l c o n c e n t r a t i o n , s t a n d a r d s h a v e b e e n a n a l y z e d o n all
ALGORITHM
177
FOR C O M P U T E R P R O C E S S I N G
TABLE 1B Orchard leaves.
Element
G a m m a energy (keV)
Dy La
Mo Sb
Sm Ba Cr Sc Cs W Au
Conc. (ppm)
94.6 108.2 328.7 487.0 815.8 1596.2 141.0 564.0 645.8 1691.0 69.6 103.1 216.0 496.3 320.1 889.0 1120.5 795.7 479.4 685.7 411.7
0.14 0.14 0.84 0.91 0.92 0.91 0.19 (N) 2.52 4.6 3.5 0.64 0.48 26. 38. 1.25 0.038 0.035 0.087 0.139 0.097 0.00358
Sample 1 Count error (%)
Fit error (%)
31.9 6.0 2.8 1.8 3.8 2.0 35.9 0.6 12.5 5.0 3.5 0.8 21.0 14.9 14.1 5.5 6.2 32.9 42.8 38.9 12.1
Conc. (ppm)
3.1 0.9 0.7 3.7 2.9 4.4 0.0 0.0 36.1 205.0 0.3 0.1 6.3 24.1 5.7 35.1 223.2 22.4 3.7 19.1 1.8
Sample 2 Count error
0.179 0.156 0.88 0.92 0.83 0.91 0.202 (N) 3.2 4.0 (N) 3.9 0.60 0.44 36. 39. 1.04 0.0038 0.0043 0.067 0.103 0.16 0.00398
(%)
Fit error
(%)
22.1 6.4 2.5 1.7 3.6 1.8 30.8 0.5 12.7 4.3 3.4 0.8 15.6 13.2 14.4 5.1 5.0 43.4 46.5 28.9 9.7
1.2 2.5 0.7 3.0 3.5 4.3 0.4 0.1 I.I 189.4 0.2 0.1 6.2 6.3 2.8 30.1 47.9 22.5 1.4 4.6 2.8
Sample 2 Count error
Fit error
(%)
(%)
1.0 0.2 0.2 8.5 10.2 0.6 1.4 27.1 2.6 5.8 1.3 5.4 8.4 16.5 34.3 0.8 0.8 11.8 10.7 7.1 5.4
4.8 0.2 0.3 211.9 129.7 0.1 7.5 1.5 22.5 38.8 1.3 14.5 8.5 23.4 37.2 8.5 1.8 319.7 187.6 0.9 0.2
TABLE 2 Bovine liver.
Element
K Na Fe Cu Zn Rb Mn
Se Ca CI Co Mg Mo
Gamma energy (keV)
1524.7 1368.5 2754.1 1099.2 1291.5 511.0 1039.0 438.9 1115.5 1077.0 846.6 1811.2 2112.6 161.6 3085.0 1642.4 2167.5 1173.2 1332.5 1014.4 141.0
NBS a conc. (ppm)
9700 2430
4-600 4- 136
270
4- 20
193
4- 10
130
4- 10
18.3410.3 4-
I I
1.1 + 0.1 (123) (2600) (0.18) (605) (3.2)
a Values in parentheses are uncertified.
Conc. (ppm)
Sample 1 Count error (%)
10320 2003 2027 234 315 275 (W) 192 89 103 19.3 11.2 11.0 8.9 1.8 74 (N) 2891 2806 0.39 0.38 795 2.6
1.1 0.8 0.3 9.8 10.4 0.6 1.5 24.7 3.0 6.5 1.4 6.2 11.8 26.7 66.2 0.9 0.9 11.9 11.8 6.2 6.0
Fit error (%)
3.8 0.2 0.5 462.2 466.5 0.8 11.0 2.6 250.4 118.8 1.4 6.3 0.5 4.4 34.0 5.8 10.6 118.1 225.2 5.9 0.1
Conc. (ppm)
10080 1920 1955 224 246 264 (W) 181 67 107 18.2 10.3 10.6 11.7 1.6 125 2640 2666 0.30 0.35 573 (N) 2.4
178
BRADLEY A. ROSCOE AND A. KEITH FURR
elements that are sought by the VPI & SU laboratory. In this manner, inaccuracies in branching ratios, cross section, neutron spectrum, and natural abundances have no effect on final data evaluation. Elemental concentrations are then calculated by standard methods from the peak area and activation times. In order to minimize computer time and eliminate confusing answers appearing in the output, numerous flags have been incorporated in the program so that the program only looks for individual elements under conditions predetermined by the activation, decay, and counting times. Although certain rough rules have been incorporated in setting these flags, empirical factors, such as the presence of interferences, have proven more important and hence the flags have been determined laboriously, element by element.
6. Program testing In order to test the accuracy of VPI & S U ' s analysis program, samples of standard reference material, such as bovine liver and orchard leaves, were obtained from the National Bureau o f Standards. In the data presented here, two samples of each material were prepared and analyzed using the same procedures that any other samples undergoing analysis would receive. In this manner, the overall accuracy of the laboratory was tested. The results of the analysis of the reference materials are given in tables 1A, 1B, and 2. Tables 1A and 2 give comparison of certified and uncertified NBS concentrations while table 1B gives elemental concentrations contained in the sample that were found at VPI & SU and not reported by the NBS. In the tables, an (N) or (W) represents a narrow or wide peak as determined by the c o m p u t e r ' s width vs energy calibration. As should be expected, the concentration for a narrow peak tends to be low while the value for a wide peak tends to be high. It should also be noted
again that the fit error is only a relative measure and not absolute. As can be seen from the data, the results obtained correlate very well with NBS certified and uncertified values except for the cobalt in bovine liver and sulfur in orchard leaves. Since N A A is not very sensitive to sulfur, the result is not surprising; however, for sulfur, an order o f magnitude answer is reasonable. The cobalt in the bovine liver disagrees by about a factor of two. We believe that the uncertified NBS value is off by this factor due to our extreme sensitivity to cobalt and the fact that the cobalt data for the orchard leaves compared reasonably well. Another factor providing confidence in our cobalt value is that cobalt is our basic standard for determining the thermal neutron flux in the VPI & SU reactor. F r o m the results obtained from the NBS reference material, it can be seen that the analysis program employed at VPI & SU is doing its intended task very effectively. It has vastly reduced data analysis time, decreased errors and increased the a m o u n t of data that can be processed at the laboratory.
References 1) T. Meyers and M. Denies, Longterm and Peakscan: neutron activation analysis computer programs, Michigan: Museum of Anthropology (1972). 2) Radionuclide program: Peakfind and Isoquan, Ortec Incorporated (1975). 3) T. Harper and T. Inouye, GAMANL, a computer program applying fourier transforms to the analysis of gamma spectral data, Massachusetts Institute of Technology, MIT-3944-2 (1968). 4) L. Varnell and J. Trischuk, Nucl. Instr. and Meth. 76 (1969) 109. 5) j. Kern, Nucl. Instr. and Meth. 76 (1970) 233. 6) j.T. Routti and S.G. Prussin, Nucl. Instr. and Meth. 72 (1969) 125. 7) p. R. Bevington, Data reduction and error analysis f o r the physical sciences (McGraw-Hill, New York, 1969).