SATAN special data evaluation procedures
| GSI | Biophysics | SATAN long write-up |

SATAN special data evaluation procedures

NOTE: the features here are not yet implemented in this version.

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 anlout
The 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

The response matrix

Theoretically the energy distribution of the radiation of a monoenergetic source is a delta-function. Subdividing the range of energy in small intervals
    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) |                       -===  <- peaks
Generally 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).

Calibration

A user-defined linear calibration
 
    x(k) = x0 + x1 X k
of 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.

The deconvolution procedure

The following section describes the principle of operation of the decon- volution procedure. For detailed information refer to the original pub- lications. The deconvolution works iteratively according to the following equation:
 
    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.

Errors (uncertainties)

Experimental errors of the input data may be specified via a particular analyzer; by default statistical errors (square roots of the observed counting rates) are assumed. An error analysis for the deconvoluted spectrum is performed if the ERRANL-keyword (with an analyzer name as argument) is specified. The familiar law of quadratic propagation of errors is not applicable since correlations between the variables involved (original spectrum elements y(k) and the elements of the previous iteration fi-1(k)) are neglected. An estimation for the error DELTAf(k) of a channel contents f(k) is given by variing f(k) holding all other channel contents constant until chi2 increases by 1.0 (i.e. one standard deviation). This means, the chi2-equation is solved for DELTAf(k) after replacing f(k) by f(k)+DELTAf(k); the procedure is repeated throughout the whole spectrum.

Problems

The proper execution of deconvolution is strongly dependent on the input data. The following difficulties may occur.

Peak search

The commands DFPEAK (to mark and list peaks) and FPEAKS (to define Gaussian or Lorentzian peaks as a fit function) are provided for an automatic search of any number of peaks within a displayed spectrum. The utility is also available as procedure $DFPEAK to be called from a user written program. The displayed data are smoothed internally over the specified width and convoluted with a bipolar function to a spectrum, in which zero crossings correspond with relative minima and maxima of the original spectrum. The net peak area estimated from the convolution function is divided by the square root of the corresponding contents of the smoothed spectrum to define the statistical peak significance. Peaks with a significance smaller than the specified value are ignored; if more peaks than desired are found, only the most significant ones are considered. A reasonable guess for the full peak width at half maximum (FWHM) is essential. A small value may overestimate statistical fluctuations; if chosen too large, sharp peaks or peaks located in a shoulder may be smoothed off.
| GSI | Biophysics | SATAN long write-up |

Last updated: M.Kraemer@gsi.de, 3-Aug-1999