TSpectrum

Background estimation

Goal: Separation of useful information (peaks) from useless information (background)

• method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping algorithm

• new value in the channel “i” is calculated

where p = 1, 2, …, numberIterations. In fact it represents second order difference filter (-1,2,-1).

Function:

const char* Background(float *spectrum, int ssize, int numberIterations, int direction, int filterOrder, bool smoothing, int smoothingWindow, bool compton)

This function calculates background spectrum from the source spectrum. The result is placed in the vector pointed by spectrum pointer. One can also change the direction of the change of the clipping window, the order of the clipping filter, to include smoothing, to set width of smoothing window and to include the estimation of Compton edges. On successful completion it returns 0. On error it returns pointer to the string describing error.

Parameters:

spectrum-pointer to the vector of source spectrum

ssize-length of the spectrum vector

numberIterations-maximal width of clipping window,

direction- direction of change of clipping window

- possible values=kBackIncreasingWindow

kBackDecreasingWindow

filterOrder-order of clipping filter,

-possible values=kBackOrder2

kBackOrder4

kBackOrder6

kBackOrder8

smoothing- logical variable whether the smoothing operation in the estimation of

background will be incuded

- possible values=kFALSE

kTRUE

smoothWindow-width of smoothing window,

-possible values=kBackSmoothing3

kBackSmoothing5

kBackSmoothing7

kBackSmoothing9

kBackSmoothing11

kBackSmoothing13

kBackSmoothing15

compton- logical variable whether the estimation of Compton edge will be incuded

- possible values=kFALSE

kTRUE

References:

[1] C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the quantitative analysis of PIXE spectra in geoscience applications. NIM, B34 (1988), 396-402.

[2] M. Morháč, J. Kliman, V. Matoušek, M. Veselský, I. Turzo.: Background elimination methods for multidimensional gamma-ray spectra. NIM, A401 (1997) 113-132.

[3] D. D. Burgess, R. J. Tervo: Background estimation for gamma-ray spectroscopy. NIM 214 (1983), 431-434.

Smoothing

Goal: Suppression of statistical fluctuations

• the algorithm is based on discrete Markov chain, which has very simple invariant distribution

being defined from the normalization condition

n is the length of the smoothed spectrum and

is the probability of the change of the peak position from channel i to the channel i+1. is the normalization constant so that and m is a width of smoothing window.

Function:

const char* SmoothMarkov(float *spectrum, int ssize, int averWindow)

This function calculates smoothed spectrum from the source spectrum based on Markov chain method. The result is placed in the vector pointed by source pointer. On successful completion it returns 0. On error it returns pointer to the string describing error.

Parameters:

spectrum-pointer to the vector of source spectrum

ssize-length of the spectrum vector

averWindow-width of averaging smoothing window

Reference:

[1] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 (1996), 451.

Deconvolution

Goal: Improvement of the resolution in spectra, decomposition of multiplets

Mathematical formulation of the convolution system is

where h(i) is the impulse response function, x, y are input and output vectors, respectively, N is the length of x and h vectors. In matrix form we have

• let us assume that we know the response and the output vector (spectrum) of the above given system.

• the deconvolution represents solution of the overdetermined system of linear equations, i.e., the calculation of the vector x.

• from numerical stability point of view the operation of deconvolution is extremely critical (ill-posed problem) as well as time consuming operation.

• the Gold deconvolution algorithm proves to work very well, other methods (Fourier, VanCittert etc) oscillate.

• it is suitable to process positive definite data (e.g. histograms).

Gold deconvolution algorithm

where L is given number of iterations (numberIterations parameter).

Boosted deconvolution

1. Set the initial solution

2. Set required number of repetitions R and iterations L

3. Set r = 1.

4. Using Gold deconvolution algorithm for k=1,2,...,L find

5. If r = R stop calculation, else

a. apply boosting operation, i.e., set

i=0,1,...N-1 and p is boosting coefficient >0.

b. r = r + 1

c. continue in 4.

Function:

const char* Deconvolution(float *source, const float *respMatrix, int ssize, int numberIterations, int numberRepetitions, double boost)

This function calculates deconvolution from source spectrum according to response spectrum using Gold deconvolution algorithm. The result is placed in the vector pointed by source pointer. On successful completion it returns 0. On error it returns pointer to the string describing error. If desired after every numberIterations one can apply boosting operation (exponential function with exponent given by boost coefficient) and repeat it numberRepetitions times.

Parameters:

source-pointer to the vector of source spectrum

respMatrix-pointer to the vector of response spectrum

ssize-length of the spectrum vector

numberIterations-number of iterations (parameter l in the Gold deconvolution

algorithm)

numberRepetitions-number of repetitions for boosted deconvolution. It must be

greater or equal to one.

boost-boosting coefficient, applies only if numberRepetitions is greater than one.

Recommended range <1,2>.

References:

[1] Gold R., ANL-6984, Argonne National Laboratories, Argonne Ill, 1964.

[2] Coote G.E., Iterative smoothing and deconvolution of one- and two-dimensional elemental distribution data, NIM B 130 (1997) 118.

[3] M. Morháč, J. Kliman, V. Matoušek, M. Veselský, I. Turzo.: Efficient one- and two-dimensional Gold deconvolution and its application to gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.

