A baseline estimation procedure to improve MDA evaluation in gamma-ray spectrometry

The evaluation of minimum detectable activity (MDA) for a radionuclide in a gamma-ray spectrum is generally carried out through the computation of a suitable background count. This task is sometimes difficult for complex spectra for the presence of many photopeaks which make the trend of continuum extremely variable due to multiple dispersion effects and interference factors. It follows that the MDA assessment must be take into account the contributions of all gamma emissions of radionuclides contained in a sample and its value can be significantly higher than that determined by considering only the background of the spectrometric system due to the overlapping of other peaks. A procedure or an algorithm to determine, each time, the count values to be used for the calculation of MDA is interesting and useful. In this work, some of the more recent algorithms proposed for background subtraction in a gamma-ray spectrum have been examined, applying them in an inverse way for the evaluation of baseline trend in the whole energy range. Among the algorithms examined, particular attention was paid to the application of SNIP (statistical sensitive nonlinear iterative peak clipping) algorithms, which are the simplest to adopt and implement in an application procedure. The results obtained in the analysis of test gamma-ray spectra are satisfactory and allow to quickly determine the MDA values with a formulation based on the ISO-11929 standard.


Introduction
Gamma-ray spectrometric analysis with semiconductor detectors is an increasingly used technique for various experimental applications such as environmental monitoring, radiation protection surveillance, food qualification, material characterization, and others.This development is mainly due to its intrinsic characteristics, energy resolution, high efficiency of the detectors and the availability of automatic analysis codes that makes it easy to use in many applications and scientific research activities.
An important parameter for spectrometric analysis used in radiation protection monitoring activity is the Minimum Detectable Activity (MDA), i.e., the smallest amount of a radionuclide which in one measurement can be distinguished from the background with some degree of confidence.For this aim, it is necessary the knowledge of the main measurement parameters such as an efficiency value and the determination of a relative "background" count.This term refers to the continuous base (baseline) on which the photoelectric peaks are superimposed, mainly due to Compton interactions of higher energy gamma photons, incomplete collection of charge in the detector, photons scattered by the surrounding materials and other external factors, such as emissions from radionuclides in the environment, cosmic rays and noise from electronic components.
Generally, MDA is a quantity that can be evaluated "a priori" [1,2], through the measurement of a "blank", i.e., a matrix having the same dimensions, composition and density as the sample under examination in which the activity of the radionuclide of interest is not significant.However, this MDA value cannot always be adopted for all measurements due to the variability of the background components.In routine measurements, often characterized by small quantities of radioactivity in the sample, MDA values substantially does not vary: conversely, in spectrometric measurements of variously composed samples, MDA values can be influenced significantly by the radionuclide activities in the sample and their associated components (continuum and photopeaks).Evaluation of MDA in complex gamma-ray spectra characterized by multiplets composed of more overlapping photopeaks can be a very difficult task.It is then interesting and useful to have a simple and suitable algorithm to provide the trend of the gamma-ray spectrum baseline for evaluating the MDA value for a radionuclide in each spectrometric measurement.
In past decades, several algorithms have been proposed for the evaluation of the gamma-ray spectrum baseline trend, although most are oriented to optimize the gamma peaks automatic search procedure.A concise but also detailed review of the main background subtraction methods, listed on the basis of different used algorithms, is given in [3].Some of these baseline evaluation techniques were previously analyzed to verify the reliability of the procedures.The results were acceptable for relatively simple spectra, i.e., with isolated and spaced peaks, while the presence of complex multiplets makes the results unreliable.The main conclusion of the analysis is that a valid method applicable to all situations could not be easily identified, leaving the choice of the most reliable method to the judgement of an expert researcher.Among all methods, a particular interest has been directed to SNIP (Sensitive Nonlinear Iterative Peak-clipping) technique due to its easy computation of values and the possibility to implement the corresponding routines in a normal spreadsheet.
The SNIP algorithm, initially proposed in [4], is considered the most efficient method for background subtraction in order to optimize the search for the photopeaks present in a gamma-ray spectrum.In this work, it is applied in an inverse mode, i.e., for the estimation of the baseline trend in the whole gamma-ray spectrum energy range.The results of the evaluation depend significantly on the spectrum parameters and can lead to different results in relation to the assumed values [4].For this reason, several improvements have been proposed, almost all concerning the evaluation of the width of the peak region in relation to its variability with calibration parameters [5][6][7][8].
In this work, after a brief description of the SNIP iterative technique, some applications of the algorithm are presented, taking into consideration various gamma-ray spectra detected under different conditions.The results obtained confirmed the goodness and reliability of the technique which is therefore proposed for baseline estimation, useful for calculating MDA value for a radionuclide in a fairly complex spectrum.

