Nuclear Instruments and Methods 190 (1981) 89-99 North-Holland Publishing Company
89
SAMPO80: AN ACCURATE GAMMA SPECTRUM ANALYSIS METHOD FOR MINICOMPUTERS Markku J. KOSKELO, Pertti A. AARNIO and Jorma T. ROUTTI Helsinki University of Technology, Department of TechnicalPhysics, SF-02150 Espoo 15, Finland Received 28 January 1981 and in revised form 9 June 1981
In this paper we report a minicomputer version of the well known gamma spectrum analysis program SAMPO. The algorithms an the general structure of this new version called SAMPO80 are explained. We also report results of statistical tests showing that the peak locations and areas given by the new version are both consistent and accurate.
1. Introduction The SAMPO family of gamma spectrum analysis programs is a result of systematic program development at the Helsinki University of Technology. The starting point for our work is the original version coded by Routti at the Lawrence Berkeley Laboratory of the University of California in 1969 [ 1,2]. The original SAMPO was designed to meet the stringent requirements of nuclear spectroscopy work and it features algorithms for automatic peak search and peak fitting as well as peak shape, energy, and efficiency calibration procedures. The original version was coded for a CDC-6600 computer using Fortran IV. At the Helsinki University of Technology the original SAMPO is run on a UNIVAC 1108. The code was soon found out to be applicable also in many routine applications and it has gained a large number of users throughout the world [3]. The original code is rather time consuming, however, and with modifications and improvements a reduction o f the execution times to about one tenth to one fortieth of the original is possible without seriously affecting the accuracy of the results. Our first modification to achieve this was the SAMPO76 which was developed in 1976 in collal~oration with the Radiation Protection Group of the European Organization for Nuclear Research, CERN [4]. The analysis methods employed in this version are simplified from those in the original code. Good precision is still obtained, however, and a number of features have been added. These include a versatile File management system for handling large amounts of data and two identification algorithms. The peak 0029-554X/81/0000-0000/$02.75 © 1981 North-Holland
shape calibration procedure is not included. This version is run both on the CDC-7600 of CERN and our UNIVAC 1108 and is coded in Fortran 5. Our latest version of the SAMPO is the SAMPO80 which is discussed in this paper. SAMPO80 is coded in Fortran 5 for a Data General Nova 2 minicomputer with a 32 kiloword central memory and a supporting disk unit. A PDP-11 minicomputer version of the code also exists. The SAMPO80 consists of three separate parts, the shape calibration part SAMPOSHAPE, the peak search and fitting part SAMPOFIT, and the nuclide identification part SAMPOID. The basic algorithms are imporved and modified from those of the previous SAMPO versions. 2. The mathematical algorithms A more thorough explanation of the algorithms used in SAMPO80 can be found in ref. 5. Thus only the main features of it and particularly where it differs from the original SAMPO are explained here. The most characteristic feature of the original SAMPO is the adoption of a separate peak shape calibration procedure. The peak shape function has been chosen in such a way that its parameters vary smoothly with energy. Hence the shape parameters are determined separately only for a few isolated peaks and interpolated from the values thus found for the rest of the peaks. The shape calibration procedure only needs to be done once for a given measurement configuration ff the conditions remain reasonably stable. This feature is therefore especially well suited for routine applications where the speed of the algorithm is also of importance.
M.J. Koskelo et al. / SAMP080
90 -
K
[] =CL-PFIRFIMETER )( =CH-PFIRFIMETER
._d CI2 cv t.l.J p--o
~(E Old ClZ~
g N
J
o
~00
8.00
IG
O0
CHAMMEL
24,00
32,00
MUMBER
40 O0
-!0'
Fig. 1. The dependence of the peak shape parameters on the channel number for the gamma spectrum of the IAEA G-2 source no. 2. In SAMPO80 the peak shape calibration procedure has been separated into part I, SAMPOSHAPE, of the code. The functional representation for the peak is still assumed to be a Gaussian with exponential tails on both sides. The exponentials are joined to the Gaussian in such a way that both the function and its first derivative are continuous. The peak shape is characterized by 7 parameters, two background parameters, the peak height and its location, the peak width (CW) and the distances from the peak centroid to the starting point of the exponentials on either side (CL and CH). The last three are used as calibration data for subsequent analyses since they vary rather smoothly with energy (see fig. 1 for an exampie). Fig. 1 as well as the printouts in this paper are
all calculated for the IAEA intercomparison G-2 source no. 2. The complete spectrum is shown in fig. 2. The minimization of the least-squares expression to solve for the peak parameters is performed by a subroutine package with an iterative gradient algorithm utilizing the variable metric method [6]. This package has been taken from the original SAMPO as such without any significant changes. The minimization is thus terminated if all components in the next step change by less than 10 -8 , or if four succeeding values of X2 are the same, or if 100 iterations have been completed. An example of the printout produced for the shape calibration of one peak is shown in fig. 3. The upper part shows numerical information on the preliminary linear fit and the subsequent non-linear fit. The lower part is a semilogarithmic printer plot of the fit. For checking that the performed shape calibration is acceptable a few error parameters are also calculated and printed out. The original goodness-of-fit indicators are the values CHI/DEG.FREE and SIGMA. The former is the conventional X2 per degrees of freedom. The latter is its square root. For SAMPO80 an additional parameter called ERROR CORRELATION in the printout has been developed [7]. As a rule of thumb, SIGMA values below 5 and ERROR CORRELATION values from - 1 to +1 can be considered to be acceptable. Another characteristic feature of the SAMPO is the use of the smoothed second difference to find the peaks in the spectrum. The method is particularly good in detecting small peaks on a high or low background. Multiplets are detected with a reasonable accuracy. The original peak search routine of the SAMPO is based on the work by Mariscotti [8] and defines the tentative location for a peak as the integer channel where the smoothed second difference ssi = ddi/sdi
(1)
reaches a minimum while the final location is found in the non-linear fitting. The ddi and the sdi are defined as the generalized second difference expression k dd/= ~ c/ni+I , /=-k and its standard deviation
(2)
M.J. Koskelo et al. /SAMP080
91
~o
_j~= Ld :-
C_] r
7.
73 C) ~
~E
°I CD
I
I 4.~I
.01
I
I 6. 41
I
i 12.G1
I
1 --I IG, BI
I 21o01
CHAHHEL
I --
Io 25 21
HUMBER
I
I 2S.41
I 33.$I
i
: 3~. 81
7---I
42.01
.I0 Z
Fig. 2. The g a m m a s p e c t r u m o f t h e IAEA G-2 source no. 2.
sdi =
cTni+/
1
(3)
.
The counts per channel are denoted by n i and the summation is done over 2k + 1 channels where the number k depends on the coefficients c/. The coefficients are defined in SAMPO80 as Cj
- 100(]2 - °~v) exp Oav 2
[:1 ~ 2Oar
,
(4)
where Oar = fwhmav/2.355 of the peaks in the search interval, co is thus - 1 0 0 and the set is terminated at k, where ICk+l I < 0.01 Icol- The second coefficient is then adjusted so that the sum of the coefficients is zero. The normalization of the coefficient c/is different in SAMPO and SAMPOS0. In SAMPO80 the peak search algorithm defines the centroid channel of the found peak with the expression i • ssi
i
ch =
~ ssi i
(5)
The summing is done over the channels i where the significance value ssi is negative, i.e. the location is found as the weighted average of those channels. The method has been adopted for our code from the ideas presented in e.g. refs. 9 and i0. This modification is of crucial importance for linear SAMPO80 where the peak location parameter is not allowed to move during the fitting. In SAMPO80 as in the original SAMPO there are two checks to make sure that only the right peaks are accepted for fitting. One is based on the absolute value of the significance ssi and the other on the width of the peak. The effect of these acceptance tests can easily be seen from the printout of the peak search routine as depicted in fig. 4. The different method of finding the peak location in SAMPO80 as compared to the original SAMPO can also be seen from this printout. In the original version only the approximate integer peak channel was given in this printout. In SAMPOS0 the channel is given as the exact location of the peak and the value given at this stage is f'mal. The peak areas are calculated by fitting the precali-
1076.00
lla
5, SS
1067.00
1088.00
.DO
1.00
.IBSE 04 -.14SE Ol
=
FINRL GRRDIENT
.BTE 04 .B4E 04 .BIE 04
.2 -.7 I.I
lOBG 1087 1088
I
3°663 S.330
I
I I
+/~
¢I+/~/-
.003
.003
.028 .IBB
FIT
•
=
3,G72
SUMCRLC=
=
SIGMA=
÷
4..
(. +.
.I12E-03
-.185E 02
.198E 01
.235E 04
(
=
lAND4
)
=
.232828E 07 ERRORCORRELRTION=
.AND+
IRND.
-.174
$
=
IAND+RND.
-.841E 02
.lOBE 04
.I#TE 03
+,
YMA~= .478E OS
.5SIE O0
.231E 01
-.47SE 04
IIIIllIlIIIIIl
-.267E 01
.191E Ol
.739E 04
I I I I I I I I I I I I I I S E M I L D G R R I I H M I CS C R L E
-.834E-02
OG
,307E-01
Fig. 3. T h e p r i n t o u t o f t h e s h a p e c a l i b r a t i o n f o r o n e o f t h e p e a k s f r o m t h e g a m m a s p e c t r u m o f t h e s o u r c e n o . 2
1,978
10~B.122
CP =
CL = CH = CW =
SHRDE CRLIBRRTION RESULTS CHI/DEG,FREE = . 1 3 S E02
5YMBOLSs SUMDRTR=
I I I)+ $
I I
I I
I I
[ I
I [
. 2 4 G ED0
~/,3102SE-OZ -.309E 03 .4GGE
.BODE 04 .(
-.IGBE-02
,20223E 0 3 . 1 l I E 05
I = CONTINUUM = DRTR .232B47E 07
.~OE 0 4
-1.0
05 0S 05 05 05 05 OS 05 05 OS 05 OS 04 04 04
lOBS
1079 I080 1081 1082 1083 1084
1078
IqE 13E 13E 13E 12E 12E 12E lIE lIE IIE IIE IOE 99E 96E ,53E
4.3 .9 -3. B -3.6 5.1 .4 -2.S -5.3 2.S Z.B 1.1 -.9 -2.4 -.4 .~
1070 1071 1072 1073 1074 1075 1076 1077
CHRM, BEg, CUMTINUU~ 1067 -3.2 15E OS 1068 -.3 I SE OS 10S9 1.G 14E O S
= =
FI~RL F OR CHISBR FI~RL PRRRMETERS
TR[E~ IO ~OUE MORE THR~ 12TIMES TRIED TO MOVE MORE THR~ 12TIMES
I~ITIRL F OR CHISQR= I~ITIRL GRRDIE~T =
-m-
4-•
FINRL F OR CHISQR = .131GTE04 + / .5S22SE-02 FINRL PRRRMETERS = .464E OG .IOBE 04 .19BE Ol .119E OS - . 3 7 0 E 03 FI~RL GRflDIE~T .870E-04 -.21SE 0 2 .2BoE 02 . B 2 1 E - O ~ -.129E-02 TRILED GRUSS.?flRAM,= CO~T,CONST. CONT.SLO?E PERK HEIGHT GRUSS.WIDTH LOW~ Tfl[L HIGH TAIL PERK CE~TROID .1990 O1 .lORD 04 [~ITIRL GUESSES = .9880 Oq - . 3 7 0 0 03 .4GODOG .1990 01 .1990 Ol
CRLI3RRTION FIT GRUSS[R~ U I T H EXPONENTIRL TRILS CHRNNEL MUMBER= 1078 GAUSSIRN FIT ?RRRM.= ?EBK HEIGHT ?ERK CENTR. GRUS.WIDTH CONT.CONST. CONT.SLOPE INITIAL GUESSES = .4670 OG .1080 04 ,1800 01 .1110 OS -.2BBD 03 I~ITIRL F OR CHISQR= . 4 9 7 E OS I N I T I f l L GRflDIENT 5 9 3 E DO - . 9 7 4 E 05 -,41SE OG - . G I 3 E Ol .241E O1 TRIED TO MOVE MORE T H f l N - ' 1 2 T I M E S TRIED TO MOVE MORE THRB 12TIMES
SHRTEDO
~E~T CONTROL CRR3,
DRT~ .143E ,144E .147E .ISSE .IGPE .209E .33SE .G94E .14GE .Z73E .40SE .478E .~3SE .308E .172E .77SE .315E .IS3E .104E .911E .843E .821E
OS OS 0S OS OS OS OS OS OG OG OG 06 OG OS OG OS OS OS OS 04 04 04
FIT • 146E OS • 14SE 0S • 14SE OS • 145E OS • ISGE OS 214E OS 341E OS 670E OS 14GE OG 274E DG 40BE DG 477E OG 433E DG 308E OG 172E OG 2 8 2 E OS 315E OS • 1G2E OG • lOSE 0S • 909E 04 • 849E 04 • 811E D4
co
M.J. Koskelo et al. /SAMP080 NEXT CONTROL C A R D ,
PEAKFIND
.00
40
408S
SEARCH BETWEEN CHANNELS
RUERRGE PEAK WIDTH (SIGMA) = 2.24 GENERALIZED SECOND DIFFERENCE RRE
PERK CHRNNEL
I 97. S63 2 I06.670 3 IGI.627 4 246.201 5 370.6GI G 414.223 7 487.25~ @ GPG.SPS 5 723.097 10 837.377 II 894.319 12 917.376 13 1078.130 14 1162.302 IS 1322.606 16 1 6 4 6 . 8 6 3 17 1766.268 18 1844.472 19 2626.842 20 2768.191 21 2 5 3 3 . 3 1 4 22 3 3 7 6 . 8 1 9 23 38S8.G02 24 4 0 3 4 . t 7 2
PERK ENERGY
31.9 34.9 $3.0 81.0 122.1 136.S IG0.7 223.4 238.8 276.G 296.S 303.1 366.3 384.1 437.1 Sll.3 $83.S 609.7 836.4 511.8 569.6 II16.0 fZ~S.1 1333.1
I0.00
AND THE COEFFICIENTS OF THE -100.00 - 7 1 . 8 6 - 1 3 . 4 8 32.47
SIGNIFICANCE LIMIT FOR 7HE SEARCH =
[
.00
93
I0.00
SIGNIF C H E C K - I OF PERK SIGMIF
432.739 396.272 126.62S 946.628 IS24.198 489.G14 74.329 48.170 26.274 348.675 IO.GG2 SSG. I@3 566.126 321.937 13.049 107.413 14.695 21.216 672.328 12.220 11.233 S23.403 22.78S 11.S3S
AND FOR FITTING =
I00.00
44.43
32.89
IS,
.00
17.00
G.S9
Ss34
1.97
i00.00
CHECK-2 RCCEPT SHRPE CHANNELSMUM
CHEC CNEC
57 10G 161 2O,G 370 414
l 2 3 4 S G
837
7
517 1078 1162
8 9 10
IS4G
ll
2626
12
3376
13
SMAL SMAL SMRL
SMAL
SMAL SMAL SMAL SMRL SNflL
SMRL SMAL
Fig. 4. T h e printout of the peak search routine for the gamma spectrum of the source no. 2. Note that the printing and fitting thresholds have been set quite high to reduce the number of peaks to the most significant ones that are of interest in this case.
brated modified Gaussian to the data with a weighted least squares formula using a parabolic background. The fitting intervals are determined automatically by the program. Peaks separated by less than four times the average fwhm are fitted together. The minimization of the X2 is done with respect to the background parameters, the peak height and the peak location in the original SAMPO with the same non-linear least squares routine that is used for shape calibration. In SAMPO80 the peak centroid is obtained directly from the peak search routine and hence only the background and the peak height are fitted by a linear algorithm. The errors for the peak areas are calculatec[ from the errors given for the peak height parameters. In the non-linear original SAMPO the errors were obtained as the diagonal elements of the inverted matrix used in the minimization. In the linear algorithm the inverse of the matrix of the least squares minimization is not calculated explicitly. Rather the equation is solved by Gaussian elimination and back substitution. The diagonal elements of the would be inverse of this matrix are easily calculated, however, and this has been done in SAMPO80 in a new subroutine
called DIAG. The errors for the peak heights are obtained from these diagonal elements. The errors for the peak locations in the original SAMPO are obtained directly from the errors calculated for the centroid parameters. In SAMPO80 the error for the peak location is calculated in the peak search routine using the normal equation for the standard deviation of a mean and multiplying that by a factor of 1.3. The mean is the one calculated in eq. NUCLIDES 1
2
3
4
S
X u~ bd
2
L O n,, b.J Z bJ
3
b,c rr b.l p0 i-o "r p-
XX X XX
4 S
X X
S
X B
Fig. 5. The structure of the working matrix of the identification routine.
94
M.J. Koskelo et al. / SAMP080
(5). The multiplying factor was found necessary to give consistent and accurate values [11]. The calculated peak locations and areas are finally corrected with energy and efficiency calibration data, respectively, to yield peak energies and intensities. For the energy calibration linear interpolation on a linear scale and for the efficiency calibration linear interpolation on a log-log scale are used in the SAMPO80. In SAMPO also a functional fit on these calibration curves was possible. Calibration errors are added to the peak location and intensity errors to give the final result. The method of identification in the SAMPO80 is the technique of associated lines which is based on the work of Gurmink and Niday [ 12]. Several changes to the original method have been made, however. The forming of the matrix equation has been retained but IntercompBplson
IAEQ
G-2 s o u p c e
no.
the tests to ensure that only the right nuclides are accepted in the final printout are different. In our method a matrix which contains the branching intensities is formed by comparing the computed peak energies in the analyzed spectrum to those in the reference library. The tolerance value at which a nuclide is accepted for a tentative match can be given in the input data. The result is a matrix of the form depicted in fig. 5. This matrix is then searched for interference sets which are each solved separately using the weighted linear least squares technique. The weighting factors are the recipocals of the errors of the calculated peak intensities from the fitting program. The solution is accepted if the calculated activities are non-negative and their errors less than 50%. If this is not the case the nuclide failing in either one or
22.2.75
2
20000s
THE
FOLLOWING NUMBER II
PEAKS ARE A T T R I B U T E D TO CHANNEL ENERGY IS4G.BS SII.30
ISOTOPE
MR-22
THE
FOLLOWING NUMBER 12
PEAKS ARE A T T R I B U T E D TO ISOTOPE CHANNEL ENERGY 2S2G.84 B3S.3?
MN-S4
PEAKS ARE A T T R I B U T E D TO CHANNEL ENERGY 320.66 122.12 414.22 13G.$4
C0-$7
THE F O L L O W I N G NUMBER S E
ISOTOPE
THE FOLLOWING PEAKS ARE ATTRIBUTED TO ISOTOPE NUMBER CHANNEL ENERGY 13 3326.82 1116.01
ZN-GS
THE FOLLOWING PEAKS ARE ATTRIBUTED TO ISOTOPE NUMBER CHANNEL ENERGY 3 t G t . G3 $3.0S 4 24G.20 80.99 P 837.38 ZTG.61
BA-133
8
91P.38
303.09
9 I0
1078.13 1162.30
3SG.30 384.14
IAEA I n t e r c o m p a r l s o n G-Z source no. 2 2 7 . 2 . 2 9 20000s THE FOLLOWING ISOTOPES WERE IDENTIFIED NUMBER NUCLIDE COMF.URLUE SAMPLE ACT PC-ERRORSAT. ACT. PCI PCI I Nil-22 .592G 3.893E 04 4.22 7.240E I0 2 MN-S4 .5767 1.774E OG 3.70 I.IS3E 12
3
C0-S?
.9993
1.308E 0G
3.13
7.382E lI
4 S
ZN-GS BA-133
.9808 .8142
3.GSIE OG 1.506E 0S
2.10 2.07
1.8~IE 12 I . S P S E 13
TOTAL ACTIUITY ACCOUNTED FOR B.G79E OG PCI IHE FOLLOWING PEAKS WERE MOT IDENTIFIED NUMBER CHANNEL ENERGY INTENSITY 1 2
37.563 IOG.GPO
31.880 34. 885
PC-ERROR
I . 14E 08 t.51E 09
Fig. 6. The printout of the identification routine for the source no. 2.
3.83 3.91
M.J. Koskelo et aL / SAMP080
both of these respects is discarded and the interference set formed anew. A final list of results (see fig. 6) contains both information about which peak is assumed to be associated with which nuclide and the calculated activities for each nuclide along with the associated errors and the values for the confidence index.
3. Program structure A very typical feature of the SAMPO programs is the use of special SAMPO control commands to direct the execution of the program. These control commands are identified by English code words in
"•I1
SAMPOSHAPE
SHAPEIN r m - - - -
I I t I SAMPOFIT
7
,L i SAMPOID
Fig. 7. The structure of the SAMPO80 set of programs with their associated input and output files.
95
the main program. This feature is also used in the SAMPO80 although the number of recognized control words is more limited than in the previous versions. Each control word causes an appropriate subroutine call and the calculations are then done in the subroutines. The original SAMPO as it is installed on e.g. our UNIVAC 1108 consists of just main program and 27 subroutines. For all operations the whole program has to be initiated. In the SAMPO80 a different approach has been taken. The different tasks done by SAMPO are divided into three separate programs which are connected into a complete analysis system via common input and output files (see fig. 7). The different parts are called SAMPOSHAPE, SAMPOFIT, and SAMPOID and they are designed to do the peak shape calibration, the search and fitting, and the nuclide identification, respectively. The change into separate programs was necessitated by the limited central memory capacity of our minicomputer. This arrangement also makes the use of the system more flexible. The data and intermediate results are stored on a disk unit. Fig. 7 shows that I/O file structure. The measured spectrum can be read in from any input device. Usually it is first stored onto the disk as a file and then just the appropriate file name is given to the program. The same applies for the control command sequences and the gamma reference library, too. The output can be directed to any output device recognized by the operating system although usually the line printer is used. The SHAPEIN shape calibration data created by the SAMPOSHAPE can be stored onto a file and then transferred into the input sequence of SAMPOFIT using the editor of the operating system. The fitting results are always printed into a file with the same name (currently INFOUT) and are read directly into the SAMPOID as intput data. Each of the program parts is overlayed sufficiently to fit into the 32 kilowords of core we have available and their memory requirements can, if needed, be reduced even more. That would necessitate some structural changes in the subroutines, however, and also lead to increased computing times. Typical running times are currently approximately 25 s to find the peak shape parameters of one peak, 2 s per peak for the peak search and fitting run, and 3 s per identified nuclide for identification. The running times are dependent on the complexity of
96
M.J. Koskelo et al. /SAMP080
the problem, of course, and slightly longer times may be needed for more difficult cases. The different parts of the system are compiled by a Fortran 5 compiler. The code is basically Fortran IV however. The Fortran 5 compiler has been used because the overlay structure is easier to implement with it since it does not need any additions to the code. The code thus becomes more computer independent. The Fortran 5 compiler also makes a more efficient machine code and is able to utilize the hardware floating point unit we have available in our computer thus making the computations faster and more accurate.
4. The performance of the SAMPO80 4.1. The testing method The performance of the peak search routine and the peak fitting algorithm of SAMPO80 has been tested using the method suggested by Heydorn [13] and the results are compared with those obtained for the original SAMPO. The method is a standard procedure of statistics where the estimated standard deviations of replicate analytical results should fully account for the variability of the results themselves. The appropriate test is based on the Analysis of Precision and the set of six reflicate spectra provided by the IAEA in their intercomparison G-1 in 1977 was used as input data. The test parameter T is calculated from the equation
m 6 T = ~_J ~ (Yi/-- ~tj)2 j i 0~. '
(6)
where m is the number of peaks in each spectrum (=22), Yij is the area or location of the jth peak in t h e / t h spectrum, /3] is the weighted mean of the peak areas or locations of the jth peaks in each spectrum, 0/j is the calculated standard deviation in the programs for the peak area or peak location of the jthe peak in the ith spectrum. This parameter is chi-squared distributed with 5 X m degrees of freedom. Since the underestimation of the precision for some photopeaks may be compensated by the overestimation of others the Kolmogoroff-Smirnov test as suggested by Heydorn has also been used. In this method the vartial sums
6
(7) are computed. These sums should be random samples of the cumulative chi-squared distribution with five degrees of freedom. The test parameter d is defined by plotting the function lii Sm(T) =
N
for T < T 1 , for Ti <~ T <~ Ti+l , i = 1 ..... m - l , for T > T m
,
(8)
where m is the number of peaks in each spectrum on the same graph with the cumulative chi-squared distribution. The maximum deviation of the function Sm(T) from the ×2 (5) cumulative distribution is the test variable 3 which can then be tested. If the program passes both tests it is in statistical control and gives results with reliable standard deviations. The programs are then tested for significant deviations from the true values to ascertain accuracy. To the test of accuracy the absence of significant deviations between the true ratios of the replicate spectra to the reference spectrum as given by the IAEA are compared to the ratios of the weighted means of the peaks in the replicate spectra to each comparator peak in the reference spectrum. The test parameter X is computed using the equation
m X = ~ ~ i -/a/) 2 J
07
(9) '
where ~] is the weighted mean of the jth peak areas or loca tions of each spectrum, /a] is the true value for the ]th peak area or location calculated from the reference spectrum with a known attenuation coefficient or shift, Oj is the calculated standard deviation for the weighted mean ~j, m is the number of peaks in each spectrum (=22). The test parameter X is chi-square distributed with 22 degrees of freedom and can thus be tested against the theoretical distribution. The average precision of SAMPO and SAMPO80 can further be compared to each other by calculating the mixed term
X12 = ~ J
(ft~l)--]d-j)2 0.(2)2 ,
(10)
where the squared deviations from program 1 are
M.J. Koskelo et al. / SAMP080
97
Table 1 The results of the overall precision test for peak location calculations
Table 4 The results of the test for accuracy for peak location calculations
Program
T
Degrees of freedom
P (x 2 > T)
Program
X
Degrees of freedom
P (x 2 > T)
SAMPO SAMPO80
12300 128
110 110
<10 -7 0.11
SAMPO SAMPO80
1.06 × 1019 21.5
22 22
<10 -7 0.49
Table 2 The results of the overaU precision test for peak area calculations Program
T
Degrees of freedom
p (×2 > 7")
SAMPO SAMPO80
357 112
110 110
<10 -7 0.43
weighted with the reciprocal variance of the mean from program 2.
4.2. The results of the performance tests The values for the test parameter T [eq. (6)] were first c o m p u t e d for b o t h p r o g r a m s a n d for peak locat i o n s and areas. The results are summarized in tables
tables 1 and 2, respectively. It can clearly be seen from these results that the SAMPO80 passes both tests but the SAMPO does not. To check that the underestimation of the precision for some photopeaks is not compensated by the overestimation of the errors for others the Kolmogoroff smirnov test was applied to the SAMPO80. The graphs of the SIn(T) and the X~ (5) cumulative distribution for the peak locations and peak areas are drawn in fig. 8 and 9, respectively. The results are also summarized in table 3. The results for both cases are random samples of the x2-distribution and the SAMPO80 can thus be said to be in statistical control. The test of accuracy was then performed to verify the absence of significant deviations from the known true values. The test values for the peak locations are
Table 5 The results of the test accuracy for peak area calculation Program
X
Degrees of freedom
P (x 2 > T)
SAMPO SAMPOS0
1510 22.2
22 22
<10 -7 0.45
given for both the SAMPO and SAMPO80 in table 4 and the values for the peak areas in table 5. The SAMPO80 again passes the test on both accounts and thus gives us both accurate and reliable results. The two programs were then compared by calculating the mixed term given in eq. (10). The results of this test are summarized in tables 6 and 7. It can be seen from the tables that the two programs do give different results. Since both programs use exactly the same peak shape function and both are given the
100
5O
..J
0
Table 3 The results of the Kolmogoroff-Smimov test of cumulative distributions for the SAMPO80 Case
d
rn
P (d - 0)
Peak locations Peak areas
0.14 0.12
22 22
>0.20 >0.20
Q. ^
5
10 UNITS OF T
15
Fig. 8. The Kolmogoroff-Smirnov test curves for the peak locations obtained for the IAEA G-1 spectra 300-305 with SAMPO80.
98
M.J. Koskelo et al. / S A M P 0 8 0
100
I
a suitable constant as is done in the case of the peak location error estimate of the SAMPO80. The problem of incomplete iterations would not be removed by this method, however.
,
5. Conclusion 50
_J cD o 0~
D_
I
I
I
5
10
15
UNITS OF T
Fig. 9. The Kolmogoroff-Smirnov test curves for the peak areas obtained for the IAEA G-1 spectra 300-305 with SAMPO80. same input data one can draw the conclusion that the difference lies in the method of solving the least squares equation in the fitting process. Inspecting the input data, i.e. the peak locations and areas and their associated errors for the individual peaks it can easily be seen that the method of SAMPO has some inherent problems in the iteration process. Some peaks have been determined with excessive errors whereas most of the error estimates tend to be too small which is also indicated by the large T and X values. Because the errors are systematically too small the errors could be made larger by multiplying them with Table 6 Comparison of precision and accuracy for peak locations given by the SAMPO(1) and SAMPO80(2) Parameter
Value
Degrees of freedom
Significance
X12 X21
7.5 X 1017 1.06 X 1019
22 22
P < 10-7 P < 10-7
Table 7 Comparison of precision and accuracy for peak areas given by the SAMPO(1) and SAMPO(2) Parameter
Value
Degreesof freedom
Significance
X12 X21
69.3 78.7
22 22
P < 10-6 P < 10-7
The SAMPO80 described in this paper is an updated and improved version of the original SAMPO coded for gamma spectrum analysis. The basic philosophy of the SAMPO has been retained including the peak shape calibration procedure and the use of the SAMPO control commands. To make the code more useful to the numerous minicomputer users the original code has been made into three separate programs in SAMPO80 thus making the implementation of the code on a minicomputer possible. This process has also necessitated some changes in the code mainly to speed up the computations. Also a nuclide identification method has been added. The accuracy and precision of the new SAMPO80 has been tested using the statistical method suggested by Heydorn. The results show that the SAMPO80 gives both consistent and accurate results. The same tests were also done with the original SAMPO and the results are not quite satisfactory. This is largely explained by the use of mere fitting uncertainties obtained from the nonlinear minimization algorithm for the error calculations in the original SAMPO, rather than the proper statistical estimates. The overall performance of SAMPO80 as compared to the original SAMPO is improved in peak search, significantly accelerated in peak fitting with comparable results except for unresolved multiplets, much improved in the error estimates, and extended to nuclide identification. The analysis can be completed in a minicomputer and the fitting times are reduced by an order of magnitude. Thus for most applications, including activation analysis and radioactivity studies, but still excluding stringent nuclear spectroscopy applications, the new SAMPO80 can be recommended over the original SAMPO. We are grateful for the help of S. Toivonen and J. Palom~'ki for doing some of the recoding o f the original code for our minicomputer.
M.J. Koskelo et al. / SAMP080 References [ 1 ] J.T. Routti, University of California, Lawrence Berkeley Laboratory Report, UCRL-19452 (1969). [2] J.T. Routti and S.G. Prussin, Nucl. Instr. and Meth. 72 (1969) 125. [3] R.M. Parr, H. Houtermans and K. Schaerf, Proc. ANS Topical Conf. on Computers in activation analysis and gamma-ray spectroscopy, Mayaquez, Puerto Rico (1978) U.S. Dept. of Energy Symposium Series no 49 (CONF-780421) (1979) p. 544. [4] G.C. Christensen, M.J. Koskelo and J.T. Routti, CERN, Health and Safety Division Report HS-RP/015/Rev. (1977 revised 1978). [5] M.J. Koskelo, P.A. Aarnio and J.T. Routti, Helsinki University of Technology Report TKK-F-A427 (1980); to be published in Comp. Phys. Comm. [6] W.C. Davidon, Report ANL-5990 (1959). [7] P.A. Aarnio, M.J. Koskelo and P. Zombori, Nucl. Instr. and Meth. 184 (1981) 487.
99
[8] M.A. Mariscotti, Nucl. Instr. and Meth. 50 (1967) 309. [9] M. Giannini, P.R. Oliva and M.C. Ramorino, Comitato Nazionale Energia Nucleare Report CNEN-RT/FI (72) 14 (1972). [10] W.C. Schick, Jr., PEAKFIND, a Fortran program for locating and fitting peaks in semiconductor detector spectra, Iowa State University, Ames Laboratory Report IS-3636 (1975). [11] M.J. Koskelo, J.T. R0utti and S. Toivonen, Annual Meeting of the American Nuclear Society, Atlanta, Georgia, USA, 1979 (Helsinki University of Technology, Department of Technical Physics Report TKK-FA390) (1979). [12] R. Gunnink and J.B. Niday, Computerized quantitative analysis by gamma-ray spectrometry, vol. 1, University of California, Lawrence Livermore Laboratory Report, UCRL-51061 (1972). [13] K. Heydorn, in ref. 3 (1979) p. 85.