GRAVEL unfolding of gamma-ray spectra of CeBr 3 detectors and related uncertainties

,

A : Unfolded spectra are needed for the application of some spectrum analysis methods. Frequently, iterative deconvolution methods are used, especially if the spectra have a high resolution. However, the uncertainties on the obtained results are difficult to judge. In general, the user of an unfolding code gets no information on the influence of chosen parameters, especially the choice of the first estimate of the unfolded spectrum. General studies regarding the uncertainties of unfolding algorithms are difficult to perform, if not impossible, and therefore not available. In this work, uncertainties of the results on the iterative GRAVEL algorithm, which is based on SAND II, are closely investigated by using medium resolution spectra recorded by using a CeBr 3 detector. Relevant parameters are varied systematically to gain some information on the associated uncertainties on the basis of examples. The uncertainty budget of the calculation of total fluences and doses of unfolded spectra is compared with the uncertainty budget of the results obtained by applying the conversion method. The results allow, at least, the estimation of the magnitude of typical uncertainties.

Introduction
Every measured gamma-ray spectrum is influenced by the characteristics of the utilized detector. The energy response of a spectrometer depends on the detector material, its dimensions and on the properties of the casing of the instrument. The absorption of photons in the detector material is energy-dependent in a characteristic way. The detector size governs the photo peak to Compton ratio and the detection probability, and the casing will absorb especially low-energy photons depending on its material and construction. In many cases, the detector will not fully absorb the energy of an entering photon. Scattering inside the detector may lead to a partial detection of the photon energy according to the Compton formula. When electron/ positron pair production happens, the escape of 511 keV photons is possible, so that artificial escape peaks will appear in addition to photo peaks. Therefore, a measured spectrum is shaped by detector properties. If the undisturbed fluence spectrum, which would prevail at the location of the detector, is needed for further evaluations, the reconstruction of this spectrum by unfolding is a standard method [1].
After measured spectra have been unfolded, the obtained fluence spectra can be converted to dose spectra in terms of air kerma or ambient dose equivalent, * (10) [2]. The latter quantity is -1 -

Methodical procedure
Back in 1967, the SAND II algorithm was introduced by McElroy et al. [6]. It was published as an empirical enhancement of simpler approaches which are, however, based on the same idea. An estimated fluence spectrum is multiplied by the response matrix of the applied detector and the obtained spectrum, in the following called test spectrum, is compared with the measured spectrum. The SAND II code optimizes the estimated fluence spectrum iteratively and thus converges to a solution automatically (in contrast to earlier approaches). If the agreement is good, the unfolding was successful. A theoretical justification of the SAND II algorithm was published by Sekimoto and Yamamuro in 1982 [7]. However, this argumentation cannot fully explain the actual structure of the iterative SAND II formula. Nevertheless, this method has been used until today, because in many applications the results were satisfying and even better than those of other approaches. Especially, they turned out to be satisfactorily reasonable if the boundary conditions were observed. However, the SAND II algorithm does not take uncertainties into account. This disadvantage was removed by introducing the very similar GRAVEL algorithm [8]. In 1994, a theoretical rationale for this algorithm and SAND II, based on a maximum likelihood approach, was published by Matzke [9]. The idea behind is the 2 minimization of the difference between the measured and the reconstructed spectrum (using the detector response matrix) with logarithmic weighing. Afterwards, an experimental verification has been performed by different users. To mention just one example, in one case it turned out that the results of the more complex MAXED algorithm, which is based on a maximum entropy method, and those of GRAVEL were nearly identical [10]. Another verification of the GRAVEL algorithm on the basis of examples is found in [3]. The suitability of the GRAVEL algorithm is partly based on the fact that, like SAND II, it allows for the inclusion of a priori information. Nonetheless, a systematic prediction of the results of both codes is not possible, because differing results will be obtained depending on the input parameters. The quality and detailedness of the a priori information has an influence on the final results.

Unfolding
The goal of unfolding is the reconstruction of the undisturbed fluence spectrum which would prevail at the location of the detector if the detector was taken away. To remove detector effects from a measured spectrum, the photon response of the detector has to be known as precisely as possible. Mathematically, the algebraic equation ì = M · ì has to be solved after measuring the spectrum ì and determining the detector response matrix M, so that the vector ì , the undisturbed fluence spectrum, is found as a result.

Matrix inversion
If M is quadratic, the simplest algebraic approach would be ì = M −1 · ì . All matrix elements, which are calculated by Monte Carlo (MC) simulations, are linear independent. However, the calculation of the inverse matrix of a large matrix width e.g. with 1000 × 1000 elements is not possible. The reason is that the matrix elements have at least small statistical uncertainties as a result of the MC simulations. As a consequence, the inversion algorithm makes the inverse elements grow or shrink to extremely large or small positive and negative values (often only limited by the precision of the variables in the computer code). Hence, matrix inversion does not work and the original fluence spectrum has to be reconstructed in a different way, i.e. by unfolding.