Materials and methods
The SNIP procedure consists of an iterative approximation of the content of a channel of a gamma-ray spectrum through substitution with the minimum value between the average value of the gradually equidistant channel contents and the channel content itself.
To reduce the influence of statistical uncertainties and to take into account also weak peaks in a gamma-ray spectrum, the experimental data, i.e., the channel number and its content (i, y i ), are normally transformed in a data couple (i, z i ) using an operator LLS (two log operators plus square root operator) with the relation Therefore, in addition to the vector of the experimental data with N the maximun number of channels of the gamma-ray spectrum, another working vector is available To apply the iterative procedure described below, a region of interest (ROI) in which one or more photopeaks are present must be identified.Different methods allow to identify the start and end channels of a ROI and, even in the absence of a calibration curve, to determine its width in terms of channels.In this work, a suitable "zero area" filter possibly integrated with a peak search technique is used.For most ROIs, only a rectangular filter with a mathematical description formulated as follows is needed [9]with and n 2 × INT(FWHM(i) + 1), where FWHM(i) is the "Full Width at Half Maximum" at channel i.The presence of a probable photopeak is detected if with δ a suitable threshold value (generally selectable in a range 0 ÷ 4) and σ (F d, i ) the standard deviation evaluated by applying the error propagation as The choice of a very low δ value extends the application of the algorithm also to very small peaks, while a very high value limits the application to photopeaks with higher channel content.Other digital filters (triangular, Gaussian,….) can be used without significant differences.
The application of a filter in some cases is not sufficient, especially in the presence of complex multiplets.In these cases, the ROIs selection procedure can be integrated by a photopeak search based on the study of the first derivative y i .Its value can be compared, channel by channel, with a threshold given by the product of the square root of y i by a factor f i in which the choice of the value of s, named "sensitivity factor", is important for the recognition of small peaks in the spectrum.Verifying the relationship (5), it is possible to identify the channels i' and i' + 1 where the indicative conditions of the presence of a photopeak are realized and determining a value of an hypothetical center channel All the peaks with center distance less than 5 × FWHM(i) away from each other, considered partially overlapping, are grouped in a single ROI.The extreme channels of the interval represent the limits of the group and their distance in terms of channels represents the width w i to be considered for the ROI.This analysis leads to identify ROIs with high probability of the presence of photopeaks with respect to others representative of local continuum, information that can be encoded in another auxiliary vector R 0, 0, 0, . . .w 1 , w 1 , w 1 , . . ., 0, 0, . . .w h , w h , w h , . . ., w q , w q , w q , . . .0, 0, 0 where 0 is assigned to a continuum channel and w h the width of the h-th ROI among those identified.The representative parameter of the width of the region, m, can therefore be evaluated via the width of the ROI w h , as 2m + 1 w h (8) with h 1,.., q, if q is the number of identified ROIs.
When the gamma spectrum is simple, i.e., composed of single and spaced peaks, the width of the ROI can be placed in relation to the FWHM(i) value for the channels of interest through an appropriate proportionality factor t which generally assumes values in the range 2-3 [7].
After establishing the value of m for each ROI, the procedure computes the value z i for the p-th iteration using the expression: where p increases iteratively from 1 to the value m and r i is the element of the vector R. Once the final z i value is obtained, the corresponding baseline value, b i , is computed with the anti-transformation the trend of which can be assumed as representative of the spectrum baseline (a fourth working vector B).Generally, b i values are subjected to a further smoothing process to attenuate statistical fluctuations.Furthermore, to avoid unnecessary calculation, the iterative process can be interrupted when the baseline value calculated in the p-th iteration with respect to the value of the previous one differs by a percentage factor ε as small as desired [7] Figure 1 shows a graphical representation of the SNIP iterative procedure and its implementation with the four working vectors used.
For the study of the application of the SNIP algorithm, some gamma-ray spectra detected on samples of different nature, various experimental conditions have been considered.The variety of detectors, listed in Table 1, due to their characteristics and fields of application, seems to constitute a sufficiently valid test basis for evaluating the method application.
Finally, for the calculation of the MDA, the relationship in the form elaborated in [10] with reference to the ISO-11929 standard can be used, having indicated with B the background area subtended to the peak of interest and k a value relative to the chosen degree of confidence.Considering the "no-peaked background" formulation, we have that with ω ε(E)I (E) T c , ε(E) photoelectric efficiency (FEPE, full-energy-peak efficiency), I(E) emission intensity for the energy E of the radionuclide for which MDA has been computed, T c the live counting time.The variance of ω can be easily calculated taking (2023) 138:700  into account the propagation of the uncertainties on FEPE value and emission intensity of the reference gamma line as reported in radionuclide data libraries (i.e., [11]).Assuming k 1.645 (confidence interval 2σ ) we also obtain the relation similar to that proposed by [1] for "paired" measurements, i.e., carried out for the same live counting time for the sample and the "blank", except for the introduction of the variance of the parameter ω.Finally, the relationship is modified by accepting Brodsky suggestion [12] to replace the value 2.71 with 3 and calculating the background value B in a channel interval of 2.5 × FWHM (± 1.25 FWHM) around the channel corresponding to energy E of interest.

