Study of scintillation spectrometry: Unfolding methods

NUCLEAR INSTRUMENTS AND METHODS

54 (1967) TOj-Il5; ©

NORTH-HOLLAND PUBLISHING CO.

STUDY OF SCINTILLATION SPECTROMETRY: UNFOLDING METHODS Y. S. SU· Institute of Nuclear Science, National Tsing Hua University, Hsinchu, Taiwan

Received 13 January 1967 and in revised form 28 March 1967 Analysis and comparison of the inverse matrix method, Scofield's iteration method and the least-squares method are made by utilizing an ideal "measurement" of an assumed spectrum. It is found that Scofield's method is, in general, a fairly good method while the simple least-squares method can only be used for a small number of energy intervals. Some new approaches to use integral quantities instead of differential quantities in the response

relation and to assume a polynomial approximation for the integral responses and the spectrum are proposed for further study, in order to get better interpolation and to use polyenergetic sources as well as monoenergetic sources for constructing the response relation. A modified technique for Scofield's iteration is also presented for faster convergence.

1. Introduction In measuring a spectrum by scintillation spectrometry, the pulse-height distribution y(V) and the true spectrum x(E) are related to each other through

approximate solutions which satisfy the physicallimitation resulting in an even better answer, have to be used. One of the difficulties in comparing the different unfolding methods and therefore to improve unfolding techniques is due to the lack of information of the true spectrum as standard so that the calculated spectra by different methods can be compared. Our other paper4) proposed a method of ideal "measurement" so that y(V) can be obtained by integrating eq. (1) with the assumed A (E, V) and x(E). With the right point approximation, Y, A and the true X can be determined. By utilizing the Y and A with different unfolding methods, the calculated X's can be obtained and then be compared with the true X. IBM 1620 has been used to calculate X for the different dimensions of A, n = 11, 21,31,41 and m = 10,20,30 with the condition m< n. For application of scintillation spectrometry to handle a bunch of data in radiation shielding, the computing time has to be taken into consideration. Modification of Scofield's iteration method for saving computer time is also presented.

y(V) =

f:

ma •

A(E,V)x(E)d£,

(1)

where A (E, V) is the response function of the spectrometer produced by a single photon of energy E. The usual approach is to approximate eq. (1) by a matrix equation 1)

Y=AX,

(2)

where Y and X are distribution vectors and A is the response matrix which can be obtained by interpolating the experimentally obtained responses of the spectrometer resulting by several monoenergetic radiation sources. The purpose of this paper is to analyze and to compare some current methods which have been used to unfold the pulse-height distribution Yor to solve spectrum X basically from eq. (2), and to propose some new approaches to reduce the errors of obtaining a response relation. 2. Comparison of methods There have been published many kinds of unfolding methods. Among them, the direct inverse matrix method 1), Scofield's interation method 2) and the simple least-squares method 3 ) are most frequently used. . Mathematically, to solve X, from a known Y by Inversing A e.i.,

(3) is the most simple and exact way, but, since A cannot be accurately determined in an actual case, the exact solution becomes physically meaningless, and some more complex methods, e.g. an iteration method, with

2.1. INVERSE MATRIX METHOD1) If A is a square matrix, the spectrum can be directly solved by inversing A according to eq. (3). As shown in fig. 1, the calculated spectra obtained from eq. (3) with Yprovided by an ideal measurement always oscillate along the assumed or true spectrum. For small matrix (10 x 10), no negative element appears, but for the other two cases, oscillation is very violent so that negative values are developed. It is therefore concluded that the inverse matrix method can not be used for spectrum with peak even if no negative value is developed. Since no experimental and statistical errors are involved in an ideal measurement, oscillation must be due • Present address of the author: Department of Chemical Engineering, University of Salford, Salford, Lancashire, England.

109

110

Y. S. SU

60

v

50

v

40

-

o n .. 11

m.. IO

v n=21 IIn=31

m=20 m=30

ideal "measurement". Physically, this method simply means to choose m good data from n data as basis for calculation. By this method, eq. (3) is replaced by

TRUE SPECTRUM