Description of the GRAVEL algorithm
The iterative algorithm GRAVEL is selected because it can handle high-resolution spectra very well in acceptable computation time. The iterative GRAVEL formulas, originally derived by Matzke [9], can be rewritten in the form [3]: where ∈ {1, 2, 3, . . . , } and ∈ {1, 2, 3, . . . , }, In these formulas, the so-called test spectrum simplifies the representation of the equations by defining: In the following lines, the meaning of the variables is explained (element denotes the number of counts in one bin): : th element of unfolded spectrum of iteration : response matrix element 1 : 1 st estimate of an element of the unfolded spectrum : normalisation matrix element : th element of test spectrum : normalisation matrix element : Element of measured spectrum : iteration step (starting at 1) : Estimate of measurement error of bin -3 -

JINST 18 P07005
The GRAVEL formula can be modified in a simplified form as follows: where ∈ {1, 2, 3, . . . , } and ∈ {1, 2, 3, . . . , }, defining a recalculated normalization. If Poisson statistics is valid (i.e. 2 / 2 = ) one can write: (2.4) In the latter notation of the GRAVEL formulas, it is obvious that zero entries in the measured spectrum and in the test spectrum are not allowed, namely the code has to omit zero bins in both spectra. The GRAVEL algorithm iteratively adapts the test spectrum ì , corresponding to the unfolded spectrum ì times the response matrix R, to the measured spectrum ì . In the case of a perfect agreement, the unfolded spectrum is not changed, any more, by further iterations ( +1 = for all ). The basic idea behind is to maximize the likelihood of the probability ( ì | ì ), so that the expectation of the vector ì expresses the most likely solution. But this is not strict maximum likelihood estimation because the uncertainties and the choice of the first estimate can lead to slightly different solutions in a space of solutions.
The adaptation of the test spectrum to the measured spectrum is controlled by the normalization matrices W in eq. (2.2) and V in eq. (2.3), which are calculated from the response matrix R. This normalization is needed because the unfolded spectrum has to depend on the absolute number of counts in the bins of the measured spectrum, but it must not depend on the absolute size of the matrix elements. The normalization matrices include a normalization of the scalar product of the columns of the response matrix times the unfolded spectrum to unity (in the case of V, the scalar product of the numerator was not written explicitly, because the actually cancel out in the exponent of eq. (2.1)).
The formulas above describe the SAND II algorithm, if the terms in the rectangular bracket are left out. The fact that this difference is only small explains the similar behaviour of both algorithms. As a result of own investigations based on example measurements, the following tendency was found: if the uncertainty on the counts in the bins of a spectrum is governed by Poisson statistics, both algorithms will run into similar solutions, but GRAVEL converges slightly faster.
In case of Poisson statistics, the term in the rectangular brackets in eq. (2.1) and (2.3) are just equal to . This means that the normalization by applying is proportional to the product of the element of the measured spectrum (i.e. ) times the product of element of the unfolded spectrum (i.e. ) with the matrix column element . In eq. (2.2), this has the effect of a drastic amplification of the logarithm in the exponent of by a high normalization factor when the counts of the measured and unfolded spectrum are high, so that the unfolded spectrum is strongly modified -4 -in high count regions of the spectrum. The entries of the measured and of the estimated spectrum as well as all response matrix elements have to be positive, otherwise unphysical solutions will be obtained. As a consequence, strategies are needed to smooth or correct measured difference spectra (foreground minus background spectra) in a way that the generation of artefacts is excluded.
The selection of a reasonable first estimate of the iterative process is important. If the energy resolution of the applied detector is high enough (as is true for examples used in this work), either the measured (after applying a smooth spline to reduce statistical variations) or an artificial spectrum can be used. If there is good knowledge about all peak positions, the latter can consist of a flat background with superimposed peaks featuring a width adapted to the detector resolution. The latter method can only be realized if there is good knowledge about all peaks and their positions. Otherwise, unphysical artefacts will be produced. In this work, where the FWHM (full with half maximum) value of the photo peaks is in the order of a few percent, a ratio of the peak area to the background level of 1000:1 is at least needed to have a visible impact on the unfolding (in comparison with a flat first estimate).
Because a perfect agreement of the measured spectrum and the test spectrum will never be reached in reality, a criterion is needed when to stop the iterative process and avoid infinite oscillations of the unfolded spectrum. As a helpful indicator, the parameter is defined as [3]: where is the number of degrees of freedom, i.e. the number of spectrum bins. If this parameter falls below a predefined value or passes a minimum, the iterations are to be stopped. Unphysical oscillations of the unfolded spectrum may lead to a rise of the parameter after several iterations. A second observable is the relative change in the integral over the unfolded spectrum, which is defined as: If the test spectrum is gradually adapted to the measured spectrum, the parameter Δ becomes small. The parameters and Δ are independent from each other, so that only one of them may undergo a minimum (or none).
Because the elements of the response matrix are related to (implicit) energy scales of the vectors in eq. (2.1), the use of a quadratic matrix is highly favourable, with the consequence that only one energy scale exists. The high channel number implementation of GRAVEL in the UMG package only works with quadratic matrices [8], whereas the implementation in this work can also handle non-quadratic matrices. However, a general requirement is that the number of bins of the unfolded spectrum (especially the first estimate) is identical with the number of columns, , and the number of bins of the measured spectrum is identical with the number of rows of the response matrix, . In the case < , the unfolded spectrum has a higher granularity than the measured spectrum. A typical example is the evaluation of the data of a few-channel spectrometer. In this case, the measured data are not sufficient to explain the complete structure of the unfolded spectrum. Hence, the use of the measured spectrum as a first estimate of the unfolded spectrum is not possible. Instead, it is inevitable to introduce additional information when constructing the first estimate. This will have a strong influence on the final result.
-5 -In principle, unfolding with a matrix with more rows than columns, > , is possible. In such a case the question has to be raised, how the energy scales of the measured spectrum and the estimated spectrum match, as one and the same element of the matrix is related to different scales in horizontal and in the vertical direction. A possible application would only be = · , where is an integer number. This would lead to a case where one bin of the estimate is related to several bins of the measured spectrum. This means that the response matrix would spread one bin of the unfolded spectrum to several bins of the measured spectrum, introducing an additional uncertainty to the solution. Larger problems can be expected if the multiplying factor was not an integer. Now the rounding of the mean bin energies to the two energy scales would lead to artefacts and a higher uncertainty. Therefore, in this work only quadratic response matrices are calculated and used.

