Nolin: Nonlinear analysis of complex alpha spectra

Nuclear Instruments and Methods 185 (1981) 261-269 North-Holland Publishing Company

261

NOLIN: NONLINEAR ANALYSIS OF COMPLEX ALPHA SPECTRA E. GARCIA-TORANO and M.L. ACENA

Divisi6n de lnvestigaci6n Bdsicay Aplicada, Junta de EnergiaNuclear,A venida Complutense 22, Madrid (3), Spain Received 8 July 1980 and in revised form 25 November 1980

An equation has been deduced which appropriately represents the form of an experimental alpha peak obtained with Si surface barrier detectors. On this basis a nonlinear alpha spectra analysis method is proposed. The results are checked by applying the method to synthetic and real spectra.

1. Introduction Despite the great progress in alpha spectrometry with semiconductor detectors, the resolution obtained is insufficient to separate very close lines, and so it is highly important to develop spectral analysis methods. These methods have not yet reached the level obtained in gamma spectrometry due to the great complexity of an experimental alpha line, which depends on factors such as type of detector, geometry of the measurement, sample, electronic equipment, e,tc. Recently, several methods of alpha spectra analysis have been presented together with computer programs [ 1 - 3 ] . Many of these are based on modified gamma spectra analysis methods thus enabling strong asymmetry effects to be taken into account. Generally, peak shapes which represent very closely the real shape of an experimental alpha line are not used. NOLIN is an alpha spectra analysis program based on a nonlinear least-squares fitting. For its application it has been necessary to draw up a new peak shape which can be sufficiently adapted to real spectra so as to allow analysis. Firstly, this function will be studied, then the method of analysis and the results obtained.

2. Description of the line shape

To formulate an equation which closely describes the alpha line shape it is necessary to bear in mind the following experimental facts: 0 029-554X/81/0000 -0000/$ 02.50 © No rth-H olland

1) The region between the maxamum" intensity and the end of the line is Gaussian. 2) The line is asymmetrical, the high and low energy regions being very different. This asymmetry is a characteristic of the spectrum, and depends on such diverse:factors as the detector, the sample used, geometry of measurement etc. 3) There is a low energy tail which extends to the origin and may be expressed as a percentage of the peak area, according to the type of detector used [4]. The decrease towards this low energy tail is asymptotic. In order to study the alpha line structure, numerous spectra from monoenergetic emitters have been obtained. The region between the maximum and the end of the line has been adjusted to a Gaussian and this was substracted from the spectrum. It was observed that the residual spectrum has a shape similar to the initial one; a Gaussian with identical fwhm and its maximum displaced by a distance equal on average to the fwhm. Consequently, the proposed curve consists (fig. 1) of two regions: 1) The high energy region, comprising the sum of the two Gaussians with heights A and B, with identical fwhm, the smaller of the two being responsible for the asymmetrical effects, the equation being:

y(x)=Aexp(

(x-x°)2~

_ (x - Xo + po) 2) + B exp

202

,

E. Garcia-Torago, M.L. Acega / NOLIN

262

0z HYPERBOLIC RANGE

-~ ~

GAUSSIAN RANGE

8

LU m -J

I

I

I

I

~

t ENERGY

Fig. 1. Proposed line shape.

2.1. Determination o f the high energy region parameters

where: Xo

= abscissa of the maximum,

Y(x) = value of the function at the abscissa point x, p

po

= 2x/2 in 2. = fwhm of the Gaussians.

2) The low energy region, which follows a hyperbolic equation;

Y(x) -

K

(a - x)' t- C,

where K, a and t are characteristic parameters and C the low energy tail, constant to the origin. The abscissa of the joining point of the two curves is that o f the maximum of the least energetic Gaussian. The general expression of the curve is, therefore:

Y ( x ) = l +sgn(x-x°+Pa)IAexp(-(x-x°)2~2

-2~f

Fig. 2 shows a region of a typical monoenergetic alpha line spectrum. To characterize the asymmetry of any real peak we define a parameter called "asymmetry factor", the value of which is given by:

f = Ri/R d . Evidently, the minimum value of this parameter will be f = 1, corresponding to a perfect Gaussian. The upper practical limit of this factor is around 1.4 or 1.5. With the given definition of the asymmetry factor and the condition that the maximum height of the

]