30

v

20 10

where Y' = AT Y, A' = AT A and AT is the transpose of A. The calculated results are shown in fig. 2. Some points missing in the figure are out of scale. If we make a close examination of the missing points, it is not difficult to conclude that only for the cases 20 x 10, 30 x 10 and 40 x 10, all points are complete. In other words, the least-squares method is only good for a small number of energy intervals. Furthermore, for these cases, the simple least-squares method gives better results than the results of Scofield's method as shown mfig. 3.

v

v

V

7.0

V

II

II

II

0.8 II

in 0.6 z ~

0.4

1&.1

0.2

~

>

iii

0

II 0

1Ivfl'1

0

-J

V

II

1.0

1&.1

II

II

II

3.0

>~

II

II

~O

v

II

II

~ -2.0

II

-4.0 -6.0

V

-8.0

II

/I

II II

2.3. SCOFIELO'S ITERATION METH002) In order to obtain a reasonable spectrum with physical meaning, Scofield proposed an iterating scheme to modify the conventional incremental correction by a factor correction. Actually the modificat1on imposes a limitation or condition

II

/I

II

V

II

II

-10 -20

V V V

-40 V

-50 -60

0

0.2

0.4

0.6 ENERGY

v 0.8

1.0

1.2

1.4

Mev

Fig. 1. Effect of size of response matrices on the calculated spectra by the inverse matrix method (different scales are separated by the symbol ~).

to the error of approximation 4 ) of matrix representation of the real continuous functions which are largely magnified by the inversing process. The magnitude of oscillation, therefore, depends on the size of matrix. Large matrix possesses closer approximation to real response function but also induces larger magnification of error in matrix inversing process. For a rather smooth spectrum as the one we assumed, the small matrix gives better result.

2.2.

j= 1,2 ... m,

Xj;;;;O,

V

-30

(4)

LEAST-SQUARES METH00 3)

One of the improvements in unfolding technique is to couple the least-squares criteria into the inverse matrix method by using rectangular response matrix. The simplest way is to assume the statistical weights as constants. In this paper we omit the statistical weight since no statistical error is taken into consideration in

in solving eq. (2) and find out an approxima. solution as close as possible to the exact solution of eq. (2) but simultaneously satisfying the condition x j ;;;; 6If the Xl obtained from certain known Yand A by eq. (3) are positive or zero, the imposed condition xl;;;; 0 has no effect on the solution, hence the result of the inverse matrix method should be the same as that of Scofield's method with an infinitive number of iterations. Table I shows the convergence of iteration TABLE 1

Convergence of Scofield's iteration as

-Assu~:l y

I

xJ~

O.

from Scofield's iteration

J

J

120 160 40 80 eq. (3) \ iterations. iterations iteration~. it~ratio~ ___.--c

_ _ _ _ _. __

2.7977 4.0581 2.7374 1.9944 1.8150 2.0578 2.2823 1.8460 1.6498 1.7630 1.8769 1.6328 0.6063

XJ

~

0.3889 0.4671 0.0000 0.0000 0.0000 0.9889 1.5903 0.6780 1.7503 2.4341 3.9222 3.2864 0.0000

0.3749 0.4481 0.0649 0.0326 0.0542 0.9398 1.6356 0.6203 1.;934 2.3892 3.9482 3.2520 0.0524

---

...

~-

0.3821 0.4571 0.0362 0.0153 0.0290 0.9644 1.6095 0.6552 1.7646 2.4153 3.9291 3.2708 0.0352

0.3844 0.4606 0.0251 0.0097 0.0199 0.9742 1.6029 0.6636 1.7585 2.4223 3.9253 3.2768 0.0252

0.3855 0.4623 0.0192 0.0070 0.0151 0.9762 1.5997 0.6674 1.7561 2.4254 3.9241 3.2794 0.0194

111

STUDY OF SCINTILLATION SPECTROMETRY

,

, 0.5

II

N-41

M-IO

,

N=21

M=IO

o

N-41

M-20

•

N-31

M .. IO

v

N-41

M-30

,

N=31