Implementation of the unfolding procedure
In order to have a sufficient energy resolution at low energies, 4k spectra are always used in this work. Therefore, detector response matrices with a size of 4k × 4k are to be calculated by Monte Carlo (MC) simulations. The detector geometries were implemented in the MC code MCNP 6.2. Matrices were calculated for a simple geometry, which can be generalised based on the following scenarios. A detector is positioned in a distance of 1 m from the source, so that the divergence of the incoming beam can be neglected. An energy range up to 10 MeV is covered for a small CeBr 3 detector (1 × 1 ) and a large CeBr 3 detector (1.5 × 1.5 ). Dealing with 4k spectra, 4096 simulations have to be performed to obtain one response matrix. Each MC spectrum is folded with the energy-dependent detector resolution, so that the simulated spectra resemble measured spectra [3]. Folding with the detector resolution also smoothes the ripples of the MC spectrum due to statistical uncertainties. Each obtained spectrum represents one column of the response matrix.
The GRAVEL algorithm is implemented in the homemade SpecConvert Windows code [3], capable of handling non-quadratic matrices. This code is also used to convert spectra formats, to stretch or shrink measured spectra to the energy scale of the response matrix, and to calculate smooth splines of measured spectra to generate first estimates for unfolding. The code can also invert matrices and unfold spectra using the SAND II algorithm. A well-known implementation of GRAVEL is found in the UMG software package, which is distributed for free by PTB [8]. However, the FORTRAN code of this package was not used because it cannot read standard file formats and because it cannot cope with 4k×4k response matrices and 4k spectrum vectors due to limited memory capacities. In addition, the high channel number implementation of GRAVEL in the UMG package only works with quadratic matrices [8].

Investigation of the GRAVEL algorithm
First of all, the GRAVEL algorithm will be applied to measured and simulated spectra to investigate the behaviour of this algorithm and its typical results. This will be done using spectra with different complexity. In a second step, systematic investigations will be performed to work out typical uncertainties on the results of the method. This is also done using examples, because a general theoretical investigation of the propagation of errors in the code is not possible. The GRAVEL algorithm itself does not deliver uncertainties on unfolded spectra. All simulations and measurements were performed with a 1.5 ×1.5 CeBr 3 detector in a distance of 39 cm from the source with the exception of the spectra related to the R-F field, where the same detector was placed 1 m in front of the target.