Results and discussion
The algorithm has been tested using a common Microsoft Excel spreadsheet and some routines implemented in Visual Basic.The value of the representative parameter m is particularly critical for the results of the iterative process.A rather small m value with respect to the FWHM causes the overestimation of the baseline and consequently an increase in MDA value, while a too large m value introduces distortions in the baseline estimation.If the FWHM calibration curve is known, the value of m can be set equal to a value approximately equal to 1.25 × FWHM, except for multiplets, otherwise the value m must be evaluated from the width of the ROIs as previously defined.In Fig. 2 is represented, as an example, the variation of the baseline for a calibration peak at 1332 keV of 60 Co as function of m values.It is noted that the m values from 4 to 12 are small for representing photopeak width and therefore In Fig. 3 is represented a complex multiplet in which the definition of the peak region has been achieved using both the digital filter and the above-indicated center search procedure.While the filter fails to accurately identify the limits of the regions, the second procedure identifies several regions with more or less complex multiplets with values of m 18, m 24 or m 30 to be compared with the values for a single photopeak, m 7 or m 8 with respect to their different energy ranges.Fig. 3 Baseline trend in a gamma spectrum of a sample with natural radioactivity.Live Counting time: 3000 s.In the insert is reported, as a comparison, the linear background trend obtained with the use of Gamma Vision routines (green line) by ORTEC [13] Table 2 Comparison of background counts evaluated with Gamma Vision and SNIP routines for peaks and multiplets is identified in the gamma-ray spectrum reported in Fig. 3 ROI Channel range m Background counts (Gamma vision) Background counts (SNIP-type) Difference (%) The analysis of multiplets in the spectrum reported in Fig. 3 allows to highlight the main differences with the results obtained through the use of the commercial software available today for the analysis of spectrometric data.As the most programs, the evaluation of the continuum trend is carried out in the peak fitting procedure.A comparison is realized using the fitting routines provided in Gamma Vision version 9 code furnished by ORTEC [13], the program used inside our laboratory.The continuum trend is obtained through a linear approximation considering the average contents of the first three and the last three channels of a ROI.The linear function estimates the continuum under the peaks as a trapezoid and is a simple, straightforward equation that is adequate when the continuum in the spectrum is relatively flat.To optimize the fitting parameters included continuum data, the "Interactive in Viewed Area" option has been used.The comparison between the values deduced with the Gamma Vision software and the ones of SNIP-based procedure is shown in Table 2.It can be seen that the differences are not significant for ROI related to a single photopeak since the linear continuum approximation may be sufficient (except for the peak of ROI #4).Significant differences occur in correspondence of multiplets or for peaks mounted on Compton edges of gammas with higher energy, as ROI #5, whereby the continuum values computed with Gamma Vision is significantly lower than the one evaluated with SNIP procedure as highlighted in the insert of Fig. 3. Other programs, for example, Genie 2000 [14], provides the use of a step function as an alternative to linear function, while a more modern program, HyperLab 2023 [15], can adopt a polynomial function or a step or tail function with respect to baseline trend.
Regardless of the type of function used, which make more or less complicate the peak fitting procedure, the advantage of using the proposed SNIP procedure is the easy and fast evaluation of baseline trend over the whole gamma-ray spectrum.The results are comparable with the ones of other codes in case of simple multiplets analysis, in particular for single peaks, but its use can be important in case of complex multiplets as highlighted in Fig. 3.
Figure 4 shows the comparison between the background trend and the one of a spectrum detected on an atmospheric particulate filter measured with the GLP planar detector particularly suitable for the measurement of low energy photons [16].The presence of 7 Be in the measurement of the air particulate filter modifies the MDA of 210 Pb at the energy of 46.5 keV by 38%, confirming the influence already highlighted in [17].
Fig. 4 Comparison baseline estimation for gamma-ray spectra detected on an atmospheric air particulate filter and on the corresponding "blank" sample.Live counting time for both the spectra: 80,000 s In Fig. 5, the corresponding trends of three spectra measured under the same conditions with a "low-background" ORTEC GEM50195S detector are reported [18].The first spectrum is that of the background of the system, the second is related to the measurement of a Marinelli beaker filled with water drawn from a public distribution system and the third of an atmospheric particulate filter measured after a short decay period.The same figure shows the baseline estimation trends, from which it is possible to evaluate an increase in MDA in the 137 Cs area of 4% for the Marinelli beaker and 8% for the filter compared to the background gamma-ray spectrum, and with reference to 122 keV energy the variations are of the order of 4.5% and 12%, respectively.
It should be also noted that the baseline trend in correspondence with the annihilation photopeak at 511 keV presents in all spectra a wider peak with an higher FWHM with respect to the calibration value, and therefore the value of m for routine analysis is always underestimated.To avoid the effect highlighted in Fig. 5 a special routine must be adopted to increase m value only for 511 keV annihilation peak region.
Finally, in Fig. 6 is reported a portion of the spectrum related to an atmospheric particulate filter taken on May 3, 1986 and measured with the ORTEC GEM18180 detector after a delay of 440 h from the end of suction.Many photopeaks related to the radionuclides present in the radioactive cloud following Chernobyl accident are highlighted [19].Up to April 1986, sampling of atmospheric particulate matter on a paper filter was carried out almost daily with a high-volume sampler (about 15,000 m 3 per day).The spectrometric measurements were mainly aimed at determining the air concentrations of 7 Be and 137 Cs.MDA values were computed taking into account the gamma-ray spectrum obtained with the same detector on a blank filter for a counting time of 80,000 s, and resulted 0.4 Bq and 0.053 Bq for 7 Be and 137 Cs, respectively, After Chernobyl accident, the baseline trend of a gamma-ray spectrum detected on a filter, as the one reported in Fig. 6, leads to a variation of MDA up to 235 Bq for 7 Be and about 21 Bq for 137 Cs with a counting time of 7900 s, corresponding to an air concentration of 19 mBq m −3 and 1.7 mBq m −3 for 7 Be and 137 Cs, respectively.In fact, while gradually decreasing the value of MDA, 7 Be was not recognized for a long time after detection of arrival of Chernobyl accident cloud due to its low activity on the filter, while for 137 Cs the activity on filter was several orders of magnitude higher than MDA until about the year 1990.
MDA variations also occurred for other radionuclides that have not been identified for a long time after the Chernobyl accident such as the decay products of the radioactive chains of 238 U and 232 Th besides 40 K.
Fig. 5 Comparison between baseline trends obtained with an ORTEC GEM 50195S detector, a "low-background" type detector, with cryostat in HJ configuration and low-background shielding [15].(1) measurement of the background of the spectrometric system; (2) Gamma-ray spectrum of a Marinelli beaker with water withdrawn from the public supply system; (3) particulate air filter measured after a short time from the end of suction.Live counting time for all spectra: 500,000s

Conclusions
The evaluation of the baseline trend is an important step in the analysis of a gamma-ray spectrum, both for the identification of the peaks and for a correct evaluation of the MDA.The SNIP-type iterative procedure has proved to be effective and easy to apply to experimental data of a gamma-ray spectrum of any complexity as well as allowing easy implementation in commercial software such as those already available for the analysis of gamma-ray spectra.
Particular attention must be paid to the selection criteria of the parameter which defines the different operating windows, correlating its variability with the FWHM calibration curve and/or with the number of components of a complex multiplet.The correct adoption of the analysis parameters, with the suitable operating procedures for selecting the peak regions and determining the parameter m, lead to reliable results throughout the whole energy range.
Although originally oriented toward optimizing the photopeak search procedure, by subtracting the background values thus determined from the experimental data, the application of the same procedure for the evaluation of the MDA for a given radionuclide allows to obtain in a short time an indication of the sensitivity of the spectrometric measurements with the use of simple algorithms.

Fig. 1
Fig. 1 Schematic representation and functional diagram of the iterative procedure and vectors used for the computations

Fig. 2
Fig. 2 Representation of the baseline trend change as a function of the variable m.Photopeak at 1332 keV of 60 Co.Unsmoothed baseline values

Fig. 6
Fig. 6 Portion the gamma-ray spectrum detected on the atmospheric particulate sample taken in Palermo on May 3, 1986 and measured on May 21, 1986.Live counting time: 7900 s