[4] Morháč M., Matoušek V., Kliman J., Efficient algorithm of multidimensional deconvolution and its application to nuclear data processing, Digital Signal Processing 13 (2003) 144.

Richardson-Lucy deconvolution algorithm

· for discrete systems it has the form

· for positive input data and response matrix this iterative method forces the deconvoluted spectra to be non-negative.

· the Richardson-Lucy iteration converges to the maximum likelihood solution for Poisson statistics in the data.

Function:

const char* DeconvolutionRL(float *source, const float *respMatrix, int ssize, int numberIterations, int numberRepetitions, double boost)

This function calculates deconvolution from source spectrum according to response spectrum using Richardson-Lucy deconvolution algorithm. The result is placed in the vector pointed by source pointer. On successful completion it returns 0. On error it returns pointer to the string describing error. If desired after every numberIterations one can apply boosting operation (exponential function with exponent given by boost coefficient) and repeat it numberRepetitions times (see Gold deconvolution).

Parameters:

source-pointer to the vector of source spectrum

respMatrix-pointer to the vector of response spectrum

ssize-length of the spectrum vector

numberIterations-number of iterations (parameter l in the Gold deconvolution

algorithm)

numberRepetitions-number of repetitions for boosted deconvolution. It must be

greater or equal to one.

boost-boosting coefficient, applies only if numberRepetitions is greater than one.

Recommended range <1,2>.

References:

[1] Abreu M.C. et al., A four-dimensional deconvolution method to correct NA38 experimental data, NIM A 405 (1998) 139.

[2] Lucy L.B., A.J. 79 (1974) 745.

[3] Richardson W.H., J. Opt. Soc. Am. 62 (1972) 55.

Unfolding

Goal: Decomposition of spectrum to a given set of component spectra

Mathematical formulation of the discrete linear system is

Function:

const char* Unfolding(float *source, const float **respMatrix, int ssizex, int ssizey, int numberIterations, int numberRepetitions, double boost)

This function unfolds source spectrum according to response matrix columns. The result is placed in the vector pointed by source pointer. The coefficients of the resulting vector represent contents of the columns (weights) in the input vector. On successful completion it returns 0. On error it returns pointer to the string describing error. If desired after every numberIterations one can apply boosting operation (exponential function with exponent given by boost coefficient) and repeat it numberRepetitions times. For details we refer to [1].

Parameters:

source-pointer to the vector of source spectrum

respMatrix-pointer to the matrix of response spectra

ssizex-length of source spectrum and # of columns of the response matrix

ssizey-length of destination spectrum and # of rows of the response matrix

numberIterations-number of iterations

numberRepetitions-number of repetitions for boosted deconvolution. It must be

greater or equal to one.

boost-boosting coefficient, applies only if numberRepetitions is greater than one.

Recommended range <1,2>.

Note!!! sizex must be >= sizey After decomposition the resulting channels are written back to the first sizey channels of the source spectrum.

Reference:

[1] Jandel M., Morháč M., Kliman J., Krupa L., Matoušek V., Hamilton J. H., Ramaya A. V.: Decomposition of continuum gamma-ray spectra using synthetized response matrix. NIM A 516 (2004), 172-183.

Peaks searching

Goal: to identify automatically the peaks in spectrum with the presence of the continuous background and statistical fluctuations - noise.

The common problems connected with correct peak identification are

non-sensitivity to noise, i.e., only statistically relevant peaks should be identified.
non-sensitivity of the algorithm to continuous background
ability to identify peaks close to the edges of the spectrum region. Usually peak finders fail to detect them
resolution, decomposition of doublets and multiplets. The algorithm should be able to recognize close positioned peaks.
ability to identify peaks with different sigma

Function:

Int_t SearchHighRes(float *source,float *destVector, int ssize, float sigma, double threshold, bool backgroundRemove,int deconIterations, bool markov, int averWindow)

This function searches for peaks in source spectrum. It is based on deconvolution method. First the background is removed (if desired), then Markov smoothed spectrum is calculated (if desired), then the response function is generated according to given sigma and deconvolution is carried out. The order of peaks is arranged according to their heights in the spectrum after background elimination. The highest peak is the first in the list. On success it returns number of found peaks.

Parameters:

source-pointer to the vector of source spectrum

destVector-resulting spectrum after deconvolution

ssize-length of the source and destination spectra

sigma-sigma of searched peaks

threshold- threshold value in % for selected peaks, peaks with amplitude less than threshold*highest_peak/100 are ignored

backgroundRemove- background_remove-logical variable, true if the removal of background before deconvolution is desired

deconIterations-number of iterations in deconvolution operation

markov-logical variable, if it is true, first the source spectrum is replaced by new spectrum calculated using Markov chains method

averWindow-width of averaging smoothing window

Fig. 27 An example of one-dimensional synthetic spectrum with found peaks denoted by markers

References:

[1] M.A. Mariscotti: A method for identification of peaks in the presence of background and its application to spectrum analysis. NIM 50 (1967), 309-320.

[2] M. Morháč, J. Kliman, V. Matoušek, M. Veselský, I. Turzo.:Identification of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000) 108-125.

[3] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 (1996), 451.

TSpectrum

class TSpectrum : public TNamed

Data Members

Class Description

Inline Functions