SATAN provides means for special processing of analyzer
spectra. Among these are spectrum deconvolution and
automatic peak search.
Deconvolution
The SATAN command ADECON provides the interactive deconvolution of spectra
/MOL61,GOL64,KRA82/. A measured spectrum 'anlin' is expected to be
a convolution of an ideal spectrum 'anlout' with a calculated detector
response function 'anlres'.
anlin = anlres X anloutThe spectrum 'anlout' is calculated by an iterative procedure. The region of the analyzer 'anlin' which is to be deconvoluted must be displayed; the names of 'anlres' and 'anlout' are given within the com- mand. A 'check-analyzer' may be specified to contain the reconvoluted resulting spectrum for comparison with the input spectrum
DELTAE = E(k+1) - E(k)one expects an ideal spectrum f(k) with the observed counts in one channel:
f(k) | +-+ (counts) | | | | | | | | | | | | | | | | | | +------------------------------ channels k E0 (energy E(k))In practice the spectrum y(k) (e.g. gamma-rays) is measured:
tail peak y(k) | (counts) | +-+ | | | | | | |---------+ +-+ +-+ | +-+ | | | +-----+ +-+ +------------------------------- channels k E0 (energy E'(k))Mathematically this can be expressed as:
y(k) = SIGMAl (r(k,l) X f(l))where the summation is performed over index l. r(k,l) is the response matrix. In case the ideal spectrum f contains one discrete energy E0 in channel l0 (delta-function), one column of the response matrix is projected onto vector y:
y(k) = r(k,l0)Hence to construct the matrix the detection system has to be calibrated with monoenergetic sources; a calibration spectrum y(k) measured at the energy E(l) is stored in column l of the response matrix r(k,l) as illustrated in the two-dimensional representation below.
+--------------------------- E(l) |====-- --------------- | -====-- ------------- | -====-- ----------- <- tails | -====-- --------- | -====-- ------- | -====-- ----- | -====-- --- | -====-- - | -====-- | -====-- | -====-- | -====- E'(k) | -=== <- peaksGenerally a sufficient number of sources of different energies to fill the whole matrix is not available. This means that the measured response spectra have to be interpolated. Finally the columns of the response matrix must be normalized to unity, because only on this condition the detector efficiency can be mathematically separated from the total response function. Generating the response function it should be considered that a matrix needs a lot of storage (4 bytes per element).
x(k) = x0 + x1 X kof the input spectrum always corresponds to channel numbers k. If a non-zero offset x0 is specified, an internal bin shift is performed in order to fit to the corresponding response matrix elements which are not allowed to have an energy offset. (I.e., the energy at the low bound of the first bin is expected to be zero. This is due to the fact, that spectra are measured quantities with a 'natural' offset while response functions are generally calculated or fitted and can be provided without offset.) If errors of the input data are given by an analyzer, the same linear calibration is on effect. The returned deconvoluted spectrum and - in case an error analysis is invoked - the output error analyzer imply the same linear calibration as specified for the input spectrum. The calibration of the response matrix may differ from the spectrum calibration with the only restriction that the energy per bin of the matrix elements has to be smaller than that of the spectrum. In case of different calibrations internally a redistribution of the matrix elements is performed by splitting and adding bins in order to fit to the spectrum. This feature allows the user to create the response analyzer once for the detection system used and then apply it to every spectrum measured with it independent of the actual calibration which in general is run-dependent . The constructed response matrix finally is normalized to unity in each column.
fi(k) = c(k) X fi-1(k)where k labels the elements of spectrum y and i is the number of the current approximation to the exact solution. Each element is corrected individually by a correction factor c(k). As zeroeth approximation the input spectrum is taken:
f0(k) = y(k)In every iteration step i a check spectrum
di(k) = SIGMAl (r(k,l) X fi(l))is computed to test the quality of the solution fi(l) by folding it with the matrix. The result should be a good approximation to the input spectrum. An expression for chi-square is evaluated:
chi2 = SIGMAk ((di(k)-y(k)) / DELTAy(k))2 / N(DELTAy = errors of input spectrum, N = number of spectrum elements) which should be of the order of 1 or less. If the relative change of chi-square becomes less than the value in the accuracy parameter the iteration process is stopped even if the maximum number of iteration steps defined by the ITER-parameter is not executed . The algorithm to evaluate the correction factor is based on the formula
c(k) = SIGMAl (r(l,k) X y(l) / di(l)) / SIGMAl r(l,k)where the summation is performed over the index l extending from lmin to lmax which are implicitly determined by the range parameter n to
lmin = k-n , lmax = k+n(lmin and lmax are truncated when k reaches the spectrum bounds). For n>0 this procedure has a smoothing effect, since the neighbouring points are taken into account and a weighting with the matrix elements is done ('response weighted method'). This option works best with continuous spectra where no sharp peaks are expected. The default value of n=2 is an optimum based on experience in deconvoluting beta+- and continuous gamma-spectra. For n=0 the correction factor is equal to
c(k) = y(k) / di(k)This is known as the 'quotient method'. It is useful for spectra with a lot of peaks where the intensity from the tail contributions shall be stored into separate peaks. It is possible that even small 'hidden' peaks hooked on the shoulder of dominating peaks (e.g. conversion electron spectra where an L-line may be located very close to the K-line) appear after applying this deconvoluting technique. On the other hand the quotient method may interprete statistical fluctuations in the spectrum as nearly hidden peaks and emphasize them during deconvolution . If the conformity between refolded and measured spectrum is found to be too poor (chi2>1) an iteration step is calculated from the derivatives of chi2 with respect to the f(k) ('gradient search algorithm'). Since the evaluation of the gradient is time-consuming for large spectra (more than 100 bins) the number of steps is limited automatically depending on the spectrum size. Negative elements of the input spectrum are treated as zero after con- firmation.
Last updated: M.Kraemer@gsi.de, 3-Aug-1999