Investigation of Cs-137 spectra
The simplest case of a spectrum is represented by a (quasi) mono-energetic spectrum, here a spectrum of a 137 Cs source, which only includes a photo peak at 662 keV and the related Compton continuum at lower energies as well as a photo peak at 32 keV (actually two unresolved X-ray lines slightly below and above this energy). In figure 1, the unfolding of a measured 137 Cs spectrum is visualized. The measured spectrum is used as a first estimate. After one unfolding cycle, the pattern of the Compton continuum is already drastically reduced. When the number of iterative loops increases, this structure further decreases, while the emission peaks are increased. The peak at 662 keV is heightened, because the detection efficiency of the detector at this energy is about 0.5. By unfolding, the loss of counts by incomplete energy deposition is compensated for. Here and in the following, the -axes use arbitrary units, whose absolute values are not of interest, unless otherwise indicated.   Some bins above the 662 keV peak have non-zero entries caused by statistical uncertainties, because the unfolded spectrum is based on a difference spectrum (after background subtraction) and negative entries had to be removed before unfolding. However, as a basic characteristic of GRAVEL,  structures with counts lower than 1% of those of the peaks cannot be regarded as reliable, especially if neighbouring bins have zero entries. An uncertainty is caused by the fact that "over-unfolding" is possible, i.e. the iterative process is not stopped before artefacts are produced. After more than 10 loops, the unfolded spectrum starts to oscillate more and more due to the misinterpretation of statistical uncertainties in the measured spectrum. A similar outcome was found in all other spectra that were investigated in this work. This is a general behaviour that can be explained by the fact that the normalisation of the applied response matrix leads to a preferred alteration of high-count bins. This even leads to a periodic self-amplification of statistical fluctuations. Another general characteristic of GRAVEL unfolding can be seen in figure 1. The FWHM of the unfolded peaks is slightly smaller than that of the measured peaks ( figure 1, top), whereas the FWHM of the test spectrum peaks is slightly larger ( figure 1, bottom). This is also a consequence of the construction of the GRAVEL algorithm. If a spline of the measured spectrum is used as an initial estimate, the unfolding algorithm converges to the final solution slightly faster.
The unfolding of a simulated 137 Cs spectrum is presented in figure 2. The simulation was done by using MCNP, also including air as a surrounding material of the detector, and folding the result -8 -

JINST 18 P07005
with the detector resolution. The statistical uncertainty of the simulation leads to results similar to those shown in figure 1. The apparent peak at about 200 keV is not completely removed by the unfolding because this structure is produced by photons which underwent Compton scattering outside the detector and deposited electrons in it. Such detailed geometrical effects are not included in the response matrix because the response matrix is intended for general use, preferable for the evaluation of data recorded in a large room, where there are no walls in the close vicinity of the detector.
To demonstrate the dependence of the final result of unfolding on the input spectrum, an artificial spectrum is constructed as a first estimate. This spectrum includes both known emission peaks of 137 Cs and has a background of 1 in all other bins, so that the unfolding algorithm can adapt all bins in the region of interest (up to 750 keV). The area under the photo peaks in the first estimate is identical, so that the adaptation of the final peak height of the unfolded spectrum is left to the unfolding procedure. Again, the simulated spectrum depicted in figure 2 is unfolded, so that the results (figure 3) can be compared directly. Under the ideal conditions on which figure 3 is based, already the first iteration leads to a good result, in which the Compton continuum is drastically reduced to even less than a per mille of the 662 keV photo peak. Further iterations almost do not

JINST 18 P07005
alter the picture. Only the structure at 200 keV caused by Compton scattering outside the detector is maintained, which is physically reasonable. Due to the fact that the response matrix includes the actual detector response, the 662 keV peak is elevated much more than the low-energy X-ray peak, because the peak detection efficiency is much lower at higher energies. In a further attempt, a constant function is used as a first estimate of the unfolded spectrum ( figure 4). After a number of iterations, results very similar to those displayed in figure 2 are reached. The main difference is that the GRAVEL algorithm converges slower to the final result. The reason is that the peaks in the simulated spectrum are so pronounced that GRAVEL (mainly looking for high bin entries) is forced to bring theses peaks up. However, in a situation where only small peaks are present, such a simple approach will not lead to a satisfying result, as will be shown below. To trigger the appearance of clear peaks in the unfolded spectrum, either the peak-to-background ratio of the measured spectrum or of the first estimate (or of both) must be high. Therefore, care has to be taken if peaks are defined in the first estimate. Otherwise, artefacts could be produced.