M=211

•

0.4 - - - TRUE SPECTRUM

>iii

I-

zUJ

0.3

I-

~ UJ

~

I-

0.2

v

,

i

II

el

...J UJ

a:

•

v

0.1

O~~

o

__~~__~__~~~~~~~__~~__~__~~ 0.2

0.4

0.6 ENERGY

0.8 MeV

1.0

1.2

1.4

Fig. 2. Effect of dimensions of response matrices on the calculated spectra by least-squares method.

O.S

•

• •

0.4

I-

j! /'.

20 • 10 30. 10 40 • 10

,

'I

>iii

I-

zUJ

,

LEAST SCOFIELD SQUARE ITERATION

I

:"'---'-"I~:

0.3

~

1

~

I

;::: C

...J

UJ

a:

•

•

0.2

I

0.1

•

I

ENERGY MeV

Fig. 3. Comparison of least-squares method and Scofield's method.

112

Y. S.

TABLE 2 Convergence of Scofields' iteration as negative elements appear. ------

itcr~o~:~r~;!~on-:-Her~~~on~,-~i~e~~~ons_ ~ter~~~~ 0.3883 0.4678 0.0000 0.0000 0.0049 0.9830 1.6006 0.6589 1.7697 2.4184 3.9381 3.2743 0.0000

--

0.3889 0.4671 0.0000 0.0000 0.0001 0.9889 1.5904 0.6775 1.75\0 2.4335 3.9230 3.2857 0.0000

0.3889 0.4671 0.0000 0.0000 0.0000 0.9889 1.5903 0.6780 1.7503 2.4341 3.9223 3.2864 0.0000

0.3889 0.4671 0.0000 0.0000 0.0000 0.9889 1.5903 0.6780 1.7503 2.4341 3.9222 3.2864 0.0000

0.3889 0.4671 0.0000 0.0000 0.0000 0.9889 1.5903 0.6780 1.7503 2.4341 3.9222 3.2864 0.0000

--_._-

-~-~--------

towards the spectrum calculated by the inverse matrix method. If the X obtained from the inverse matrix method has some negative elements, the calculated spectrum of Scofield's iteration will never be the same as that from the inverse matrix method. By several numerical tests performed by computer, it can be easily shown that the calculated results are still convergent to certain finite values. Table 2 shows an example. The question we wish to answer is whether the calculated result converges into the true spectrum. By using the ideal measurement approach we found that the Scofield's iteration does not converge to the true spectrum as shown in fig. 4. In

su other words, an increasing number of iterations does not mean increasing of accuracy of the calculated spectrum. According to our experience, 30 to 60 iterations are good enough to get the best approximation. The case with rectangular response matrix (30 x 20) coupled with the least-squares is not as good as the case with the square matrix (20 x 20). The comparison of different dimensions of the response matrix is shown in fig. 5. It seems that the case 30 x 30 is the best one. Another conclusion which can be drawn is that the coupling of the least-squares technique in Scofield's iteration does not help the results.

2.4. COMPUTER TIME A modification of Scofield's iteration is proposed in order to increase the rate of convergence. By close observation of a real response matrix, we noted, .. that besides the diagonal line, the next strongest line is the line below the diagonal line. Hence if we use both lines instead of the diagonal line only in each iteration, '\he rate of convergence may be largely increased. Table 3 shows our modified scheme and the comparison of our scheme with Scofield's scheme. An arbi trary spectrum is assumed and the corresponding pulse-height distribution obtained according to eq. (2) can then be used to calculate the spectrum by both methods as shown in table 3. It was found that 90 itej'ations or 29 min is required to obtain a calculated spectrum with accuracy up to 1% for all elements

TABLE

3

COMPARISON OF CALCULATING SCHEMES OF BOTH METHODS.

Our modified method

----

'~'---I-~~-

1. Let Xo = Y 2. Put y(m) = AX(m), m = 0, 1,2 ... 3. Solve B(m) from y(m) = B(m)x(m)

