Automatic analysis of γ-ray spectra including multiplets

246

Nuclear Instruments and Methods in Physics Research A258 (1987) 246-249 North-Holland, Amsterdam

A U T O M A T I C A N A L Y S I S O F T-RAY S P E C T R A I N C L U D I N G M U L T I P L E T S H. M A C H N E R

Institut fiir Kernphysik, Kernforschungsanlage, Jiilich, D-5170 Jiilich, FRG Received 2 October 1986 and in revised form 16 February 1987

A method is developed which searches for photopeaks and identifies multiplets. This identification is based on three different procedures which are partly complementary. The method is applied to complex spectra including weak peaks on a large and strongly fluctuating background. 1. Introduction

The advent of high resolution -/-ray detectors like Ge diodes influenced the development of numerous computer codes which allow for the determination of the energy and area of photopeaks. A rather large number of such codes is cited in ref. [1]. The usual procedure is to fit a background function and a model function to measured pulse height spectra. Because the background may be rather complex due to processes like Compton edges and because the peak shape may be distorted from a pure Gaussian form depending on the actual detector, mathematical formulae were chosen for the description of both with different degrees of sophistication. While not too complex spectra with separated peaks are well suited for automatic peak search and fitting algorithms, complex spectra are not. Here, very often different peaks overlap each other. It has, therefore, become customary to analyze such spectra interactively, i.e. background, -/-lines and the corresponding multiplicity are selected by the experimenter on a monitor. This is a rather time consuming procedure and, therefore, automatic algorithms which allow also for fitting of multiplets are called for. In this contribution we will present a peak search algorithm which allows for the detection of multiplets in cases where there is still some separation. This is most often the case. Stronger overlapping peaks are fitted as single peaks. Their signature is then a broader peak width than those for real single peaks. Then start parameters are derived for the following minimization of the square sum of the deviations between data and the model function. 2. Peak search algorithm

As already stated above, some photopeaks may be superimposed on a fluctuating background [1]. A possi0168-9002/87/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

ble peak search algorithm may, therefore, be based on the inspection of the frequencies in the Fourier transformed spectrum. Because the transformation processes need rather long computing times in cases of larger channel numbers, here another way was chosen. To get rid of the strongly fluctuating background some smoothing procedure may be appropriate. For this purpose the convolution method of Savitzky and Golay [2] was employed. In this method a data point Y, in the channel i is replaced by a number Y~ obtained by fitting a polynomial to an interval around channel i. Fitting a polynomial of order k to 2k + 1 data points is, for equidistant data, equivalent to multiplying the data with constant numbers aj. The (smoothed) fit value is k

~=

E

ajYi+j"

(1)

j~-k

The constants fulfil the relation aj = a_j. If we assume the smoothed background to be hnear, the second derivative will show almost no background. We will further make use of the fact that the photopeaks are approximately Gaussian. The second derivative of a Gaussian is a peak with a width approximately a factor of 2 smaller that that of the Gaussian itself. The smoothed second derivate can be calculated again by the convolution method described above. If one approximates the derivative of a curve in a point by the difference of the two neighbouring data points divided by their distance, one again obtains constant numbers cj, so that the smoothed second derivative Y," in channel i is derived from the counting rate Y; by k

r,"= ~ cjr,+,.

(2)

j_-~-

The counting rate numbers Y,. become smaller for increasing channel number i due to the smaller detection efficiency of the counter for higher T-energies. To get a nearly energy independent quantity we define a filter

H. Machner / Automatic analysis of y -ray spectra 100

[

'

I

I

b

247

I

i

-6E+04 Q) >

60-

°_

-5E+04

"6 > k-

~0 "I]

2 0 - ~

_~-4E+04

cq "I3

-20

N

._

-3E+04

-6

c ct3 cO