Investigation of Bi-207 spectra
The reconstruction of fluences using the GRAVEL algorithm in a more complex case can be studied on the basis of an evaluation of a 207 Bi spectrum. The measured spectrum exhibits four prominent peaks at 74 keV (two unresolved X-ray lines), 570 keV, 1064 keV and 1770 keV, whose heights differ by orders of magnitude (figure 5, dark blue curve).  In a first attempt, the measured 207 Bi spectrum is unfolded by using a spline of the measured spectrum as a first estimate of the unfolded spectrum ( figure 5). In the result, the Compton continuum of scattered photons is drastically reduced in the region of high counts (below 500 keV), but it is less reduced where the photo peaks are less pronounced. Such behaviour is expected because GRAVEL will not significantly alter bins with a low number of counts.
In figure 6, the unfolding of the same spectrum, using a) an artificial spectrum, which includes predictions of the four prominent peaks (peak-to-background ratio 100:1 -see insert), and b) the measured spectrum as a first estimate are compared. The results look very similar, if the obtained spectra of case a) (after 15 iterations) and case b) (after 20 iterations) are compared. The background consisting of low-count bins in the right part of the spectrum is even less reduced, if the first estimate is artificially composed. The rise of the unfolded spectrum (in orange) at high energies is an artefact produced by the cutting of the unfolding at 2.5 MeV. In general, a truncation of the spectrum close to a peak-like structure may lead to artefacts. This effect, however, is exaggerated by the half-logarithmic plot. It hardly has an influence on the calculation of the total fluence and dose. Conversely, the background of the unfolded spectrum is almost completely removed (not depicted), if the artificially constructed spectrum includes higher peaks (peak-to-background ratio 1000000:1). This observation verifies the expected behaviour. Integrated fluence and dose values differ only slightly (by less than 5%). The peak heights and derived total fluence and dose values hardly depend on the first estimate. The unfolding reveals physical facts. This is a great advantage, if spectra of an unknown source have to be evaluated.

Investigation of spectra recorded in the ISO R-F photon field
A spectrum recorded in an R-F photon field has a higher complexity and covers a larger energy range than the spectra investigated before (figure 7). In this context, the question of "unknown" or "hidden" peaks in the spectrum is addressed. The example above shows that unfolding is especially powerful Reconstructed spectrum n=10 Figure 7. Unfolding of a measured R-F spectrum using the measured spectrum as a first estimate of the unfolded spectrum.

JINST 18 P07005
when a physically reasonable first estimate is used as an input for the GRAVEL algorithm. This leads to the question whether peaks may be suppressed when an unfortunate or incomplete estimate is chosen.
First of all, a measured R-F spectrum is unfolded by using a smooth spline of the measured spectrum (figure 7, dark blue curve). In the unfolded spectrum the physical peaks at low energies (<2000 keV) are raised, but at high energies the unfolding does not remove the escape peaks completely. Especially the escape peaks at 5.1 MeV, 5.6 MeV and 6.6 MeV remain clearly visible. In this energy region, the escape peaks are larger than the photo peaks. The peak at 6.1 MeV is ambiguous, because the photo peak at this energy is superimposed by the second escape peak of the 7.1 MeV peak.
To solve this issue, the low-energy structures associated with high-energy bins, here escape peaks, had to be removed from high energies to low energies, which is the basic idea of the stripping method. However, the GRAVEL algorithm adapts the whole estimate of the unfolded spectrum at the same time so that the measured spectrum is reproduced ( figure 7, green curve). Especially, the very pronounced escape peak at 6.6 MeV is not removed, because the test spectrum can be well adapted if this peak is "(mis)interpreted" as a photo peak. However, mathematically, the parameter stays high, so that a problem is indicated. Unfolding cannot separate escape peaks from photo peaks, unless a sound first estimate is provided. Thus, in the next example the first estimate of the unfolded function is constructed artificially so that only peaks are stressed, which are expected as a consequence of the nuclear excitations or reactions in the target (figure 8). The artificially constructed first estimate includes the expected peaks, all with the same area (10000), and a background set to 1 in every other bin (figure 8, insert). Under these conditions, unfolded spectra do not show the escape peaks any more, while the physically expected peaks are strongly enhanced. Interesting is the structure which looks like a very wide, flat peak between 1.5 MeV and 4 MeV. This corresponds to the expected distribution of the energy of electron/ positron pairs that are emitted by an excided state of 16 O with a maximum energy of 4.5 MeV, peaking at about 3 MeV. This is a hint that unfolded spectra still can reveal additional physics, which was included neither in the detector response matrix nor in the first estimate of the unfolded spectrum. In figure 8 the

JINST 18 P07005
three high-energy peaks have almost identical peak height. This is not in agreement with literature data [13] and is a consequence of the first estimate of the unfolded function. This unwanted effect will be further investigated below.