Subject to the constraint b~":}I)b~":J = aJ - 1.J/ aJ,J and find out [B(m)] - 1, where bi~) 0 0 ... 0 b~~) b~m/ 0... 0 B(m) = 0 b~m/ b~'"J... 0

o 0 0 ... b~:) and a/.J is the element of A 4. Put X(m+ 1) = [B(m)r 1 y 5. Iteration

Scofield's method

1. Let Xo ::: X VI.) ='~ - , m = 0, 1, 2 ... 2 . Put Y (m) 'i;:,l,;A 3. Solve o(m) from y(m) = o(m)X(m)

and find out [o(m)] - 1, where di~) 0 0 ... 0 o d~~) 0 ... 0 o(m) = 0 0 d~'":/ ... 0

o 4. Put

x(m+l)

5. Iteration

o = [o(m)r 1 y

STUDY OF SCINTILLATION SPECTROMETRY 0.6

•

30 TIMES

N-21 M=20

•

60 TIMES

o

90 TIMES

0.5

o

30 TIMES

N ... 31

•

1'01 .. 20

60 TIMES

•

0.4

j .

TRUE SPECTRUM

r t-t-t/ I - , !-l-·-·/~ ~. ! •

0.3

I.&J

>

~

..J I.&J

It

0

l

Z

I.&J I-

~i

,I:

90 TIMES

>le;;

~

113

•

•

, o

:/!

0.2

I

0.1

o

, 3

•

•

°OO~~~--~n7--~~--~~--L-~--L-~~~~ 0.2 0.4 0.6 0.8 1.0 1.2 1.4 ENERGY MeV

Fig. 4. Convergence of Scofield's iteration. 0.6 j

N-II M-IO N-21

0.5

x 1'01-10

•

N-31 M-IO

"

N-31

o N-41

v M-20

N-41 M-IO

M-20

- - TRUE SPECTRUM 0.4

..

03

,, I

0.2

I

• •

y

0.1

O~~__~~~~~__~~~~__~__~~~~__~~~~

o

0.4

0.6

ENERGY

0.8

1.0

1.2

1.4

MeV

Fig. S. Effect of different dimensions of response matrices on the calculated spectra by Scofield's iteration method.

114

Y. S.

compared with the assumed original spectrum, while 45 iterations or 14 min is required for our proposed method for even be'~ter result. This shows that the use of a complicated form of approximation increases the computing time of each iteration, but the gain obtained from fast convergence is more than the loss required for compensation. 3. New approaches As an experimental factS), one of the bottle-nec,ks of scintillation spectrometry is the difficulty to obtain a good response matrix because of the limited number of monoenergetic sources and uncertainty in interpolating the irregularly and rapidly varying response function. We therefore propose a new approach to establish a new relation between integral quantities of pulse height distribution and responses so that interpolation is made in integral quantities rather than in differential quantities. Furthermore by the assumption of a polynomial approximation of the responses and the spectrum, and under well conditions, polyenergetic sources as well as monoenergetic sources may both be used in calibrating the relation between pulse height distribution and the spectrum. Both approaches are proposed for further study and justification. 3.1. INTEGRAL FORMS OF RESPONSE AND DISTRIBUTION

su as compared to A (V,E) which varies much rapidly with respect to E for each fixed V. Thus, interpolation becomes much easier and more accurate if B(V,E) is used instead of A (V,E). 2. The statistical errors in B(V,E) can be largely reduced, especially for large values of V. 3. Since B( V,E) is less varied and almost monotonic decreasing with respect to E, a numerical solution to the integral equation is easier to be handled and therefore, the accuracy of the solved x(E) can be largely increased. 3.2. POLYNOMIAL APPROXIMATION OF RESPONSES AND SPECTRUM

Tn order to construct a response matrix to a great degree of accuracy a large number of calibrating sources is needed, but the actual available monoenergetic sources are less than, Itay, ten. Hence'the attempt to use polyenergetic sources is desirous. Suppose there are n sources (not necessary to ~be monoenergetic) S1' S2"" Sn. Denote the spectrum of SI by xi(E). If each xl(E) emits a number of gammarays of different energies Ef, ... Ell where k I ~ 1 and if the fractional intensity of the component gamma· ray having energy E} where I ~j ~ kl is O(~, then xl(E) can be represented as:

m-

As the differential response 1 A( E, V)d V and the pulse-height distribution y( V)d V, the integral form of the response g A (V,E)d V for a diltcrete set of E and V and integral form of pulse-height distribution f~Y( V)d V for a discrete set of V can be obtained from a modern multichannel analyzer. If

k,

Xl( E) =

In-I

J: J:

B(V,E) = and

Z(V) =

=

f:

m ••

where b(E - Ej) is the delta function corresponding to

EJ and Let

A(V,E)dV

then

y(V)dV,

B(V,E)x(E) dE.

O(~b( E - E~),

j=1

BI(V)

eq. (1) can be reduced to

Z(V)

L

(5)

The functions Z( V) and B( V,E) are both monotonic increasing with respect to V. Moreover, as an experimental fact, for fixed V, the function B( V,E) is almost monotonic decreasing and sufficiently smooth with respect to E. There are at least three advantages in using eq. (5) instead ofeq. (I) to solve x(E): 1. B( V,E) is much better behaved with respect to E

=

k,

L

j=1

IX~B(V,Ej)'

Each BI( V) can be determined experimentally up to a discrete set of V, V= VI 'V2 ••• Vm for each source. We now approximate each B( Vk,E) by a polynomial of degree n-l, e.i., B( Vk,E) = Pk(E). Then k,

BI(Vk)

= L

j=1

IX~Pk(Ej),

i=I,2, ... n.

This is a system of n linear equations, the n polynomial coefficient!. for Pk(E) may then be determined if the coefficient determinant is not equal to zero.

STUDY OF SCINTILLATION SPECTROMETRY

Let the unknown spectrum be approximated by a polynomial of degree m - I, e.i. m

x(E) =

L: c E, -

1=1

I

1

(6)

•

Substituting eq. (6) into eq. (5) for a discrete set of values V 1 , V 2 ••• Vm, we have

(7) Let

f:mn.

Et - 1 Pi( E) dE

= ri and Z(l'/) = Zi'

then eq. (7), becomes a matrix equation

Z= RC, where Z1

ru

r12··· r 1m

C1

Zz

r Z1

r22 ••• r2m

C2

C=

R=

Z= Z",

rm1

r",2

rmm

Cm

C can be solved by inversing Rand x(E) is then determined by substituting the c/s into eq. (6). There are also three advantages of the approach of polynomial approximation: 1. We may use polyenergetic sources as well as monoenergetic sources for establishing a better relation between pulse-height distribution and energy spectrum. 2. Interpolation is not necessary in this approach,

115

there is no connection between the number of calibrating sources and the number of energy or channel intervals. 3. The final obtained spectrum is in the form of continuous polynomial function which can be used directly in derivation or calculation.

4. Conclusion The result of this paper concludes that there is no method available so that we can calculate the real spectrum very accurately, and the uncertainty of the calculated spectrum is mainly due to the error caused by matrix approximation to the original response equation and the error caused by lack of sufficient original data to calibrate the relation between pulseheight distribution and the energy spectrum. The former can be reduced by iteration or the least square technique, but never be neglected. The latter can be improved by the new approaches we proposed and our recent proposed approach of a linear measurement 6 ). It seems that a method to transform eq. (1) into another form of integral equation which can then be solved by iteration would be a promising route to avoid the matrix approximation.

References 1) J. W. Motz, Phys. Rev. 100 (1955) 1560. 2) N. E. Scofield, NAS-NS 3017, Proc. Symposium (1962) p. 108. S) J. I. Trombka, NAS-NS 3017, Proc. Symposium (1962) p. 183. 4) Y. S. Su and S. M. Shih, Study of scintillation spectrometry: Error analysis, to be publ.

J. C. Chou and S. C. Lin, Nucl. Sci. (China) 4, no. 3 (1965) 39. 6) Y. S. Su and T. C. Wu, Nucl. Instr. and Meth. 41 (1966) 129.

6) Y. S. Su,