(0 0

E k-

O t-

~

160--

-100 400

4½0

Ao

&o

480

_-2E+04

500

channel hr. Fig. 1. Part of a y-ray spectrum (lower curve and right abscissa) is shown together with the filter function eqs. (3) and (4) (upper curve and left abscissa). Also shown is a broken line indicating a threshold at - 5. function

F, = ~ " / ~ r , " ,

(3)

with

AY,":

y" cZY~+s.

(4)

j=-k

This filter function F i is shown in fig. 1 together with the original spectrum. A line is detected at position i if the ~ fulfil the following three conditions: F i ~< - K,

(5a)

F,+I >/F,, F,_ 1 >/F,,

(5b)

and

F,_, + F,+, ~
(5c)

with K a positive number chosen to discriminate between photopeaks and the fluctuating background. Typical values employed are between 2 and 5. The second criterion is due to an assumed maximum in channel i of the line while the third criterion was empirically found by N y m a n [3]. Suppose a peak was found at channel number i. Now, the interval has to be fixed to which background and model function should be fitted. Therefore, we look for the 2nd channelil fulfilling Ell ~ 0 and F i l + l ~ 0.

(6)

We assume il to be the left boundary. N o w a similar procedure is applied to find the right fitting interval boundary at i2: Fi2-1 >/0 and F/2 ~< 0.

(7)

It may, however, happen that in the region from i to i2

there is one more line. If this is the case the multiplicity is increased by 1 and we proceed to find a new interval boundary. The so defined loop is repeated until no further line is found in the interval between the last found line and the last defined i2 value. Furthermore, we have limited the size of the interval to have a straight line as an acceptable approximation for the background. The case of a doublet is shown in fig. 2.

3. Fit to photopeaks There are two more ways to identify peaks. These methods can be employed only after fitting a model function M to the data. For the sake of simplicity only pure Gaussians and a linear background are considered here:

M(i, a))=ao+ali+ EaK1 exp[aK2(aK3--i)2], Ki

(8) with parameters a which are varied so as to minimize the quantity i2

Q= Y~ Gi[Y,- M(i, aj)] 2.

(9)

i--il

Here, G i denotes the statistical weight of the channel content ~ . At this point we would like to stress the fact that the fit is performed to the original data and not to smoothed data. For a discussion of the influence of data smoothing to peak areas the reader is referred to ref. [4]. It was found in the analysis that in the case of multiplets the fit results in very large values for aK2 corre-

H. Machner / Automatic analysis of 7 -ray spectra

248

I

1

1E+06

I

I

20-

i2' i2.

//~\

>

0

^v - ^~----

>

1

I -20-

-40-

I I

152Eu

- 6 0 ~ 8 t"-

I I

1

I

-80

l

I

I

I1) N

C

0 U

I

32E40 32J60 ' 32180 chonnel nr.

32J20

-1E+04

I t

I

-100 , 3200

-1E+05 Cc

3300

Fig. 2. Same as fig. 1 except for the threshold. The left boundary il is given by eq. (5). The first channel fulfilling eq. (6) is indicated as a dashed line without a label and the second channel, being a tentative right boundary by i2'. Because between the peak maximum and i2' there is a second peak the next channel fulfilling eq. (6) is denoted by i2. The linear background being used as a start value in the fitting procedure is also shown, as a solid line in the spectrum.

sponding to small widths for the small peaks and vice versa. This can be overcome by assuming simply the same value a,r2 for all lines in one multiplet. Another possibility to detect overlapping peaks which are not identified as such by the peak search algorithm is based on the fwhm f given by fK=

2,/T£2 V aKz"

(10)

If this fitted width fK of a single peak is larger than the one in the calibration spectrum f ( E ) , i.e.

fK>~f(E)

and

fK-f(E)

>~AfK,

a doublet is assumed to exist. The uncertainty in f is denoted by zlf. Starting values for the fit are than taken from counting rates at positions a°3 + f ( e ) / 2 with a°3 being the position of the previously assumed singlet. The third possibility to find peaks is based on the fit itself. If the fit results in a Q-value (eq. (9)) with