Incomplete first estimate of the unfolded spectrum
Using a physically well justified first estimate of the unfolded spectrum is helpful to obtain a physically meaningful result by unfolding. This is especially the case if non-diagonal elements in the response matrix are not small. On the other hand, the physical knowledge about any measurement may be limited. What is the consequence of an incomplete estimate of an unfolded spectrum? This question cannot be addressed in general. One example will be studied, instead.
A small 137 Cs peak (at 662 keV) is artificially added to the measured spectrum recorded in the R-F field (last section). This peak is so small that its top does not even double the counts of the affected bins of the continuous background. The net area below this peak is much smaller than that of all other peaks. The modified R-F spectrum is unfolded again, as described in the section before. Both approaches, by using (a spline of) the measured spectrum as a first estimate (figure 9) and using an artificially constructed spectrum as a first estimate (figure 10) lead to unfolded spectra where the 137 Cs peak is clearly visible, though the artificial first estimate does not include this peak. The peak area under the 662 keV peak is nearly identical in both cases. One can conclude that unidentified peaks will not necessarily be removed by unfolding, even if they are small. This is at least true for a spectral region where no pronounced escape peaks are involved, so that the response matrix is highly diagonal. Nonetheless, a calculation of a parameter like a detection limit is not possible, because the shape of the unfolded spectrum also depends on the whole spectrum and on the absolute number of counts in any region of interest. It has already been pointed out that structures in the percent range (related to the spectrum bins with the highest counts) are not unfolded according to physics, but mainly to the first estimate. Unfolding cannot reveal new information in bins with a low number of counts.

Systematic investigation of uncertainties
Because a general, systematic investigation of the propagation of uncertainties by the GRAVEL algorithm is not possible, a detailed study of the consequences of a systematic variation of input parameters is presented in this section based on examples. Assuming typical scenarios, this will show how sensitively calculated fluence and dose results will depend on input parameters.

Effects of supposed detector resolution and first estimate choice
The energy resolution of a specific spectrum may not perfectly agree with the supposed detector resolution implemented in the response matrix. Apart from the pure detector characteristic, the actual detector resolution also depends on the electronics including the amplifier and the digitizer as well as on the temperature stability. The latter will especially have an influence on the gain of photomultipliers and analogue main amplifiers. To investigate the consequence of a non-ideal energy resolution, the energy resolution of the detector response matrix of a CeBr 3 detector is varied relatively to the energy resolution of a spectrum, i.e. the simulated 137 Cs spectrum (figure 2). For this test a simulated spectrum is chosen, because it has exactly the same energy resolution as that used for the calculation of the response matrix, so that a precise relative change in the energy resolution can be modelled. Furthermore, the implication of the use of different first estimates of the unfolded spectrum is investigated. Two examples of unfolded spectra are plotted in figure 11. The peaks at 32 keV and 662 keV are similar. But when the first estimate is identical with the spectrum to be unfolded, the region of the Compton scattered photons is not properly removed. If an artificially composed curve is chosen as a first estimate, the background of the unfolded curve oscillates with a wave length identical to the detector resolution. Only physically relevant cases are chosen (table 1). In all cases, the iterative process is stopped when a most similar of about 20 is obtained. When the energy resolution of the matrix is too high, a in the order of 20 is never reached. Concerning the fluence (table 1) and the dose (table 2) in terms of * (10), the differences resulting from different unfolding input parameters are listed in the right two columns, respectively.   Table 3 lists fluence data, table 4 dose data in terms of * (10). The total fluence and the fluence of the high-energy peaks are hardly influenced by the input parameters, if the iterations are stopped at a loop where the parameter is similar. The same is true for the dose results. An exception is the fluence and the dose in the high-energy region, when the peaks in the first estimate of the unfolded spectrum are strongly emphasised (last lines of table 3 and table 4), where deviations of 9% and 8% are observed, respectively. Table 3. Fluence data derived from a measured R-F spectrum. The fluence data of the total spectrum, of the high-energy region ( > 5 MeV) and of the peak at 6.1 MeV are listed in the arbitrary unit cts (counts). The relative differences, normalized to the 2 line are listed in the last 3 columns. Legend: matrix: nw: normal (regular) width; 1 st estimate: calc: calculated from the simulated spectrum, s pks: small peak prediction (peak-to-background ratio 200:1), sm pks: small/medium peak prediction (peak-to-background ratio 1000:1), m pks medium peak prediction (peak-to-background ratio 10000:1); n: number of unfolding cycles; : according to eq. (2.5).  -17 -On the contrary, the peak area of the 6.1 MeV peak strongly depends on the input spectrum. This means that peak areas may diverge by an order of magnitude after unfolding. Hence, a quantative analysis of this peak by unfolding is not possible. It was already pointed out that this peak includes the sum of two peaks, the double escape peak of the 7.1 MeV emission peak and a further emission peak at 6.1 MeV. The GRAVEL algorithm does not include physical information so that it cannot separate the two problems, though the height of escape peaks is included in the detector response matrix. The overall best adoption of the test spectrum does not guarantee the absence of local problems. A poor reproduction of the peaks, like that shown in figure 7, indicates that unfolded peak areas are hardly reliable. Global parameters like and the total integral do not reveal the fact that the measured spectrum has been locally incompletely unfolded. In this spectrum, the area of 6.1 MeV peak represents a worst case scenario.

