NUCLEAR
INSTRUMENTS
A N D M E T H O D S 98 (1972) I 3 I - I 3 4 ;
© NORTH-HOLLAND
PUBLISHING
CO.
ANALYSIS OF T H E M A T H E M A T I C A L M E T H O D S F O R C O R R E C T I N G C O N T I N U O U S X-RAY SPECTRA O B T A I N E D BY S C I N T I L L A T I O N S P E C T R O M E T R Y M. P A G A N I N I FIORATTI and S. RICCI PIERMATTEI Laboratorio di Dosimetria e Standardizzazione, C N E N - CSN, Casaccia, Roma, Italy
Received 2 June 1971 An analysis is made of the analytical methods for interpreting the Bremsstrahlung spectra. Owing to the nature of the response matrix, the iterative method is the most adequate to obtain correct information.
1. Introduction The handling of data from X-ray spectra measurements has been widely discussed in the literature as far as the line spectra are concerned; the information available for Bremsstrahlung spectra are, on the contrary, rather p o o p - 4 ) . The aim of this paper is to analyse the mathematical techniques to be used to handle data from continuous spectra. The smoothing of experimental data and the correction due to the finite resolution of the detector will be examined in detail. 2. Smoothing of the experimental data The smoothing of the experimental data can be accomplished by hand. This procedure is rather subjective and, when using a multichannel analyser with more than one hundred channels, a hand procedure is impossible. An attempt is made to understand the random fluctuations recorded in the count number per channel as noise superimposed on the original data, and to filter such a noise by a Fourier analysis similar to that often u,;ed in communication theory to separate the signal from the noiseS). The experimental spectrum NdE ) can be represented as the sum of two components:
N~(E) =
Nv(E ) + R(E),
where Nv(E ) is the original spectrum and R(E)the noise. I f N£(to) and R'(to) are the Fourier transforms of Nv(E) and R(E), respectively, we have: N~(o) = N ; ( o ) + R'(to). The success of the method relies upon the possibility of choosing a filter function which eliminates, at least, part of the noise: this generally happens for line spectra where Nv(to) is spread over a certain number of channels and R'(to) all over the spectrum. 131
2O
6 E
8 15
lb
go
1~o
KeY
1~o
-
Fig. 1. X-ray spectrum- 150keVmaximumenergy. Experimental spectrum (histogram); smoothed spectrum (continuous line).
As we are dealing with a continuous spectrum, this condition is not fulfilled and, therefore, we applied a standard smoothing method using the summation formula6). In deriving this formula it is assumed the function can be represented by a third degree polynomial over the range of the observations used to calculate a smoothed value. The optimum number of points for smoothing is linked to the shape of the function we are to fit. If the number is too large the smoothing procedure flattens the peaks and fills in the valleys and
132
M. P A G A N I N I
FIORATTI
the original spectrum becomes distorted, if this n u m b e r is too small the r a n d o m fluctuations are not removed from the smoothing procedure and the significant features of the spectrum are lost. We reasonably assumed a number of channels of the same order of magnitude of the energy resolution of the detector in the region under study. In fig. 1 we report as an example an experimental spectrum and a smoothed one.
3. Resolution correction As the resolving power of the detector is not zero, the spectra are only partially resolved. If h ( E - E ' ) is the response function of the detector and N(E') the incident spectrum, the recorded spectrum is given by:
M(E) =
I S h(E-E')N(E')dE', o
a Fredholm integral equation of the first kind. The approximate solution can be obtained expressing the integral as a n terms sum:
M(E,) = E N(E'k) h(Ei- E'k)AE. k The numerical procedures for solving linear systems are direct or indirect v' 8). Direct methods would yield an exact solution in a finite number of operations if there were no r o u n d o f f errors. Indirect or iterative processes begin with an approximate solution which is improved with each step of the iteration. The steps required to obtain an exact solution would be infinite in the absence of r o u n d o f f errors. 3.1. DIRECT METr~OD The direct method could appear at a first sight as the most promising one. The availability of large computers allows practically the solution of any system of equations in any number of unknowns. The difficulties which arise depend on: 1) r o u n d o f f errors performed during the calculation; 2) experimental errors which affect the experimental data. It can be shown 8) if the eigenvalues 21 of the matrix are different f r o m zero but very close to zero an error ei affecting the counts Mi produces an error qi on the N i c o m p o n e n t o f the solution vector given by: rh = (1/,t3e~; for example if 2 ~ 1 0 -3 an error of 1% on M i gives an error of 1000% on the solution vector. This happens also when the matrix inversion is correct. This is the why sometimes the results obtained
AND
S, R I C C I
PIERMATTEI
by solving a linear system are meaningless f r o m a physical point of view. Of great importance in studying the behaviour of the solution of a system o f equations is the value of the condition number n, i.e. the ratio between 2ma× and 2m~, n =
I;-,.,~,,l;~m.d
•
When a maxtrix has all its eigenvalues close together so that I.~mad2mi,[~l, the matrix is said to be well conditioned for solving linear equations. The condition number is a measure of the effect on an approximate inverse of a matrix of the r o u n d o f f errors and of the experimental errors affecting the data. 3.2. I N D I R E C T M E T H O D As we said above the indirect method uses an approximate solution which is improved by successive iterations. Given the linear system:
Ax=b,
(1)
let us define an iteration of degree k, a function o f the form: x (k+') = Fk(A,b,x(k)), (2) which can be written as:
x(k+ 1)
=
X(k) -k c(k) p(k);
(3)
the matrices C (k) and p(k) are functions of k, A, b and
x(k). A representation of the type of (3) must leave the solution vector h = A - ~b invariant, that is, ifx (k) = A - 1b, then we should obtain x(k+ 1)= A-lb, hence c(k)p(k)=o. We define kth residual vector, r(k)= b--Ax (k) and kth error vector e (k) = h - x (~). A n iteration converges if and only if the eigenvalues the (I-CkA) are less than 1 in absolute value. As the iterative m e t h o d uses the matrix in its original form, in each iteration there will be a correction effect of the matrix on itself and, as consequence, the r o u n d o f f errors are minimized. Moreover if the matrix has some zero elements these can be used to speed up the calculations.
4. Application of the present analysis to the experimental data This analysis was applied to a typical Bremsstrahlung spectrum given by a commercial X-ray tube, as detector a NaI(T1) crystal has been used. In this case the response function h ( E - E ' ) is of Gaussian type:
exp[-(E-E')Z/2a2];
CORRECTING
CONTINUOUS
X-RAY
133
SPECTRA
TABLE 1 Comparison between direct and indirect method for correcting X-ray spectra. 2O
Energy (keV)
Experimental spectrum
Direct method
Indirect method (iterative method)
d t~
18.38 22.42 26.46 30.50 34.54 138.58 42.62 46.66 50.70 54.74 58.78 ,62.82 66.86 70.90 74.94 78.98 :33.02 87.06 91.10 '95.14 '99.18 103.22 107.26 111.30 115.34 119.38 123.42 1127.46 1.31.50 135.54
0.1748 E4 0.2881 E4 0.3728 E4 0.3773 E4 0.3100 E4 0.2344 E4 0.2135 E4 0.2806 E4 0.4577 E4 0.7553 E4 0.1134 E5 0.1496 E5 0.1738 E5 0.1819 E5 0.1772 E5 0.1667 E5 0.1555 E5 0.1449 E5 0.1341 E5 0.1222 E5 0.1089 E5 0.9541 E4 0.8283 E4 0.7159 E4 0.6145 E4 0.5174 E4 0.4195 E4 0.3224 E4 0.2330 E4 0.1572 E4
0.8160 E3 0.1765 E4 0.3784 E4 0.3676 E4 0.2346 E4 0.1961 E4 0.2056 E4 0.2743 E4 0.4068 E4 0.6199 E4 0.1153 E5 0.1562 E5 0.1852 E5 0.1874 E5 0.174l E5 0.1729 E5 0.1532 E5 0.1452 E5 0.1403 E5 0.1211 E5 0.1052 E5 0.1063 E5 0.7423 E4 0.7644 E4 0.5536 E4 0.6120 E4 0.3606 E4 0.3853 E4 0.1501 E4 0.2482 E4
0.1378 E4 0.2468 E4 0.3380 E4 0.3559 E4 0.2998 E4 0.2289 E4 0.2080 E4 0.2729 E4 0.4474 E4 0.7451 E4 0.1129 E5 0.1499 E5 0.1747 E5 0.1828 E5 0.1778 E5 0.1669 E5 0.1555 E5 0.1450 E5 0.1344 E5 0.1225 E5 0.1091 E5 0.9545 E4 0.8280 E4 0.7157 E4 0.6152 E4 0.5186 E4 0.4204 E4 0.3222 E4 0.2334 E4 0.1654 E4
depends o n the detector resolving power. The system to be solved is:
M(Ei) = • N(Ek)exp[- (El- Ek)Z/2~r~]AE.
?, o 15 ¸
10
5
I~
5~
i~o
KeV
I~o
•
Fig. 2. X-ray spectrum - 150 keV maximum energy. Corrected for resolution by iterative method. case. This result has a general validity as the matrix is quasi-diagonal, i.e. only the elements a r o u n d the principal diagonal are different f r o m zero, a n d it can be shown 7, 8) that matrices of this kind are always ill c o n d i t i o n e d matrices. Subsequently we applied the iterative method. I n this case the experimental spectrum was used as the zero order a p p r o x i m a t i o n a n d the corrected first order spectrum M I ( E ) is given by:
k
At first we try the solution using the direct method. In order to reduce at m i n i m u m r o u n d o f f errors connected with the matrix inversion, we solved the system using the C h o l e s k i - T u r i n g m e t h o d 7, 9). I n this m e t h o d the matrix is t r a n s f o r m e d in the p r o d u c t of two t r i a n g u l a r matrices and the system is solved by a substitution method. [n table 1 we report the correct spectrum. It is possible to see how the solution oscillates. A n ewaluation of the c o n d i t i o n n u m b e r of the matrix by soi[ving, with the K r y l o v method9), the secular e q u a t i o n associated to the matrix itself, gives: n ~ 1 0 3. Therefore the direct m e t h o d c a n n o t be applied in this
M(E') h(E- E')dE' ;
Mx(E ) = 0
the second order a p p r o x i m a t i o n is given by: M2(E ) =
M(E) - [M, (E) -
M(E)] = 2M(E) - M , ( E ) ;
similarly the spectra M3(E), M4(E), etc. are obtained. The process converges as the eigenvalues of the matrix are in absolute value less t h a n 1. The iterations are stopped when the difference between the experimental points a n d the correct ones, distorted according to the detector resolving power, are of the same order of m a g n i t u d e of the statistical fluctuations. I n fig. 2 we report the spectrum o b t a i n e d by the iterative method.
134
M. P A G A N I N I
FIORATTI
AND
S. R I C C I
PIERMATTEI
5. Conclusions
References
T h e c o n c l u s i o n s we m a y d r a w f r o m o u r analysis are the f o l l o w i n g :
i) G. Hettinger and V. Starfelt, Arkiv Fysik 14 (1959) 497. 2) M. Erlich, J. Res. Nat. Bur. Std. 54 (1955) 107; RP 2339. :3) L. D. Skarsgard et al., Radiation Res. 14 (1961) 261. 4) E. Casnati et al., CNEN RT/PROT (1969) 16. '~) T. Harper et al., MIT 3944-2 (1968). 6) E. Whittaker and G. Robinson, The calculus of observations (Blackie and Sons Ltd., London, 1960). 7) j. R. Westlake, Handbook of numerical matrix inversion and solution of linear equations (J. Wiley, New York, 1968). s) C. Lanczos, Linear differential operators (Van Nostrand Ltd., New York, 1961). ~'~)V. N. Faddeeva, Computational methods of linear algebra (Dover, New York, 1959).
- A p p l i c a t i o n o f the direct m e t h o d is n o t c o n v e n i e n t as the m a t r i x is quasi d i a g o n a l a n d therefore is ill c o n d i t i o n e d . T h e s o l u t i o n s o b t a i n e d c a n be, therefore, m e a n i n g l e s s . - T h e a p p l i c a t i o n o f the iterative m e t h o d is strongly r e c o m m e n d e d as one c a n t a k e a d v a n t a g e o f the zeros p r e s e n t in the m a t r i x , a n d the rate o f c o n vergence is r a t h e r fast.