Nuclear Instruments and Methods in Physics Research A 689 (2012) 35–39
Contents lists available at SciVerse ScienceDirect
Nuclear Instruments and Methods in Physics Research A journal homepage: www.elsevier.com/locate/nima
Comparison between standard unfolding and Bayesian methods in Bonner spheres neutron spectrometry G. Medkour Ishak-Boushaki n, M. Allab ´ de Physique, Universite´ des Sciences et de la Technologie Houari Boume´die ne, BP 32 El-Alia BabEzzouar, Algiers, Algeria Laboratoire SNIRM—Faculte
a r t i c l e i n f o
abstract
Article history: Received 10 October 2011 Received in revised form 30 May 2012 Accepted 12 June 2012 Available online 23 June 2012
This paper compares the use of both standard unfolding and Bayesian methods to analyze data extracted from neutron spectrometric measurements with a view to deriving some integral quantities characterizing a neutron field. We consider, as an example, the determination of the total neutron fluence and dose in the vicinity of an Am–Be source from Bonner spheres measurements. It is shown that the Bayesian analysis provides a rigorous estimation of these quantities and their correlated uncertainties and overcomes difficulties encountered in the standard unfolding methods. & 2012 Elsevier B.V. All rights reserved.
Keywords: Neutron spectrometry Bayesian analysis Unfolding Bonner spheres spectrometer
1. Introduction Neutron spectrometry measurements are often carried out to determine some integral quantities from neutron energy spectra such as the total neutron fluence and the ambient dose equivalent in the radiation protection field. Generally, the measuring apparatus does not provide a direct measurement of the quantity of interest. In the standard approach, an unfolding procedure is used to derive the neutron spectrum from the data collected from the spectrometer and the quantities of interest are computed from the unfolded spectrum. However, the solution spectrum is not unique [1] inducing some difficulties to quantify the uncertainties associated with the derived quantities. In the Bayesian analysis [2], the quantity of interest has to be inferred from the existing information on the experimental process modeled in terms of probability distributions and the predictive probability distributions of the requested physical quantities are deduced by applying the probability rules attached to Bayes’ theorem. This approach should provide rigorous uncertainties for some integral quantities derived from neutron spectrometric measurements. As an example of application of both methods, we consider the problem of estimating some integral quantities characterizing the neutron field, in the vicinity of an Am–Be source of low activity ( 37 GBq), from a set of measurements carried out using a Bonner sphere spectrometer (BSS) at the Nuclear Research Center of Algiers (CRNA) [3,4].
n
Corresponding author. Tel.: þ213 21247344; fax: þ213 21247344 E-mail address:
[email protected] (G. Medkour Ishak-Boushaki).
0168-9002/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.nima.2012.06.016
The BSS is composed of a set of six polyethylene spheres of different diameters (0, 3, 5, 6, 7 and 12 in.) and a 6LiI(Eu) scintillator located at the center of each sphere. The fluence response, expressed in cm2, has been determined in Refs. [3,4] for each of the six Bonner spheres exposed at a distance of 55 cm from the Am–Be source (cf. Table 1). This response is defined, in these references, as the number of counts measured with each sphere normalized by the neutron fluence of an Am–Be standard source. Each sphere fluence response has been estimated with a relative experimental uncertainty equal to 5.2%. The number of counts measured with each Bonner sphere is analyzed by using either the standard unfolding procedure [5] or Bayesian’s method to extract information on the energy spectrum of the neutrons incident on the spectrometer. The total neutron fluence and the associated ambient equivalent dose characterizing the Am–Be neutron field are then derived and compared.
2. The standard unfolding methods The number of counts Nk measured in each sphere k (k¼1yn, where n is the number of Bonner spheres in the spectrometer) is related to the differential energy spectrum fE(E) by the linear equation Z Nk þek ¼ Rk ðEÞfE ðEÞdE ð1Þ where Rk(E) is the response of the sphere k to neutrons of energy E (the so-called ‘‘response-function’’) and ek is a term which accounts for the discrepancies between the measured and the calculated counts due to deviations of Rk(E) from the true value of the sphere
36
G. Medkour Ishak-Boushaki, M. Allab / Nuclear Instruments and Methods in Physics Research A 689 (2012) 35–39
4
Table 1 Experimental fluence responses (cm2) of the CRNA Bonner sphere spectrometer and their experimental relative uncertainties using 10 2 mm2 6LiI scintillator (taken from Ref. [4]).
ISO spectrum Unfolded with Gravel Unfolded with Maxed
Exp. Fluence Response 7 Exp. uncert.(%)
Bare 3 5 6 7 12
0.022 7 5.2 0.086 7 5.2 0.249 7 5.2 0.314 7 5.2 0.324 7 5.2 0.228 7 5.2
Fluence (cm-2)
3 Sphere diameter (in.)
2
1
0 response. The value of ek is not known a priori, but it is expected to be of the same order of magnitude as the estimated uncertainty assigned to the value Nk. It should be emphasized that the differential energy spectrum fE(E) is not measured directly; it has to be extracted from Eq. (1) by using an unfolding procedure. The mathematical solution of Eq. (1) is not unique because of the limited number of counts measurements (6 measurements in this study) available to estimate the quasi-continuous function describing the neutron spectrum fE(E). To get an appropriate physical solution, one need to introduce some physical additional assumptions that complement the measurements information and very often an ‘‘a priori’’ knowledge about the measured spectrum is taken into account. Several unfolding methods [1,6], based on different mathematical principles (e.g, methods of maximum entropy, least square methods and iterative methods), may lead to solutions fitting the measured data. To unfold Eq. (1), it is convenient to consider for computational purposes its discrete version where the whole neutron energy range is splitted in m energy bins Nk
X
R f l kl l
ð2Þ
where Rkl are the response matrix elements and fl is the l component of the total neutron fluence ftot, (l¼1ym). The problem of unfolding consists in determining the fluence components fl using Eq. (2), namely the data Nk, the response matrix Rkl and any other available additional information relevant to the problem [6]. Many methods are proposed for the unfolding of data from neutron spectrometers [1]. In this paper, to resolve Eq. (2), we have focused on two approaches largely used in the literature: the maximum entropy unfolding method developed in the Maxed code [8] and the least square unfolding procedure developed in the Gravel code [9]. In the calculations, the reference spectrum for the Am–Be source recommended by ISO-8529 [10,11], reported in Fig. 1, has been introduced as a default spectrum to initiate both unfolding procedures and m¼50 energy bins have been considered. Fig. 1 shows the unfolded spectra from the CRNA BSS measurements by means of both selected unfolding techniques cited above. The best description of the unfolded data is obtained by using the chi-squared method. In each case the chi-2 value was nearly equal to 0.9 suggesting therefore the equivalence of the selected unfolding methods. From the obtained neutron spectra some physical integral parameters characterizing the neutron field have been deduced. These are the total neutron fluence ftot and the ambient equivalent dose H*(10) written as follows:
ftot ¼
Z
fE ðEÞdE
X l
fl
ð3Þ
0
2
4
8
6
10
12
Neutron energy (MeV) Fig. 1. Am–Be neutron spectrum unfolded using two different methods (Gravel and Maxed) compared to the typical standard spectrum (ISO spectrum).
Table 2 Total neutron fluence (ftot) and equivalent ambient dose H*(10) values evaluated at 55 cm from the Am–Be source, derived by using Maxed and Gravel codes; the ISO Am–Be spectrum was introduced as a default spectrum. The same quantities derived from the Bayesian analysis are also reported. Neutron field parameters
Maxed code
Gravel code
Bayesian analysis
ftot (cm 2) H*(10) (mSv h 1)
64.7 90.5
63.4 87.4
61.2 73.0 89.2 75.0
Hn ð10Þ ¼
Z
hðEÞfE ðEÞdE
X
hf l l l
ð4Þ
where h(E) is the appropriate fluence to dose conversion coefficient for neutrons at energy E and hl the corresponding discretized coefficient. Dose conversion coefficients are extracted from ICRP74 [12]. The values of each integral quantity deduced from both calculations are reported in Table 2 These values are in fairly good agreement, with average values equal nearly to 64 n/cm2 and 89 mSv h-1 for the total neutron fluence ftot and the ambient equivalent dose H*(10) respectively and discrepancies lower than 3%. However, the uncertainties in these values cannot be unfortunately evaluated because the unfolding methods cannot provide a mathematically rigorous or realistic estimate of the uncertainties [13,14] on the deduced spectra and therefore on any quantity computed from these spectra. Indeed, one needs to consider not only the uncertainties of the spectrum that result from the propagation of the uncertainties in the measured values and response functions, but also the uncertainties due on one hand to the fact that the solution spectrum is not unique and on other hand to the sensitivity of the solution spectrum to the prior ‘‘default spectrum’’ included to initiate the deconvolution procedure. These difficulties in the uncertainties estimation can be overcome in the Bayesian approach described in the next section.
3. The Bayesian analysis 3.1. Description of the Bayesian approach The principle of the Bayesian approach [7] differs from that of the standard unfolding methods reported in the preceding section where it has been seen that the expected spectrum is extracted from the experimental data assuming some prior criteria. In the
G. Medkour Ishak-Boushaki, M. Allab / Nuclear Instruments and Methods in Physics Research A 689 (2012) 35–39
fl ¼ fl ðEl , ai Þ
ð5Þ
where El is the discretized energy and ai represent the physical parameters related to the neutron spectrum description which, a priori, take unknown values. The second step is the use of Bayes’ theorem and the marginalization equations to determine these ai parameters [13,14]. The last step consists in using the probability rules to derive a probability density for each quantity of interest and its related uncertainty, i.e the integral discretized quantities given in Eqs. (3)and (4), functions of the ai parameters, written as follows: X X ftot ðai Þ f ðE , a Þ and Hn ð10Þðai Þ h f ðE , a Þ ð6Þ l l l i l l l l i
3.2. Parametric representation of the spectrum The parameterization choice should be guided by both the physics of the experiment and the limitations of the spectrometer. In addition, the parameterization should be as simple as possible and detailed enough to account for the fine structure that can be resolved by the spectrometer. In the present study, the complexity of the investigated spectrum and the scarce available experimental data (only six measured sphere responses, n ¼6) make the choice of the model spectrum difficult. For these reasons, the Bayesian analysis was performed by testing different models spectra to reproduce the Am–Be spectrum. The investigated spectrum has been approached by a Gaussian or a sum of Gaussian shapes. The model giving the best chi-2 value with reliable uncertainties on the interesting integral parameters has to be retained. A model spectrum consisting of a combination of two Gaussian shapes with different magnitudes, widths and mean energies, considered as free parameters, has been first tested. Unfortunately, the obtained results were not reliable probably because the number of experimental data was not sufficient to determine the six free parameters. To reduce the number of these parameters, two approaches have been considered consisting in describing the model spectrum as a sum of two Gaussians characterized by either the same magnitude or the same width. Taking into account the assumed behavior of the investigated Am–Be neutron spectrum, the last approach has been excluded. The model spectrum retained for the Bayesian analysis is thus described by the sum of two Gaussians of same magnitude A and of different mean energies (E01 and E02) and widths (s1 and s2). In this model, the fluence can be written in a discretized form as ! ! ðE E Þ2 ðE E02 Þ2 fl ¼ fl ðEl , ai Þ ¼ A exp l 201 þ A exp l 2 ð7Þ 2s1 2s2
where ai A, E01, E02, s1, s2 represent the model parameters. The starting point of the analysis is to formulate the prior probability distribution P(ai9I) that the set of the model parameters takes the ai values given the information I that is available about the experiment before performing the measurements (in our case, I is the assumption that the spectrum has approximately a shape given by the sum of two Gaussians of same magnitude). To estimate the ambient dose and the total fluence values, with their associated uncertainties, one needs to calculate the posterior probability P(ai9Nk,I) that the parameters ai have particular values given the actually measured Nk data and the information I. The posterior probability P(ai9Nk,I) is derived by applying Bayes’ theorem. It is proportional to the prior probability P(ai9I) and to the likelihood function P(Nk9ai,I) which represents the probability that the data Nk would have been measured if the hypothesis on the ai parameters was true and given the information I [13]. It can be expressed as Pðai 9Nk ,IÞaPðai 9IÞPðNk 9ai ,IÞ
ð8Þ
In this study, the analysis was performed by assigning uniform prior probability densities (with lower and upper bounds) to the ai model parameters: P(ai9I) is supposed uniform with respect to A, E01, E02, ln s1 and ln s2, for it is well established [2] that in the case of a Gaussian of width s, the reference prior is uniform in the logarithm of s. ^ k ðN ^ k ¼ P Rkl f Þ) is a funcThe number of estimated counts N l l tion of the parameters ai through the parameterized and discretized spectrum fl as given in Eq. (2). Assuming that the measured counts Nk is large enough so that ^ k and variances2 , is a good the normal distribution, with mean N Nk approximation for the likelihood function P(Nk9ai,I) 2
" # 3 ^k 2 X Nk N 1 5 PðNk 9ai ,IÞ ¼ exp4 k 2 sNk
ð9Þ
where sNk is the Nk estimated uncertainty and once the prior probability and the likelihood function are determined, the posterior probability P(ai9Nk,I) is calculated by using Eq. (8). The posterior probability is given for the ambient dose value H*(10) by [13] PðHn ð10Þ9Nk ,IÞ ¼
Z
dðHn ð10Þ
X
h f ðE , ÞÞPð i 9Nk ,IÞd i l l l l i
a
a
a
ð10Þ
4
Fluence (arbitray unit)
Bayesian approach, the spectrum is first modeled in terms of a probabilistic law described by definite parameters. The prior probability first assigned to the model parameters is afterwards modified by the information provided by the measurements procedure. The final result in Bayesian analysis is not a particular point in parameter space but a probabilistic evaluation of the parameters that includes a full description of their uncertainties. The Bayesian approach is developed in the following to determine the ambient dose equivalent value H*(10) and the total fluence ftot integral quantities estimated in the foregoing section, given the set of Bonner spheres responses reported in Table 1. The first step in this approach assumes a probabilistic parameterized model for the discretized neutron spectrum (cf. Eq. (2)) in the form
37
ISO spectrum Bayes spectrum 3
2
1
0 0
2
4
6
8
10
12
Neutron energy (MeV) Fig. 2. Modeled neutron spectrum, described as a sum of two Gaussian shapes, derived from the Bayesian analysis compared to the typical standard spectrum (ISO spectrum).
38
G. Medkour Ishak-Boushaki, M. Allab / Nuclear Instruments and Methods in Physics Research A 689 (2012) 35–39
P(fluence) 0.2 0.15 0.1 0.05 0.0
P(H*(10)) 0.15 0.1 0.05 0.0 50.0
55.0 60.0 65.0 fluence (n/cm2)
70.0
80.0
90.0
100.0
H*(10)(microSv/h)
Fig. 3. Relative probabilities of the ambient dose values (a) and of the total fluence (b) derived from the Bayesian analysis.
and for the total fluence ftot by Z X dðtot l fl ðEl , ai ÞÞPðai 9N k ,IÞdai Pðftot 9N k ,IÞ ¼
ð11Þ
where the symbol d represents the delta function. 3.3. Results of the Bayesian analysis and discussion The Bayesian analysis leads to a probability distribution P(ai9Nk,I) for the ai parameters used to model the neutron Am–Be spectrum. A good estimate of each parameter (optimal value) is given by the point where its probability distribution peaks. An estimate of the spectrum is calculated from the Eq. (7) using the optimal parameters values. The assessment of the model parameters was performed with a program written using the software WinBUGS, version 1.4 [15]. WinBUGS is an interactive Windows version of the BUGS program for Bayesian analysis of complex statistical models using Markov chain Monte Carlo techniques, developed at the MRC Biostatistics Unit in Cambridge, UK. Fig. 2 gives a plot of the selected modeled Gaussian spectrum compared to the ISO Am–Be neutron spectrum. The Gaussian parameters have been derived by using the Bayesian analysis, which leads to a chi-squared value of about 0.85. It can be seen that the model reproduces approximately the shape of the Am–Be spectrum but not its fine details. A better description should probably be achieved by using a sum of three Gaussians or more if more experimental data have been available. Fig. 3 shows the posterior probabilities for the total fluence ftot and for the ambient dose equivalent H*(10) as estimated from Eqs. (10) and (11) respectively. The best estimate for these quantities, which remain the main results of this analysis, are given by the values where the probability distributions peak. Their uncertainties values are derived from the widths (FWHM) of the respective probability distributions. Table 2 summarizes the integral quantities and their associated uncertainties derived by applying the Bayesian approach compared to the values derived from the standard unfolding methods. It can be seen that a good agreement is achieved with discrepancies lower than 10% for both ftot and H*(10) values, as required in radiation protection applications. It is important to underline that the most important advantage of the Bayesian approach, compared to the classical unfolding methods, consists in deriving mathematically the rigorous uncertainties of the investigated physical integral quantities which cannot be estimated in the different deconvolution procedures.
4. Conclusion In neutron spectrometry using Bonner spheres spectrometer, the spectra derived by using standard unfolding methods and thereby some quantities computed from these spectra (e.g,
fluences and doses) are subject to uncertainties which should be accurately estimated. However, it is not simple to quantify these uncertainties due to the fact that on one hand there are many ways to get an approximate spectrum solution (no uniqueness of the solution) and on other hand the solution may depend on the guess spectrum introduced to initiate the unfolding procedure. In this paper, it is shown that the approach using a Bayesian analysis can overcome these difficulties and can provide rigorous uncertainties for integral quantities derived from the neutron spectrometric measurements even when scarce experimental data are available. As an example, we have considered the characterization of an Am–Be neutron source field described by only six BSS response measurements. In the applied standard unfolding methods, an ISO Am–Be spectrum (close to the investigated solution) was introduced to initiate the deconvolution procedure based on the maximum entropy and on the least square unfolding techniques, developed in Maxed and Gravel codes respectively. In the developed Bayesian approach, the model source spectrum chosen consists of the sum of two Gaussians of unknown characteristic parameters. The comparison of the fluence and the ambient dose equivalent values characterizing the neutron field, derived from the unfolded spectra and from the Bayesian approach, shows a good agreement, with discrepancies between integral quantity values deduced by both methods less than the 10% percentage required in the radiation protection field. It is worth noting that the remarks on the unfolding methods and the Bayesian approach, referring here to analysis of Bonner spheres spectrometer measurements, may be applied to other neutron spectrometric data.
References [1] M. Matzke, Radiation Protection Dosimetry 107 (2003) 155. [2] P.C. Gregory, Bayesian Logical Data Analysis for the Physical Sciences: A comparative Approach with Mathematica Support, Cambridge University Press, 2004. [3] H. Mazrou, T. Sidahmed, Z. Idiri, Z. Mokrani, S. Bedek, M. Allab, Radiation Measurements 43 (2007) 1095. [4] H. Mazrou, Z. Idiri, T. Sidahmed, M. Allab, Journal of Radioanalytical and Nuclear Chemistry 284 (2010) 253. [5] S.P. Tripathy, A.K. Bakshi, V. Sathian, S.M. Tripathi, H.R. Vega-Carrillo, M. Nandy, P.K. Sarkar, D.N. Sharma., Nuclear Instruments and Methods in Physics Research Section A 598 (2009) 556. [6] M. Matzke, Unfolding of Pulse Height Spectra: The Hepro Program System, Technical Report PTB-N-19, Physikalisch-Technische Bundesanstalt, 1994. [7] D.S. Sivia ‘Jr., Skilling. Data Analysis. A Bayesian Tutorial, second ed., Oxford University Press, Oxford, 2006. [8] M. Reginatto, P. Golhagen, S. Neumann, Nuclear Instruments and Methods in Physics Research Section A 476 (2002) 242. [9] F.W. Stallmann, Lsl-m2: A Computer Program for Least-Squares Logarithmic Adjustment of Neutron Spectra, Technical Report NUREG/CR-4349, ORNL/ TM-9933, Oak Ridge National Labortory, 1986. [10] ISO 8529-1, Norme Internationale, Rayonnements neutroniques de re´fe´rence, 2001.
G. Medkour Ishak-Boushaki, M. Allab / Nuclear Instruments and Methods in Physics Research A 689 (2012) 35–39
[11] ISO 8529-2, Norme Internationale, Rayonnements neutroniques de re´fe´rence, 2001. [12] ICRP 74, International Commission on Radiological Protection, Conversion coefficients for use in radiological protection against external radiation, publication 74, Pergamon Press, Oxford, 1997.
39
[13] M. Reginatto, A. Zimbal, Review of Scientific Instruments 79 (2008) 023505. [14] M. Reginatto, Radiation Protection Dosimetry 121 (2006) 64. [15] D. Spiegelhalter, A. Thomas, N. Best, D. Lunn, WinBUGS Version 1.4 User Manual, MRC Biostatistics Unit, 2003. /http://www.mrc_bsu.cam.ac.uk/bugsS.