Uncertainty budgets
Using information of the preceding section and additional input from similar investigations, uncertainty budgets can be calculated. In table 5 and table 6 the budgets for the uncertainty calculation concerning unfolding (using GRAVEL) and the conversion method are estimated (see [3] for detailed information on the conversion method). The maximum deviation due to the detector resolution and the 1 st estimate is taken from tables 1 to 4. The results of the same MC simulations could be used to identify the uncertainty due to the number of iterations, as the output spectra of each iteration step were recorded. Concerning the influence of an energy shift, the energy calibration of spectra to be unfolded was altered in a span which may occur in realistic scenarios. Overall calibration factors were calculated in [3], so that the corresponding typical uncertainties are known. The uncertainty due to the MC physics model, which e.g. depends on implemented cross-sections and material constants, is taken from comparisons of different codes in which the same problem is implemented [3]. For the convenient use of discrete data sets in the applied spectrum to fluence/ dose conversion software, analytical functions were fitted instead of using simple interpolations. The precision of the fits is limited and the maximum deviation of single calculated values from pre-defined discrete values was investigated (uncertainty of non-MC effects). The use of interpolations (instead of fitting) would not improve the uncertainty budgets.  Total uncertainties are calculated by using a method proposed by the Guide to the Expression of Uncertainty in Measurement (GUM) [14], valid for a coverage factor of = 1. The influencing parameters are regarded as independent from each other so that all components of the uncertainty can be added quadratically. The data in table 5 are related to fluence, whereas those in table 6 are related to * (10). The results are, nevertheless, very similar.
-19 -In summary, the evaluation of unfolded spectra implies a larger total uncertainty because this method suffers from some sources of uncertainty which do not affect the conversion method [3] as the detector resolution, the choices of the first estimate of the unfolded spectrum and the number of iterations (the conversion method is not iterative and does not need the definition of a first estimate). Hence, the data derived from unfolding have an uncertainty of about 5%, which is almost twice the one obtainable with conversion (2.8% in both cases).
When the uncertainties in the total fluence and * (10) of the wide and complex R-F spectrum are calculated, again values of about 5% are obtained, which is very similar to the previous results. However, the calculation of the peak area of the 6.1 MeV peak is very uncertain, because it strongly depends on the first estimate of the unfolded spectrum, as illustrated by data in table 3 and table 4. It is very important to note that the unfolding of a spectrum which comprises two overlaying peaks, including a pronounced escape peak, is problematic. If no escape peaks are present and the response matrix is mainly diagonal in the relevant region, even the investigation of a single peak by unfolding is possible with a reasonable uncertainty. This can be concluded from data in table 1 and table 2.

Inexpedient approaches
In addition to established methods, also unfolding by matrix inversion is discussed in textbooks, articles and on internet platforms. There are also different ways how to perform iterative unfolding. Below, the non-suitability of two approaches will be demonstrated by examples.

Matrix inversion
In related literature, often the advice is given not to perform matrix inversion if the matrix in question is large. Already a matrix with 20×20 elements can be regarded as large, because an analytical inversion formula for the single matrix elements does not exist if the matrix has × elements with > 3. Nevertheless, some authors claim to have solved a problem by matrix inversion, e.g. [15]. Therefore, a matrix inversion of a 4k × 4k response matrix of a CeBr 3 detector, which was used above for unfolding, has been inverted to demonstrate the outcome by using the MATINV method published by Bevington [16].
All columns of the inverted matrix look similar except for the fact that the absolute scales are different. One randomly chosen example of a matrix column of the inverted matrix is shown in figure 12. Extreme and unphysical oscillations, even between positive and negative values, are present. The half-logarithmic plot, made after removing the negative entries, shows an underlying curve, studded with extreme spikes. This result is so unreasonable that multiplying the measured spectrum with the inverted matrix would not make sense. This unwanted behaviour is due to very small statistical uncertainties in the response matrix that are amplified during computation. Additional mechanisms would be needed to smooth the matrix, at the risk of introducing other flaws. This example illustrates why matrix inversion should be avoided, especially if the matrix in question has many elements.