PO-210PARTIALS P E C ~

+ B exp (- (x - Xo + po) 2 z

f

Ri

/

r

\

8

where sgn(x) is the function sign of x, defined as follows: sgn(x) =

i ifx>O' if x = 0 , -- i f x < 0 .

ENERGY Fig. 2. Asymmetry factor definition.

E. Garcia-Torafio,M.L. Ace~a / NOLIN line is h, the following values for A and B may be obtained (see appendix 1)

A = h I1

½exp(p2f2/8)1 21 i

'

exp(gp f )

B

L p ( - g p )[

p(gp f ) -

"

2.2. Determination of the low energy region pararn,. eters

Due to the great diversity of shapes which the region of fall towards the tail may take on it is necessary to introduce a new parameter q called the "peak/ tail ratio". Its value is given by the quotient between the maximum intensity value of the line (h) and the value of the curve at a point separated from the maximum by a distance equal to n times the resolution, n being a variable number. In appendix 2 the conditions of joining can be seen (equality of function and first derivative) and the development which leads to the fitting parameters.

2.3. Experimental part

263

the nuclides on discs one inch in diameter, mirror polished, made of platinum, silver and stainless steel. The z l°Po is deposited on a silver disc by autoelectrolysis at 80 C [5]. The isotopes 216p0 and 212p0 are deposited as daughters of a 232U hydrochloric solution practically in equilibrium with its decay products, by the Mitchell method [6] of electrolysis on a platinum disc, with a giratory anode also of platinum, in NH4C1 electrolyte. These same isotopes are deposited onto stainless steel, by isotopic interchange in a solution at pH = 7 with a sample of natural uranium, previously prepared by vacuum sublimation with electronic bombarding of a uranium metal target. Three other sources have been prepared by the Mitchell method, two of 239pu + 24°pu and another of 241Am on a platinum support. Finally by electrostatic deposition on platinum a sample of 238pu was obtained. Spectra were obtained for all the sanples under differing conditions of geometry, using a spectrometer composed of an ORTEC amphfication chain and a multi-channel SEIN pulse analyzer. It has been observed that a good fit is obtained using the deduced equation. In figs. 3 and 4 typical spectra obtained with these samples are represented, together with the fitted line.

As the sample holder has an influence on the line shape, the sources have been prepared by depositing

100

o SPECTRUM

j~

,FITTED SPECTRUM

#

ASYMMETRYFACTOR= t.21 PEAK/TAIL RATIO= 7.35{AT4 EW.H.M.) EW.H.M.~8.4 _ J

ff

r

~

10

b o°c~o o

I

I

I

I

I

I

I

I

I

I

CHANNELNUMBER Fig. 3. Spectrum of a 21°po sample. Silver support. Detector surface 50 mm 2. The fwhm is 2R d (see fig. 1).

264

E. Garcia-Torago, M.L. Acefia / NOLIN

100 0 --

SPECTRUM FITTED SPECTRUM ASYMMETRY FACTOR = 1,1g PEAK / TAIL RATIO= 250 (AT E W.H.M. = 11.5

4

E W,H,M.)

..~ 0Y ,.~

10

z l--

z

8

E

o DO ~ o f I

O I

I

O I

I

I

l

I

I

I

I

1

I

I

I

[

I

I

I

CHANNEL NUMBER

Fig. 4. Spectrum of a 216po sample. Stainless steel support. Detector surface 25 mm2.

3. Description of the method of analysis The problem of fitting the spectrum can be reduced to a study of the minimization of the expression:

the initial Xi are estimated by calibration of the spectrum. With these values a set of R i are calculated which in turn are used to recalc,ulate the Xi, thus entering an iterative process. 3.1. Minimization with respect to R i

J =~

Wx [Yx - R l f ( x , X1) . . . . .

R n f ( X , Xn)] 2 ,

x

(1)

The minimization of eq. (1) with respect to the intensities of the lines can be expressed as

in which:

aJ

Yx Ri Xi f ( x , Xi)

aRi

w~ H

= content o f channel x, = intensity of peak i, = position of peak i, = value at point x of proposed equation describing alpha line of unit height and maximum at x = Xi, = statistical weight, with value equal to 1/Yx, = number of alpha lines analyzed.

This is carried out for all points of the spectrum analyzed. For complete adjustment of the spectrum it is necessary to know the R i (intensities) and the X i (positions). Since the minimization of the expression (1) with respect to the 2n unknowns simultaneously would lead to systems with a large number of equations, the process is done in two parts. Firstly,

-0,

i=l,2...n,

and leads to a system of equations, linear in Ri, which written in matrix form is expressed as:

/.

"N

1

The solution gives values of R i for the estimated values of Xi.

E. Garcia-Torafio, M.L. Acefia / NOLIN

3.2. Minimization with respect to Xi Once the intensities are known, the positions may be obtained by solving the system resulting from minimizing the expression (1) with respect to Xi OJ

-0,

i=l,2...n.

bXi Thus, the following system is obtained, nonlinear in Xi, with equations:

I

w~ r~-~k nkf(x, Xk) R~ aX~ = ° ' i, k = 1 .... ,n . This system may be solved by the Newton method [7]. For this it will be necessary to know the Wronsklan of the transformation, the general element of which may be writen as:

wq = ~ W~R~R~af(x, x~) af(x, xi) aXi

aXj

265

In appendix 3 the procedure for obtaining the nonlinear system and the Wronskian is given.

3. 3. Description of the NOLIN program For carrying out the calculations, a program has been written in FORTRAN language to be used in two different computers: a 28k PDP 11/20 and a UNIVAC 1110. Figs. 5 and 6 illustrate the fittings obtained for two spectra, obtained with a BENSON 1100 plotter connected "on-line" with the PDP computer. NOLIN consists of a main program and 12 subroutines and occupies about 600 FORTRAN statements. In order to accelerate the convergence of the method, a simple algorithm has been used which extrapolates from the values obtained in five ,iterations of both systems (linear and nonlinear). It is considered that convergence has been attained when two groups of positions differ by less than 1 X 10 -3 in two successive iterations. If necessary, the program takes over the task of subtracting the background coming from higher energy peaks.

--=

FITTED

SPECTRUM

CHANNEL N U M B E R

Fig. 5. Fitting results of a 238pu spectrum. Circles zepresents the measured spectrum, thin lines the components and thick line their sum.

266

E. Gareia-Torafio,M.L. Acega / NOLIN

- -

=

FITTED

SPECTRUM

I

i

r

CHANNEL NUMBER

Fig. 6. Fitting results of a 239pu + 24°pu spectrum. Circles represents the measured spectrum, thin lines the components and thick line their sum.

4. Application and results In order to check the efficiency of the method, a synthetic spectrum has been used, built in accordance with the proposed line shape. Despite the initial attribution of erroneous positions to the maxima of the component peaks, the algorithm has been able to find the correct values after a reasonable number of iterations. The results obtained are shown in table 1.

The application of the method to real spectra such as represented in figs. 5 and 6 presents several problems. Due to statistical fluctuations it is difficult to determine the channel corresponding to the maxim u m energy. As the spectrum obtained in a multichannel analyzer is not a continuous curve but a histogram, at times the pulses corresponding to the maximum energy do not accumulate exclusively in one channel but are spread between two of them. All

Table 1 Results of the application of the method to synthetic spectra. Fwhm (chan.)

8.4

Theoretical

Initial

Final

Iterations

Channel

Intensity

Channel

Intensity

Channel

Intensity

10 20 28 35

0.2 0.8 0.5 1

10.5 20.5 27.5 35.5

0.236 0.796 0.546 1

10.000 20.000 28.000 35.000

0.200 0.800 0.500 1

12 21 29 37

0.293 1 0.970 0.813

9.994 20.000 28.000 35.000

0.200 0.800 0.500 1

19

E. Garcia-Tora~o, M.L. Ace~a /NOLIN these causes lead to a masking of the real maxima, and may lead to considerable errors when fixing the resolution parameters and "asymmetry factor", wich must be considered. An examination of the spectra obtained leads to the conclusion that although the "asymmetry factors" range from 1.05 to 1.4, normally, the "peak/ tail ratios" depend strongly on the sample and on the geometry of the measurement, the tail being much higher in the silver and platinum backed samples than in those of stainless steel, in which the low energy tail falls far more rapidly. The information given in figs. 3 and 4 enable these differences to be appreciated. Whereas the attribution of values to the parameters of "asymmetry factor", "peak/tail ratio" and fwhm is easy in a simple spectrum, with two or more non-overlapping lines, in the case of a multiplet problems arise which impede this from being done immediately. It is possible to estimate, knowing the characteristics of the equipment and the sample, the "asymmetry factor" and fwhm, or measure them in a region of the spectrum in which the approximate positions of the existing peaks are known. As regards the "peak/tail ratio", one possible method is to spread

Table 2 Results of the application of the method to real spectra. Spectrum

Atoms of 24°pu per 100 of 239 Pu obtained

1 2 3 4

8.470 8.355 8.483 8.290

Fwhm (keV)

real 8.65

18 17.5 17.2 18.

267

out the tail at a determined point in a way directly proportional to the approximate height of the peaks which contribute to it and inversely to their distance to the point of reference. In this way the spectra of Pu have been analyzed. The method has been applied to six different spectra, four of them coming from a sample of 2agPu + a4°pu and the other two from 238pu. The results are shown in table 2. For the case of the plutonium mixture, the number of atoms of 24°pu per 100 of 239pu has been calculated and compared with the known value.

5. Conclusions From an examination of the figures and information obtained, it can be deduced that the proposed line shape closely represents a great variety of experimetal alpha lines and mixtures of these. As regards the method of fitting, the test carried out with synthetic and real spectra have led to satisfactory results. For the fitting of spectra with a large number of lines, it is advisable to separate groups of lines which can be studied as a single problem since the number of matrix elements which have to be studied varies as the square of the number of peaks which intervene in the fitting. This is particularly important for reducing the calculation time and facilitating the convergence in the fitting. The time used to carry out the fitting depends largely on the number of lines fitted. A typical time for a 55 channel spectrum with 5 lines, as shown in fig. 6, is between 30 and 9 0 min for a PDP 11/20 computer, depending mainly on the accuracy of the assignment of initial positions to the peaks.

average 8.40 Spectrum

1

2

Branching ratios of a-rays emitted from 238 Pu obtained

real

72 28.1 0.081 72 27.9 0.13

72 28 0.1

average 72 28 0.105

Fwhm (keV)

Appendix I

Determination of the high energy region parameters 29 28

Consider a real peak, characterized by the corresponding channel at its maximum Xo (fig. 7), the maximum height h and the asymmetry factor f. In a real spectrum, the height of the least energetic Gaussian, responsible for the asymmetry effect, will reach a small value in comparison to that of the main Gaussian, and, given their separation, it may be estimated that the fwhm of the region of the line to the right of the maximum will coincide with that of the Gaussian

268

E. Garcia-Torafio, M.L. Aee~a / NOLIN

Appendix 2

//I

ASYMMETRY FACTOR = 1.2 E W.H.M. = 16.5 A = 0.9867

/~

l/

//

0 z pZ 3 0 U

]

Determination o f the low energy region parameters As seen already, the equation for the low energy region is hyperbolic, adopting the form: Y(x)

Ri

i,i

Lll a-

15

30

45

60

CHANNEL NUMBER

Fig. 7. Composition of the alpha line in its region of highest energy.

of height A. The following approximations may then be made:

R d ~ 1.180,

K

(a - x ) t I- C .

So that there are no discontinuities between the two regions of the curve, it is necessary that at the joining point of the two curves, the values of the high and low energy functions on the one hand and their first derivatives on the other, should be equalized. Let (xl, Yl) be the point of joining of the two curves, the position of which is discussed at the end of this section and (x2,y2) the point with respect to which the "peak/tail ratio" is defined (x2 = xo - npa). Then, the condition of equality of the functions gives us the equation: Z = A exp ( - (x, -- Xo)2~ \ 202 ]

+ B exp (-

(xl - Xo + po)2~ ~2

Ri =fRcl "" 1.18fa.

The conditions which should be imposed on the proposed function are as follows: a) The value of the function at x = Xo must be h

K

"J -- C - (a - XI) t '

(2)

and the equality of the derivatives gives us: Z' = - A (xl - Xo) exp 02

(xl - xo) 2 ) ~

[ (P0)2~

h =A +B expk--~-2

] •

- B (Xl - xo2 +

b) At the point corresponding to the abscissa x =Xo - 1.18fa the value of the function must be ½h by definition of the asymmetry factor we have:

=t

e x p ( - (xa - x° +

K (a -- X1) t+l"

If the following parameter is introduced:

~=Aexp(-(l~2t+Bexp(-(lP~+zPa)2 /

)

a =Z/Z', we obtain:

Starting from these equations, expressions for the values A and B are obtained for the equation in the high energ)) region. A=h[1-

½exp(p2f2/8)- I1 1 2 exp(~p f)

pl.-~p ) l

pl.~t~ s ) - 1

a--xl t = - - , c~

(3)

an equation which connects the parameters a and t, asymptote and exponent of the hyperbolic curve respectively. By definition of the peak/tail ratio, we have:

q =h/

K (a - x2) t

+C

E. Garcia-Tora~o, M.L. Acega /NOLIN the minimum condition will be:

and operating: 1

1

q

h

1

-

269

(4)

Co '

where:

aJ 2 ~ W~fyx_~k Rkf(x,Xk)1

~bi =--- ~ i =

× (-R~~I(x,~2:]0 . X , )_-]

co = (a - x 2 ) t / I ( .

Eliminating K between the equations (2) and (4), we obtain:

That is,

¢i = - ~ x

(a - x t ) ln(a - x2t = ~ ln~--h-[Zco~) . \a - X t l This equation may be solved by the Newton method, giving the following recursion equation:

= 0 Wx I g x - ~ R k ~ .fx( , Xk) I Ri -~f(x'-Xi)~Xi .

As is known, the solution to this system can be obtained by the Newton method and is given by: ¢(X p)+ W(Xp) E p = 0 ,

in(an-x2__~] an-xt o~in (_Z~_%) \an - x l / ln(an--X2~+(an--Xm]_l \an--Xl] \an--X2/

an+l=an

(s)

.

k

X p+t = X p + E p ,

where X p is the matrix of solutions corresponding to the pth iteration and W is the Wronskian of the transformation, defined as:

For the application of this equation it is necessary to consider the point (Xl, Yl) of joining of the two curves. Experience shows that the most recomendable point is that which corresponds to the abscissa

wiJ=axj

xl = xo - 1.1 8for,

wq =~ ~ ~ WxRjf(x, xj) R,. ~f(x,~Xi) ,

Any term of W will be:

that is the point at which the curve takes the value

lh. An initial value for the asymptote which has been found as valid is:

al = Xo -- ½PO . From eq. (2), we obtain:

K = Z(a - x l ) t .

since the remaining terms have zero derivative with regard to Xj. Thus we obtain:

Wi] = ~x WxRiRj ~f(~xiXi) ~f(x, X])

•

(6)

The equations (3), (5) and (6) enable the parameters for the complete fitting of the curve to be determined. Due to the high values adopt, it is advisable to work with the curve in a logarithmic way.

axj

'

which is the general expression of the Wronskian.

References Appendix

3

Construction of the nonlinear system and resolution by the Newton method Consider the expression (1):

J = ~x Wx Yx -

k R k f ( x , Xk

, k = 1 ..... n ,

[ 1] M.P. Trivedi, Report NYO-844-81 (1969). [2] H. Baba, Nucl. Instr. and Meth. 148 (1978) 173. [3] W. Watzig and W. Westmeier, Nucl. Instr. and Meth. 153 (1978) 517. [4] M.L. Ace~a and E. Garcia-Torago, Report J.E.N. 442 (1979) (in Spanish). [5] M. Frisch and I. Feldman, 1. UR-426 (1956). [6] R.V. Mitchell, Anal. Chem. 32 (1960) 236. [7] M.A. Wolfe, Numerical methods for unconstrained optimization (Van Nostrand-Reinold, London, i978).