O >~x 2 ( i 2 - il - r, 1%),

c

-

oo..o.0o

. . . . . . . . . . . 2A-

!150

,

,

11170

I ,

11~80

•

data

11~0 ' 12~00 ' 12~10 ' 12~20 chonner

,

-~

1230

hr.

Fig. 3. Part of a 7-ray spectrum. The data are fitted by a parabolic background and six Gaussians (solid curve). The fitted Gaussians on top of the background (dashed curve) are separately shown (short dashed curves).

(12)

where r denotes the n u m b e r of fit parameters, the quantity

D, = Y i - M ( i , aj)

Zo

(11)

(13)

is inspected. If D~ >/0 for a n u m b e r of neighbouring channels and this n u m b e r exceeds the value fx, then another peak is assumed to exist in the middle of the corresponding interval. The peak is accepted as such if the fit results in a Q-value smaller than the previous one. Finally, it is checked whether a parabolic background improves the fit. A n example of the method described above is shown in fig. 3. The four strong lines are detected via the peak search algorithm discussed in the previous section. The two weak lines are identified due to positive differences (eq. (13)). It should be mentioned that from single lines we know that the Gaussian is a good approximation to the photopeaks and that the two weak peaks are not tails from larger peaks. This can also be seen for example from a lower energy part of the same spectrum shown in fig. 4. In this case the peak

H. Machner / Automatic analysis of 7 -ray spectra

o

¢>

i

~,o/~

p'

,,

•

i

':<

i ,-,,_~s

270

.

\

280

290 channel

bOO

hr.

Fig. 4. Similar to fig. 3 but lower energy part of the same spectrum.

search algorithm has detected only one peak in channel 280. The fwhm of the fit was much larger than the one of the calibration spectrum and therefore two peaks were fitted. The third peak was detected by the last criterion because of the positive difference (see eq. (13)).

4. Discussion

We have presented a method for an automatic analysis of -y-ray spectra taken with high resolution Ge detectors. This method allows for the identification of multiplets. Three different criteria are applied to identify multiplets as such: - peak search algorithm, - fwhm, - fit procedure. The first criterion identifies only those peaks which are only partly overlapping or, in other words, still show some separation. The next criterion is effective only for peaks originally detected as singlets. It is thus complementary to the first one.

249

To get an idea to what extent the present results are comparable with those derived by interactive methods, we have applied an interactive version of the code SAMPO [5]. In this code, first a background is defined left and right to the peak. Then a model function Gaussian plus an optional exponential tail at the low energy side - is fitted to the differences between counts and background. The background is approximated by a polynomial. It is found that the peak positions are given at identical positions by the two codes. The error in the position is slightly larger in the present method but still sufficiently small. Typical values are 0.20 keV for weak peaks and 0.02 keV for strong peaks. The error in the second quantity of interest, namely the peak area, depends also on the peak size itself. The results of the two codes agree within the error bars. However, there is a tendency of the present code to give larger areas for small peaks, because the simultaneous fitting of peak and background gives smaller background values. This seems to be the better approximation, since in the present method there is no need for large intervals left and right to the peak region, which possibly contain very small peaks, to define a background. In summary we have developed a method for leastsquares fits to photopeaks. This method allows for the automatic identification of multiplets. The results are comparable with interactive methods.

References

[1] W. Westmeier, Nucl. Instr. and Meth. 180 (1981) 205. [2] A. Sawitzky and MJ.E. Golay, Anal. Chem. 36 (1964) 1627. [3] B. Nyman, Nucl. Instr. and Meth. 108 (1973) 237. [4] H. Tominaga, M. Dojyo and M. Tanaka, Nucl. Instr. and Meth. 98 (1972) 69. [5] J.T. Routti and S.G. Prussin, Nucl. Instr. and Meth. 72 (1969) 125; R.M. Lieder, private communication.