Neglecting the detector resolution
In a paper where a similar (but simplified) iterative unfolding algorithm was used, the author reported that quite a reasonable unfolding of a simulated spectrum was successful when the unfolding matrix did not include the detector resolution. Including the detector resolution just led to a slightly better peak resolution [17]. Under such circumstances, unfolding the spectrum that is shown in figure 1 -20 -2023 JINST 18 P07005  leads to results like those displayed in figure 13. The photo peak at 662 keV is properly reconstructed. However, a tendency is visible that spikes, which are actually artefacts, are produced by noise. The regions with a low number of count oscillate heavily. The unfolded spectrum is not improved because small statistical uncertainties, which are included in the measured spectrum and in the response matrix, are amplified in an unphysical way.  Figure 13. The same as figure 1, top, but using a response matrix which does not take the detector resolution into account. The gaps are due to the spectrum resolution of 4k bins. The -axis represents counts in arbitrary units.
If neither the simulated spectrum nor the response matrix are folded with the detector resolution, unfolding is possible, nevertheless. In figure 14, ripples are produced by the fact that the calculation of the response matrix with MCNP has some uncertainties. When the input matrix is broadened using the detector resolution this effect is hardly visible because the broadening serves like a smooth  Figure 14. The same as figure 1, top, but using a simulated spectrum and a response matrix, which both do not take the detector resolution into account. The -axis represents counts in arbitrary units.
spline. In regions of the spectrum where the number of counts is low, the statistical fluctuations are even amplified to peak-like structures. They are pure artefacts. This demonstrates that leaving out the detector resolution is generally not fortunate. In summary, the detector resolution should always be taken into account when realistic high-resolution spectra (with statistical ripples) are unfolded. Few-channel spectra (not discussed in this paper) include an intrinsic averaging (by the rough binning), so that taking the detector resolution into account is not mandatory.

Summary and conclusion
Everyone who uses unfolding methods for the evaluation of spectra has to face the question how reliable the obtained results are, especially if the used algorithm does not supply uncertainties. The GRAVEL algorithm e.g. does not deliver information on uncertainties. In general, it is unclear how uncertainties in the input data propagate into the final result. Hence, several examples of different complexity were studied.
The energy resolution of the used detector is of high importance when spectra have to be unfolded. In this work, spectra measured with a CeBr 3 detector were studied. They have a medium resolution, between that of typical NaI spectra and HPGe spectra. The CeBr 3 detector can be used to record spectra in a wide energy range. Unfolding matrices have been calculated up to an energy of 10 MeV. The energy resolution of this detector is high enough to reveal a number of photopeaks if relevant emission lines are included in a spectrum [3,18]. This allows unfolding of measured spectra with the GRAVEL algorithm just by using (a spline of) the measured spectrum as a first estimate of the unfolded spectrum, so that no prior knowledge of the expected unfolded spectrum is needed. Hence, biases in the results due to imperfect knowledge of the real gamma field are avoided.
However, especially if the measured spectrum includes pronounced escape peaks, the use of an artificially constructed first estimate may lead to physically more plausible results. The accuracy of an unfolded spectrum can be judged by a comparison of the measured spectrum and the algebraic product of the response matrix with the unfolded spectrum. The overall agreement of both spectra is quantified -22 -by the parameter, which may serve as an indicator when the iterative unfolding procedure has to be stopped. Peaks visible in the measured spectrum can be used to artificially compose a first estimate of the unfolded spectrum. Thus, no relevant peak will be overlooked. Even if one peak in the measured spectrum is not taken into account in the first estimate, this peak may be nevertheless enhanced in the unfolded spectrum, as was demonstrated by means of an example. A condition is that the bins of this peak have considerably more entries than the background. In addition, the absolute peak height has to be higher than several percent of that of the highest peaks in the spectrum, because GRAVEL will not unfold small background structures. Instead, GRAVEL amplifies bins of the measured spectrum with high counts.
When the total fluence or dose is to be calculated from a measured spectrum, the conversion method is in general superior to an evaluation of the fluence via unfolding (and a subsequent dose calculation), because unfolding brings about more sources of uncertainty. This was shown on the basis of some examples. When, however, the reconstruction of the virtual undisturbed fluence at the detector position is necessary, unfolding is inevitable. The reconstruction of peak heights by unfolding is problematic if peaks overlap in a spectrum or in case a photo peak is small in comparison with other related structures, e.g. escape peaks. For unfolding low resolution spectra, like those of liquid scintillation detectors, a good knowledge about the structure of the anticipated unfolded spectrum is inevitable. Otherwise, reliable results will not be obtained. Especially when spectra with a medium or high resolution have to be unfolded, the GRAVEL algorithm is a powerful tool, which already may reach its best results after a few iterations in the order of 10. Because a new response matrix cannot be calculated for every unfolding problem, especially when the number of bins is high, additional tools are needed to adapt the energy scale of a measured spectrum to that of the response matrix and to handle negative bin entries that may result from the calculation of net spectra.
The uncertainty budgets of several examples were calculated to obtain some information on the precision that can typically be expected. The precision of total fluence and dose data obtained by the conversion method is higher than the precision of the same data obtained by unfolding (using GRAVEL). Regarding the examples, total fluences and doses were derived from spectra with a total uncertainty of approximately 3% (conversion) and 5% (unfolding) for = 1. This corresponds to the observation that the agreement of results obtained by applying the conversion method and by unfolding did not differ by more than 10% in practical examples, at most [3].