Author’s Accepted Manuscript Evaluation of a neutron spectrum from Bonner spheres measurements using a Bayesian parameter estimation combined with the traditional unfolding methods H. Mazrou, F. Bezoubiri www.elsevier.com/locate/radphyschem
PII: DOI: Reference:
S0969-806X(17)30929-5 https://doi.org/10.1016/j.radphyschem.2018.02.014 RPC7773
To appear in: Radiation Physics and Chemistry Received date: 24 August 2017 Revised date: 14 February 2018 Accepted date: 18 February 2018 Cite this article as: H. Mazrou and F. Bezoubiri, Evaluation of a neutron spectrum from Bonner spheres measurements using a Bayesian parameter estimation combined with the traditional unfolding methods, Radiation Physics and Chemistry, https://doi.org/10.1016/j.radphyschem.2018.02.014 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Title page Names of the authors: Dr. Hakim MAZROU and Mr. Fethi BEZOUBIRI Affiliation(s) and address (es) of the author(s): Division de Physique Radiologique – Centre de Recherche Nucléaire d’Alger (CRNA), 02 Boulevard Frantz – Fanon, B.P. 399, Alger-RP, 16000, Algeria. E-mail address of the corresponding author:
[email protected] (Dr. H. Mazrou)
Evaluation of a neutron spectrum from Bonner spheres measurements using a Bayesian parameter estimation combined with the traditional unfolding methods H. MAZROU1 and F. BEZOUBIRI Division de Physique Radiologique – Centre de Recherche Nucléaire d’Alger (CRNA), 02 Boulevard Frantz – Fanon, B.P. 399, Alger-RP 16000, Algérie.
ABSTRACT In this work, a new program developed under MATLAB environment and supported by the Bayesian software WinBUGS has been combined to the traditional unfolding codes namely MAXED and GRAVEL, to evaluate a neutron spectrum from the Bonner spheres measured counts obtained around a shielded
241
AmBe based-neutron irradiator located at a Secondary Standards Dosimetry Laboratory
(SSDL) at CRNA. In the first step, the results obtained by the standalone Bayesian program, using a parametric neutron spectrum model based on a linear superposition of three components namely: a thermal-Maxwellian distribution, an epithermal (1/E behavior) and a kind of a Watt fission and Evaporation models to represent the fast component, were compared to those issued from MAXED and GRAVEL assuming a Monte Carlo default spectrum. Through the selection of new upper limits for some free parameters, taking into account the physical characteristics of the irradiation source, of both considered models, 1
Corresponding author :
[email protected]
1
good agreement was obtained for investigated integral quantities i.e. fluence rate and ambient dose equivalent rate compared to MAXED and GRAVEL results. The difference was generally below 4% for investigated parameters suggesting, thereby, the reliability of the proposed models. In the second step, the Bayesian results obtained from the previous calculations were used, as initial guess spectra, for the traditional unfolding codes, MAXED and GRAVEL to derive the solution spectra. Here again the results were in very good agreement, confirming the stability of the Bayesian solution.
Keywords: Bonner Spheres; Bayesian Method; GRAVEL; MAXED; Neutron spectrum; WinBUGS.
1. INTRODUCTION The main physical characteristics of a neutron field, in a given workplace, can be revealed, among other things, by the neutron spectrometry. Among the existing tools applied in such places, the Bonner Spheres Spectrometer (BSS) [Bramblett et al., 1960] is the preferred one and widely used in specialized institutions for radiation protection purposes [Vega-Carillo et al., 2012; Amgarou et al., 2013; Barros et al., 2013; Barros et al., 2014; Mazrou et al., 2016]. The main advantage of this system is its ability to cover the full range of neutron energies encountered in many applications. It has isotropic response, good gamma-ray discrimination, and is relatively easy to operate during in-situ measurements. However, it has few drawbacks like time consuming, particularly when a passive thermal neutron detector is used, the weight, and the need to an appropriate unfolding procedure to reveal the detailed neutron spectrum of interest and its associated dosimetric quantities for radiation protection purposes [Alevra and Thomas, 2003; Vega-Carrillo et al., 2010; Benites-Rengifo et al., 2014]. At a present day, the Nuclear Research Centre of Algiers (CRNA), through its specialized laboratory, has adopted the BSS, as a useful tool, to characterize the neutron fields encountered in some studied workplaces for diverse purposes [Mazrou and Allab, 2012; Mazrou et al., 2016]. In addition, the related unfolding codes, in combination, with a good starting Default Spectrum (DS) were used, successfully, till now to derive the neutron spectrum from the BSS measured data. These unfolding codes are based on MAXED and GRAVEL [Reginatto and Goldhagen, 1999; Reginatto et al., 2004]. However, the major drawback usually associated with those traditional unfolding codes in deriving the solution spectrum from the BSS count rates, is that the quality of the provided
2
final solution spectrum might be affected by the Default Spectrum (DS) required to initialize the computational iterative process. To correct for this shortcoming, one need to resort to alternative and/or complementary methods. Among them, novel methods based on Artificial Intelligence: Artificial Neural Networks and Genetic algorithms have been widely investigated these last years with some degree of success [Vega-Carrillo et al., 2006; Ortiz-Rodriguez et al., 2014; Shahabinejad et al., 2016; Alvar et al., 2017]. Additionally, those based on the parameterization of the spectrum as a combination of continuous functions to represent it ‘realistically’ in its different energy parts, have also demonstrated their usefulness and their proven success, in neutron spectrometry domain, according to the literature [Medkour Ishak-Boushaki and Allab, 2017; Barros et al., 2014; Amgarou et al., 2011; Reginatto, 2009; 2006]. These methods are the object of our interest in the present work. Consequently, an attempt has been made to develop, during this work, an alternative unfolding approach to the traditional ones, namely: the Bayesian method, with the purposes firstly to improve the neutron spectrum derived from the BSS count rates and ultimately to diversify our own unfolding capabilities. This method has the advantage to learn from the measured data and, particularly, when they are consistent enough. In addition, the use of a default spectrum, in this case, is not necessary to initialize the unfolding iterative process, and finally it provides proper uncertainties related to the investigated parameters. 2. BONNER SPHERE SPECTROMETER The measurements used in this study were obtained during a previous work Mazrou and Allab, 2012] where the BSS of the Nuclear Research Center of Algiers (CRNA) were used to characterize the neutron beam of the OB26 irradiation system installed in Secondary Standards Dosimetry Laboratory (SSDL) at a reference irradiation position selected at a source-detector distance (SDD) of 150 cm. The irradiator system, under consideration, is described in Figure 1 [STS, 1998]. It provides biological protection against the radiation produced by the
241
Am-Be neutron source. It
3
consists mainly of a cylindrical polyethylene (rPE=0.92g/cm ) container enveloped with thin layers of lead (ePb=10mm; rPb=11.3g/cm3), cadmium (eCd=1mm; rCd=8.6g/cm3) and stainless steel (eSS=5mm; rSS=7.9g/cm3), successively. The neutron source is positioned radially along the beam channel. The channel of 15 cm diameter is closed with a polyethylene stopper which 3
has to be removed when starting the irradiation process. An Amersham International PLC X14 capsule with its three pellets-sources [AEA, 2000] of 241Am-Be (a, n) of 185 GBq global activity (5Ci), is mounted within the irradiator. The source is doubly encapsulated in welded stainless steel of 1.2mm thickness. The outer dimensions of the capsule are 60mm long and 30mm in diameter. The radioactive source is placed in a vertically positioned channel and is lifted into irradiation position by means of motor drive. The source irradiation position is about 1.10 meter above the floor. Please insert Figure 1 here The spectrometric system in use is based on a central spherical 3He thermal neutron proportional counter with a nominal gas pressure of 200kPa. The used BSS consisted of five polyethylene spheres (3", 5", 8", 10" and 12", 1 inch = 2.54 cm) and a spherical 3He counter, type SP9, with and without cadmium. The mean polyethylene density was estimated to be (0.949±0.008) g cm–3, a value nearly equal (0.4%) to that (r=0.946 g cm–3) for which response matrix was provided in the literature [Alevra and Plostinaru, 2000]. So, since the difference, in density, is minimal, we do not need to apply a correction factor for the present response matrix. The count rates obtained for each sphere are shown on Figure 2. It is worth noting that the statistical uncertainty was maintained below 1% for all the spheres, except for the bare detector covered with cadmium where it was kept around 2%. Please insert Figure 2 here As it can be seen, the sphere diameters corresponding to 0 and 1 (inch) are labeled for the bare detector without and with cadmium (Cd-thickness, eCd=1mm), respectively. The results show that the maximum count rate is obtained for the sphere diameter of 8 inch, whereas, the lowest count rate is obtained for the irradiation configuration of bare detector covered with cadmium because of the low sensitivity of the 3He detector for epithermal and fast neutrons [Mazrou and Allab, 2012]. 3. UNFOLDING PROCEDURE 3.1 The unfolding problem: The deduction of the neutron spectrum from the Bonner sphere measurements counts is one of the typical problems encountered in multisphere spectrometry. This can be solved by applying unfolding procedures. Mathematical equations describing such systems are formally known as Fredholm integrals of first kind given by Equation (1) [Matzke, 2003]:
4
M i ± e i = ò Ri (E ) F(E ) dE
(i=1,….nd)
(1)
where Mi is the sphere i reading (counts) obtained mathematically by folding the response function Ri(E) of this sphere with the spectral fluence rate F(E ) both expressed as functions of the energy E, whereas nd is the number of used Bonner spheres and εi is the difference between the predicted and the measured fluence values for the sphere i. εi is related to the estimated standard uncertainty σi via the chi-square statistical parameter χ2 following the next expression given by Equation (2): c2 = å i
e i2 s i2
(2)
To obtain a numerical solution, the system of Equation (1) is usually written in terms of discrete system of equations described by Equation (3), where F E j (E j ) is the fluence rate in energy group j extending from energy interval Ej to Ej+1 (DEj) and Rij (E j ) represents the response function of the ith sphere to neutron of energy that corresponds to the jth energy bin: nE
M i ± e i = å Rij (E j ) F E j (E j ) DE j
(3)
j =1
Because the number of energy groups nE is usually greater, normally more than 50, than the number of Bonner spheres nd used, typically of the order 10, i.e. nE >> nd, Equation (3) has an infinite number of mathematical solutions and just few of them has physical meaning for multisphere spectrometry. In the literature, many unfolding codes, based on different mathematical concepts, were used to unfold the neutron spectrum of Equation (3). In this work, to solve this equation, we focused on three different approaches namely: the maximum entropy method programmed in the MAXED code [Reginatto et al., 2004]; the least squares method programmed in the GRAVEL code Reginatto et al., 2004], and the Bayesian method represented by the program developed under the MATLAB environment with the support of WinBUGS software [Lunn et al., 2000; Spiegelhater et al., 2003]. The strategy of work used herein will be detailed hereafter. 3.2 The computational methodology In the present work, the program developed with the support of the Bayesian software WinBUGS [Lunn et al., 2000; Spiegelhater et al., 2003] was used together with the traditional 5
unfolding codes namely MAXED and GRAVEL to extract the final neutron spectrum from the Bonner sphere measurements counts issued from the OB26 irradiation system. To reach this purpose, the Bayesian program was used here in two different ways: §
Firstly, the main features of the solution neutron spectrum were represented by a parametric model based on physical meaning. The results obtained, following an appropriate selection of the free model parameters using the Bayesian concept, were compared, thereafter, to those derived from MAXED and GRAVEL codes using a Monte Carlo (MC) default spectrum.
§
Secondly, a sensitivity study was carried out based, successively, on the selection of new limits for some free parameters in the fast neutron region, choosing another theoretical model (ex. evaporation) to represent the fast energy component for the assumed parametric model spectrum and ultimately, combining the use of the WinBUGS final results with both considered parametric models as guess spectra for MAXED and GRAVEL codes.
This methodology of work tends to demonstrate that with a judicious selection of a parameterized model together with its free parameters, describing well the main features of the desired neutron spectrum, consistent results can be obtained in both cases: by the developed program used either standalone or by using their output results as guess spectra for the tradional unfolding codes. Furthermore, the stability of the Bayesian solution is also verified in both cases. The detailed of the unfolding methods and assumptions assumed herein, during these calculations, are briefly described hereafter.
4. BAYESIAN PARAMETER ESTIMATION USING WinBUGS 4.1 Introduction to the Bayesian approach In this study, the unfolding procedure based on the Bayesian approach was applied using the software package WinBUGS (version 1.4. This later is statistical software for the Bayesian analysis using the Markov Chain Monte Carlo Methods, MCMC). It was developed through the BUGS project (Bayesian inference Using Gibbs Sampling), jointly with a team of UK researchers at the MRC Biostatistics Unit, Cambridge, and the Imperial College School of Medicine, London [Lunn et al., 2000; Lunn et al., 2009].
6
Practically, to make use of the Bayesian approach with WinBUGS, under a MATLAB environment, a new script has been developed, specially, for unfolding calculations in neutron spectrometry using the Bonner Spheres measurements. This program was based on MATBUGS [Murphy and Mahdaviani, 2010], an original MATLAB interface to WinBUGS. The Bayesian approach can be defined, in general, as a process of learning from measured data (or experience). It is based essentially on the theory of probabilities. Typically, in the Bayesian framework, one would look for the best predictive model representing a given physical system through the combination of information provided by both the probability statements made for the unknown model parameters and the knowledge of the measured data. The most frequently used quantity to estimate this model is known as the posterior density of probability. Thus, according to the Bayes theorem this quantity can be estimated by Equation (4) as follows [Sivia and Skilling, 2006]: ( ∣ ) ∝ ( ∣ ) ()
(4)
Where ∝ indicates the proportionality and ( ∣ ) is called the likelihood function (or
simply the likelihood), whether the () is known as the prior density and q represents the
unknown parameters of the model. The theorem therefore says that the posterior density P (θ | data), is proportional to the product of the prior density () and the likelihood. In the following, it will be given an insight on these quantities. 4.2 Components of the Bayesian analysis The key components of the Bayesian analysis using the new developed program, for neutron spectrometry, and supported by WinBUGS software are based on the following major steps: ˗
Introduce a parametric model of the neutron spectrum;
˗
Introduce prior distributions probabilities;
˗
Calculate the likelihood function;
˗
Use Bayes’ theorem to estimate the posterior probabilities.
The first two steps are performed using the new developed program and the last two ones related to the evaluation of the likelihood function and the posterior density of probability are carried out via the MCMC method using Gibbs sampling algorithm in WinBUGS software. The detail of each calculations step is presented hereafter. 4.2.1 The parametric model of the neutron spectrum
7
The starting point of this method is based on an introduction of a parametric model of the observable of interest based on its physical meaning. In this work, the neutron spectrum issued from the
241
Am-Be neutron irradiator system is considered as the observable of
interest. It is defined by Equation (5) as follows: Φ" = Φ (#, )
(5)
where θ is the unknown parameters vector. These unknown parameters are considered as random variables. Taking into account the physical characteristics of the desired neutron spectrum that have been previously revealed, by experiments and by MC calculations, containing nearly: 15% of thermal neutrons, 5% of epithermal neutrons and 80% of fast neutrons [Mazrou and Allab, 2012], and because there are no single theoretical models or functions capable, alone, of reproducing realistically, with reasonable uncertainties, a neutron spectrum encountered in common investigated workplaces where the energy range span over a few decades i.e. ranging from a few meV to a few tens of MeV. A model spectrum, expressed as the contribution of three components namely: thermal, epithermal and fast neutrons, has been proposed in the present study. Furthermore, it is worth to note that the same parametric model was used extensively and successfully, in neutron spectrometry, and seems to be adequate for the majority of the investigated workplaces [Tomas et al., 2004; Bedogni et al., 2007; Amgarou et al., 2009; Reginatto, 2009]. So, it can be written as Equation (6): Φ$ = Φ%" + Φ&" + Φ"
'
'
(6)
where Φ%" , Φ&" and Φ" are the thermal-Maxwellian distribution, epithermal (1/E behavior in
terms of neutron fluence) and a kind of fast-Watt fission distribution, respectively. Theirs mathematical formulations are defined by the following Equations (7-9) [Tomas et al., 2004; Bedogni et al., 2007]: "
Φ%" = % * . / -
Φ&"
01 2-
= & 81 − /
01. 1. :
# ≤ 0.1 /7 ; . # <>? . /
–1 @′
(7) 0.1 /7 ≤ # < 10 B/7
(8)
8
Φ" = ' # C / '
01 @
10 B/7 ≤ # ≤ 20 E/7
(9)
Where at, ae and af represent the magnitudes of the energy peaks located in the thermal, epithermal and fast regions, respectively. §
T0 : is the most probable energy of the Maxwellian thermal distribution (T0 = 0.025 eV) ;
§
Ed: indicates the lower energy of the epithermal region (Ed=0.0707 eV);
§
b and β': describe the slope of the ascending and descending part of the curve, respectively, and are adjusted parameters with reasonable limits (-0.5 < b < 0.5 and 0 < β' < 1),
§
α and β: describe the shape and the peak of the prompt fission spectrum, and are adjusted parameters with reasonable limits: 0<α<1 and 1<β<2.
In this work, the unknown parameters θ ≡ at, ae, af, b, β', α and β are subject to the Bayesian analysis. 4.2.2 The prior distributions probabilities As a first guess, an initial prior uniform density of probability (), with appropriate lower
and upper limits, was considered in this work for all of the free parameters [Reginatto, 2006; Amgarou et al., 2009]. 4.2.3 The likelihood function The likelihood ( ∣ ) is the information (or observation) that can be extracted from
FG (E FG = measurements. To estimate this function, one should first evaluate the count rates E ∑J IGJ KJ ) which can be obtained by multiplying the response matrix with the parameterized
and discretized fluence vector as expressed in Eq. (3).
In general, when the measured counting rates in the detectors (i.e. the Bonner spheres) are large enough it’s a good approximation to assume a normal distribution to the measured count
FG , σR ). E FG and SG represent the mean and the rates according to the formula EG ~ LMNOP (E
standard deviation of the measurements of the normal distribution, respectively. Therefore, under this assumption the likelihood can be expressed by Equation (10) as follows [Reginatto, 2006; 2010]: ?
( EG ∣ ) = /TU[− V ∑G(
FX V WX >W ) ] YX
(10)
9
Once the prior distribution and the likelihood function have been evaluated, it is straightforward to calculate the necessary posterior distribution using the Bayes’ theorem expressed in Eq.(4) above [Reginatto, 2010]. 4.2.4
The posterior probability
The estimation of posterior density of probability ( ∣ ), is performed using directly
the software WinBUGS. As stated previously, given the prior probability for all the unknown
parameters composing the parametric spectrum model and the likelihood of the measured data, one may derive the solution spectrum from the posterior density of probability obtained for the free parameters by choosing a best estimate quantity (i.e. the mean value) [Reginatto, 2006; 2009]. In this work, the neutron spectrum model is function of seven parameters. Now, once the fluence rate F E for each energy bin is determined the dosimetric quantity like the ambient dose equivalent rate, H * (10) , can be estimated given the fluence to ambient dose equivalent factors hj [ICRP-74, 1996] according to Equation (11): H * (10) = å h j F j
(11)
j
j is the energy group. 5. UNFOLDING USING THE TRADITIONNAL METHODS The unfolding process with the traditional codes used in this work, namely MAXED and GRAVEL are resumed hereafter: 5.1 Unfolding using MAXED Unfolding method in MAXED [Reginatto et al., 2004] is based on the following steps: ˗
DEF choose the best estimate of the spectrum F j to start the iterative calculations;
˗
from spectra approaching the solution, choose the F j , which is the closest as possible DEF to F j among all feasible solutions which maximize the entropy S defined by
Equation (12) [Reginatto and Goldhagen, 1999; Reginatto et al., 2004]: é æF S = - å êF j ln çç j DEF Fj j ë è
ù ö ÷÷ + F DEF -Fjú j ø û
(12)
10
where F j is the determined fluence rate in group energy j and F DEF is its corresponding j value in the discretized default spectrum. Running MAXED code always requires creating, beforehand, input files including information on: ˗
the response functions of the spectrometric system to be used;
˗
the collected spectrometric measurements;
˗
(prior information) to be used to initiate the iterative the default spectrum F DEF j calculation.
For more details the reader is referred to relevant references available in the literature [Matzke, 2003; Reginatto et al., 2004; Reginatto and Zimbal, 2008]. 5.2 Unfolding using GRAVEL At iteration k, the fit procedure is performed following the Equations (13-15) [Reginatto and Goldhagen, 1999; Reginatto et al., 2004] given hereafter: é nd k k ê å Wij ln M i / C i = F kj exp ê i =1 nd ê Wijk å ê i =1 ë
(
F kj +1
)ùú ú ú ú û
(13)
k +1 where F j is the value of the fluence rate in energy group j at the (k+1) th iteration;
M i is the count rate measured by the detector i; Cik is the count rate of the detector i calculated at iteration k; ng
(
C = å Rij exp ln F kj k i
)
(14)
j =1
Wijk is the relative contribution of energy group j in the detector response i during the run k:
Wijk = Rij F kj / Cik
(15)
The solution is obtained by minimizing the statistical parameter c2 versus (ln F i ) (i=1,…nd). To find the minimum of c2, the code uses a nonlinear least square method [Reginatto and Goldhagen, 1999; Reginatto et al., 2004]. To initiate the calculations process of equations (13-
11
15), a default spectrum is required. This later was taken from MCNP5 calculations as described in the following sub-section. 5.3 Default spectrum for MAXED and GRAVEL The default spectrum used herein for MAXED and GRAVEL to initiate the computational process was obtained from the OB26 irradiation system at a reference position of 150 cm from the geometrical centre of the neutron source using Monte Carlo (MC) calculations performed with MCNP5 [Mazrou and Allab, 2012]. The MC calculations have been performed using an energy structure consisting of 5 log-equidistant bins per decade to cover the whole interval of energy of interest (10-3 eV to nearly 20 MeV). This later was similar to that adopted for the matrix response of the BSS [Alevra and Plostinaru, 2000]. The MC calculations were performed assuming an ISO 8529/1 [ISO 8529-1, 2011] primary spectrum to represent the
241
Am–Be neutron source one. For more details on MC assumptions and
simulations, the reader can refer to the published papers in the literature [Mazrou et al., 2010a; 2010b]. Thus, the default neutron spectrum used herein is illustrated in Figure 3. Please insert Figure 3 here In this work, thermal neutrons are the particles below the ‘cadmium cut-off energy’, i.e. 0.5 eV, the epithermal ones were defined in the energy range between 0.5 eV to 10 keV and the fast neutrons were defined to be above 10 keV. Figure 3 reveals, clearly, the presence of a thermal part in the neutron spectrum due to the massive presence of hydrogenous material in the composition of the OB26 irradiator shielding and to the multiple neutron scattering occurring with concrete walls and the floor of the SSDL irradiation room. This is known as a room-return effect [Vega-Carrillo et al., 2008]. It is important to recall that our measurements were carried out very close to the floor. For more details, one can see the following references [Mazrou et al., 2010a; 2010b]. 6. ANALYSIS OF RESULTS From the Bayesian analysis performed with WinBUGS, the marginal distributions have been obtained for all (seven) free parameters considered in the parametric model spectrum. These densities of probability are illustrated in Figure 4. Three of them (at, ae and af) describe the magnitudes of the energy peaks located in the thermal, epithermal and fast regions, respectively. The parameters b and β' of the epithermal region describe the slope of the
12
ascending and descending part of the curve, respectively. While α and β describe the shape and the peak of the fission spectrum, with 0<α<1 and 1<β<2. Please insert Figure 4 here One can note that the posterior probability of b' parameter describing the epithermal region is, almost, identical to the uniform prior probability assigned earlier to this parameter. This is due, most probably, to the inconsistency of the measurements in this region. Besides, once the densities of probability for the free parameters were determined, the next step will be the evaluation of the integral quantities (i.e. the total fluence rate and the ambient dose equivalent). Figure 5 shows the posterior probabilities of these quantities. Please insert Figure 5 here ̇ = 55.4 ± The best estimate values, obtained from the median, are for the total fluence rate; ̇ ∗ (10) = 60.0 ± 0.3 (mSv/h). The 0.2 (n/cm2 s) and for the ambient dose equivalent rate; ! results are provided within 2s of uncertainty. The Bayesian results were, thereafter, compared to those derived from MAXED and GRAVEL using a Monte Carlo (MC) default spectrum as stated above. The results are shown on Figure 6. Please insert Figure 6 here From Figure 6, it can be seen that the different spectra are almost identical in shape. However, there is a slight shift appearing in the energy peak of the fast neutrons for the Bayesian spectrum (WinBUGS) comparatively to those obtained by MAXED and GRAVEL. This could be attributed to the guess spectrum and to the different unfolding methods used. The same observation was highlighted in a previous work where different unfolding codes were used [Barros et al., 2014]. In addition, one should recall that the use of Watt model with the free parameter (b) bounded between 1 and 2 MeV, is more appropriate to describe an
235
U based reactor prompt fission
neutron, since it provides the most probable fission energy peak around 2 MeV rather than representing an
241
Am-Be neutron spectrum with its mean energy around 4 MeV [ISO-8529,
2001]. Thus, to adapt the model to the present situation one should suggest a fine-tuning consisting in extending the higher limit of the Beta free parameter to 4 or 5 MeV and/or suggest a different model to represent properly this particular region (eg. Evaporation, fission-
13
Maxwellian, Madlan-Nix or Gaussian models, etc.). This will be detailed later in the sensitivity analysis. In all cases, one should keep in mind that the suggested new model should be described by no more than two free parameters as it has been stated earlier. In fact, to find a possible solution to the present system we must keep the number of all free parameters lower or equal than the number of experimental data (i.e. 7). Overall, compared to MAXED and GRAVEL [mean values for both codes are: for the total ̇ ∗ (10) = ̇ = 60.6 ± 0.1 (n/cm2 s) and for the ambient dose equivalent rate; ! fluence rate 65.0 ± 0.6 mSv/h], the Bayesian results are obtained within a difference of about 9% for both investigated integral parameters, which remains acceptable. The consistency of the unfolding procedure using WinBUGS is performed comparing the measured count rates to the calculated ones. This later is obtained by folding the solution spectrum with the response function of the Bonner Sphere Spectrometer. Figure 7 shows the ratio of the calculated to the measured readings. Please insert Figure 7 here The results obtained show that the calculated readings obtained using the Watt model for the fast energy neutrons are in agreement with most of the measurements except for the sphere diameter of 5" and 12" where the difference is around 10% and 20%, respectively. This is due probably to the inconsistency of data in these regions. The sphere diameter of one inch is for a bare detector shielded with cadmium. 7. SENSITIVITY ANALYSIS Since the fast neutrons component is the most dominate component of the neutron spectrum around the
241
Am-Be neutron source installed in the OB26 irradiator, we will try, as stated
previously, to improve the solution spectra by choosing appropriate model and related free parameters for the parametric model of the neutron spectrum proposed. Consequently, in this sensitivity analysis, a focus should be directed, successively, on the impact of the final results taking into account the following steps: ˗
Extend the beta upper limit value of the Watt fission model (Eq. 9) to 4 MeV;
˗
Choose an Evaporation model to represent the fast energy component for the parametric model spectrum [Bedogni et al., 2007];
14
˗
Combine use of the WinBUGS final results with both considered theoretical models as guess spectra for MAXED and GRAVEL codes.
7.1 Extension of the Beta free parameter upper limit of the Watt model As stated previously, the b free parameter of the Watt fission model was chosen between 1 and 4 MeV that’s to take into account the mean energy of the
241
Am-Be which is around 4
MeV [ISO-8529, 2001]. Similar methodology has been previously performed elsewhere [Bedogni, 2006]. Then, the posterior probability is newly computed, for this free parameter, with WinBUGS and illustrated in Figure 8 below. Please insert Figure 8 here Once the newly posterior probability is provided for the Beta free parameter, the following ̇ = best estimate values were obtained for the integral quantities i.e. for the total fluence rate ̇ ∗ (10) = 65.3±0.4 (mSv/h). 58.8 ± 0.3 (n/cm2 s) and for the ambient dose equivalent rate ! On can see that the final results were enhanced for both integral quantities compared to those found by MAXED and GRAVEL [mean values for both codes are: for the total fluence rate ̇ = 60.6 ± 0.1 (n/cm2 s) and for the ambient dose equivalent rate; ! ̇ ∗ (10) = 65.0 ± 0.6 mSv/h]. The difference in this case is less than 3%. 7.2 Evaporation model for the fast neutron component In this sub-section, the Watt fission model was replaced by and Evaporation one to represent the fast energy component of the parametric spectrum. The mathematical formula of this function is given by Equation (16) as follows: $
Φ# = %$
# &
+,
*'
(16)
T: is the evaporation temperature parameter expressed in MeV, 0.5 < T < 1.5 Again, similarly to the beta parameter of the Watt model and after sensitivity calculations, the free parameter T was newly bounded with a upper limit value of 4 MeV instead of the initial one set at 1.5 MeV, in order, to take into account the mean energy of the
241
Am-Be which is
around 4 MeV [ISO-8529, 2001]. The best estimate result for T was found around 3.2 MeV as shown in Figure 9. Please insert Figure 9 here
15
The influence of the variation of the evaporation temperature parameter T on the fast component of the spectrum is shown on Figure 10 hereafter. Please insert Figure 10 here As one can see, Figure 10 illustrates well the peak shift effect of neutron spectra obtained using an Evaporation model for the fast neutron component assuming an upper limit for the free parameter T of 1.5 MeV and of 4.0 MeV, respectively. As T increases, the peak tends to shift towards higher energy region. It is worthnoting that similar observations were highlighted previously in the litterature [Bedogni, 2006]. As expected, there is a much better agreement in the peak energy when considering the later best estimate evaporation temperature T (i.e. 3.2 MeV), compared to the peak energy provided by MC calculations for the SSDL's neutron irradiator. Thus, following this model, the neutron spectra using the three unfolding methods were obtained. The results are illustrated in Figure 11. Please insert Figure 11 here As it can be seen in Figure 11, the spectra obtained with different used methods are almost identical in shape, suggesting, to some extent, the reliability of the evaporation model used herein with the new upper limit imposed for the free parameter of interest. The consistency of the unfolding procedure using WinBUGS is performed comparing the measured count rates to the calculated ones. The results are illustrated in Figure 12. Please insert Figure 12 here The results obtained show that the calculated readings obtained using the Evaporation model for the fast energy neutrons are in agreement with most of the measurements except for the sphere diameter of 1" representing a bare detector shielded with cadmium. The measured data are probably not sufficient with only this irradiation configuration. Overall, the detailed results for both assumed models in the parametric neutron spectrum are resumed in Table 1. Please insert Table 1 here.
2Φ (MeV), reported in Table 1 is calculated using Equation The mean weight fluence energy / (17) as follows:
16
/2Φ =
∑5 478 #4 Φ4 ∑5 478 Φ4
(17)
Where N is the number of total energy bins (here, 53 energy bin structure i.e. 5 bins/decade has been adopted). While ℎ2 (pSv cm2) is the mean fluence to ambient dose equivalent conversion coefficient, and is determined according to Equation (18) as follows:
ℎ2 =
∑5 478 :4 Φ4 ∑5 478 Φ4
(18)
hi is the fluence to ambient dose equivalent conversion coefficient [ICRP-74, 1996] expressed in (pSv cm2) and reconstructed into 53 bins energy structure. As it can be seen, there is a good agreement for nearly all important investigated parameters for both models adopted in WinBUGS Bayesian software compared to MAXED and GRAVEL results (within 3-4% difference), except for the fluence rate at the epithermal region. This is due most probably to the inconsistency of measurements in this particular region. An additional irradiation configuration using a sphere diameter of 2" and/or 3" shielded with cadmium would be certainly an asset. Overall, these encouraging results suggest that the three unfolding methods converge on a single solution spectrum. 7.3 Bayesian Default Spectrum for MAXED and GRAVEL In this sub-section, the calculation of the integral quantities of interest has been done by the traditional codes MAXED and GRAVEL using, this time, WinBUGS results as default spectra with both considered models. The important parameters extracted from the obtained spectra are reported in Table 2. Please insert Table 2 here. The results obtained using MAXED and GRAVEL for investigated parameters are in agreement with each other for the considered cases. Therefore, the spectra unfolded with the Bayesian parameter estimation for both considered models could represent a good a priori data for the traditional unfolding codes. 8. CONCLUSION
17
During this work, an attempt has been made to develop successfully an alternative and a complementary unfolding method to the traditional ones, to evaluate a neutron spectrum from Bonner Spheres measurements at a given workplace. The suggested method is based on the Bayesian approach. To make use of this approach, a new program was developed under MATLAB environment and supported by WinBUGS software. To check the consistency of this method, we have chosen to evaluate the neutron spectrum issued from an OB26 irradiator at a source-detector-distance of 150 cm. Thus, a verification process was initiated comparing the results issued from WinBUGS program to those given by MAXED and GRAVEL codes. A parametric neutron spectrum was suggested, in the Bayesian approach, with two different models to represent the fast energy region: the Watt and Evaporation models. The thermal and epithermal parts were modeled by the Maxwellian and 1/E formulas, respectively. The results obtained, in the first step, from WinBUGS program with two models adopted herein and with a new upper limits obtained for their corresponding free parameters Beta and T taking into account the physical characteristics of the source (241AmBe with a mean energy around 4.0 MeV), were in very good agreement with those obtained by MAXED and GRAVEL. In the second step, the Bayesian results were used as guess spectra for the traditional unfolding codes, MAXED and GRAVEL. Here again the results were in good agreement confirming, thereby, the stability of the Bayesian solution. Overall, this study tends to demonstrate that with consistent measurements, the parametric model supported by the Bayesian method can constitutes an excellent alternative to the traditional unfolding codes and particularly, it can be used either in standalone or in combination with MAXED and GRAVEL, as an initial guess, to extract reasonably the appropriate neutron spectra in a given workplace. ACKNOWLEDGEMENTS This work is supported by the Nuclear Research Center of Algiers (CRNA). It is a continuation work of a previous CRNA Research Project aiming to characterize the neutron beam issued from an
241
AmBe-based source irradiator for radiation protection purposes. The
authors are grateful to the General Direction of the Nuclear Research Center of Algiers for their helpful administrative support.
18
Besides, the authors would like to express their heartfelt thanks to Dr Reginatto M. and Dr Zimbal A. from PTB Laboratory (Germany) for their help and, in particular, for providing us with their private documents. In fine, the authors are thankful to the anonymous reviewers for their very detailed expertise and for their useful comments and recommendations to improve the present manuscript. REFERENCES AEA QSA-TECHNOLOGY/GmBH, “Sources Manual” (2000). Alevra A.V. and Plostinaru V.D. (2000), Characterisation of the IPNE Bonner Sphere Spectrometry by comparison with the PTB System. PTB-6.42-00-1. PTB report, Braunschweig. Alevra A.V. and Thomas D.J. (2003). Neutron spectrometry in mixed fields: multisphere spectrometers. Radiation Protection Dosimetry, Vol. 107 No 1-3, 37-72. Alvar A.A., Deevband M.R., Ashtiyani M. (2017). Neutron spectrum unfolding using radial basis function neural networks. Applied Radiation and Isotopes, Vol. 129, 35-41. Amgarou K., Lacoste V., Martin A., Asselineau B., and Donadille L. (2009). Neutron spectrometry with a passive Bonner Sphere System. IEEE Transactions on Nuclear Science, Vol. 56 N°5, 2885-2895. Amgarou K., Lacoste V. and Martin A. (2011). Experimental characterization of the neutron spectra generated by a high-energy clinical LINAC. Nuclear Instruments and Methods in Physics Research A 629, 329–336. Amgarou K., Trocné M., Garcia-Fusté M.J., Vanstalle M., Baussan E., Nourreddine A., Domingo C. (2013). Characterization of the Neutron field from the
241
Am-Be isotropic
source of the IPHC irradiator. Rad. Measurements, Vol. 50, 61-66. Barros S., Gallego E., Lorente A., Gonçalves I. F., Vas P., Vega-Carillo H. R. (2013). Dosimetric assessment and characterization of the neutron field around a Howitzer container using a Bonner Sphere Spectrometer, Monte Carlo simulations and the NSDANN and NSDUAZ unfolding codes. Radiat. Prot. Dosim. 154(3), 346–355. Barros S., Mares V., Bedogni R., Reginatto M., Esposito A., Gonçalves I.F., Vas P., Ruhm W. (2014). Comparison of unfolding codes for neutron spectrometry with Bonner Spheres. Radiat. Prot. Dosim. 161(1-4), 46–52.
19
Bedogni R. (2006). Neutron spectrometry and dosimetry for radiation protection around a High energy electron/positron collider. Doctoral Thesis presented at the Faculty of Sciences of the University Autonoma de Barcelona (Spain). Bedogni R., Domingo C., Esposito A., and Fernández F. (2007). FRUIT: an operational tool for multisphere neutron spectrometry in workplaces. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 580(3), 1301-1309. Benites-Rengifo J.L., Vega-Carrillo H.R.; Velazquez-Fernandez J. (2014). Photoneutron spectrum measured with a Bonner sphere spectrometer in planetary method mode. Applied Radiation and Isotopes 83, 256-259. Bramblett R.L., Ewing R.I., and Bonner T.W. (1960). A new type of neutron spectrometer. Nuclear Instruments and Methods, 9(1), 1-12. International Commission on Radiological Protection (1996). Conversion coefficients for the use in radiological protection against external radiation. ICRP Publication 74. International Organization for Standardization Reference neutron radiations, Part 1: characteristics and methods of production. ISO 8529-1. ISO (2001). Lunn D., Thomas A., Best N., Spiegelhalter D. (2000). WinBUGS – A Bayesian modelling framework: concepts, structure, and extensibility. Statistics and Computing 10, 325– 337. Lunn D., Spiegelhalter D., Thomas A., Best N. (2009). The BUGS project: Evolution, critique and future directions. Statistics in Medicine. 28 (25): 3049–3067. Doi:10.1002/sim.3680 Matzke M. (2003). Unfolding procedures. Radiation Protection Dosimetry, (107) 155-174. Mazrou H., SidAhmed T. and Allab M. (2010a). Neutron field characterization of the OB26 CRNA irradiator in view of its use for calibration purposes. Applied Radiation and Isotopes, 114, 63-70. Mazrou H., SidAhmed T. and Allab M. (2010b). Neutron field characterization of the OB26 CRNA irradiator in view of its use for calibration purposes. Radiat. Prot. Dosim. 141(2), 114–126.
20
Mazrou H., and Allab M. (2012). Neutron field measurements of the CRNA OB26 irradiator using a Bonner Sphere Spectrometer for radiation protection purposes. Radiation Protection Dosimetry, 151(2), 354-364. Mazrou H., Nedjar A., and Seguini T. (2016). Neutron spectrum measurements at a radial beam port of the NUR research reactor using a Bonner spheres spectrometer. Applied Radiation and Isotopes, 114, 63-70. Medkour Ishak-Boushaki G. and Allab M. (2017). Accurate characterization of weak neutron fields by using a Bayesian approach. Applied Radiation and Isotopes, 122 (2017) 84– 88. Murphy K., Mahdaviani M. (2010). http://www.mrcbsu.cam.ac.uk/software/bugs/callingwinbugs-1-4-from-other-programs/. Ortiz-Rodriguez J.M., Reyes Alfaro A., Reyes Haro A., Cervantes Viramontes J.M., VegaCarrillo H. R., (2014). A neutron spectrum unfolding computer code based on artificial neural networks. Radiation Physics and Chemistry 95, 428-431. Reginatto M., and Goldhagen P. (1999). MAXED, a computer code for maximum entropy deconvolution of multisphere neutron spectrometer data. Health physics, 77(5), 579583. Reginatto M., Wiegel B., Zimbal A. and Langner F (2004). Analysis of data measured with spectrometers using unfolding techniques. NEA-1665/03 UMG 3.3. Reginatto M. (2006). Bayesian approach for quantifying the uncertainty of neutron doses derived from spectrometric measurements. Radiation Protection Dosimetry, 121(1), 6469. Reginatto M. and Zimbal (2008). A Bayesian and maximum entropy methods for fusion diagnostic measurements with compact neutron spectrometers. Rev. Sci. Instrum. 79, 023505 (2008). Reginatto M. (2009). What can we learn about the spectrum of high-energy stray neutron fields from Bonner spheres measurements? Radiation Measurements, 44, 692-699. Reginatto M. (2010). Overview of spectral unfolding techniques and uncertainty estimation. Radiation Measurements, 45(10), 1323-1329.
21
Shahabinejad H., Hosseini S.A., Sohrabpour M., (2016). A new neutron energy spectrum unfolding code using a two steps genetic algorithm. Nuclear Instruments and Methods in Physics Research Section A. Vol. 811, 82-93. Sivia D.S., Skilling J., 2006. Data analysis – A Bayesian Tutorial, second ed. Oxford University Press, Oxford. Spiegelhalter D., Thomas A., Best N., and Lunn D. (2003) WinBUGS User Manual. Available: http://www.mrc-bsu.cam.ac.uk/software/bugs/ STS, STEUERUNGSTECHNIK & STRAHLENSCHUTZ GmBH, “Neutron Calibrator OB26/2: Technical Specifications” (1998). Tomas M., Fernandez F., Bakali M. and Muller H. (2004). MITOM: A new unfolding code based on a spectra model method applied to neutron spectrometry. Radiation Protection Dosimetry, 110(1-4), 545-548. Vega-Carrillo H. R., Hernández-Dávila V. M., Manzanares-Acuña E.,.Mercado Sánchez G. A., Iñiguez de la Torre M. P., Barquero R., Palacios F., Villafañe R. M., Arteaga T. A., Ortiz Rodriguez J. M., (2006). Neutron spectrometry using artificial neural networks. Radiation Measurements 41(4), 425–431. Vega-Carrillo H.R., Gallego E., Manzanares-Acuna E., Lorente A., Barquero R., MartinMartin A., and Gutierrez-Villanueva J.L., (2008). Neutron room return effect. Revista Mexicana de Fisica S 54 (1) 63–66. Vega-Carrillo H.R., Ortiz-Hernandez A., Hernandez-Davla V.M., Hernandez-Almaraz B., Rivera Montalvo T., (2010). H*(10) and neutron spectra around linacs. Journal of Radioanal. Nucl. Chem. 283: 537-540. Vega-Carrillo H. R., Ortiz-Rodriguez J.M., Martinez-Blanco M.R., (2012). NSDUAZ unfolding package for neutron spectrometry and Dosimetry with Bonner spheres. Appl. Radiat. Isot. 71, 87–91.
22
FIGURE CAPTIONS
Figure 1
Schematic (a) and general (b) view of 241Am-Be neutron calibrator OB26 installed at the Algerian Secondary Standard Dosimetry Laboratory (SSDL) calibration room. (Dimensions are in mm, in Fig. 1a) [STS, 1998]
Figure 2
BSS measured count rates obtained from the OB26 neutron calibrator at a reference irradiation position selected at a source-detector distance (SDD) of 150 cm (the statistical uncertainty was below 1 % for all the spheres, except for the bare detector covered with cadmium where it was kept around 2 %).
Figure 3
MCNP5 fine energy (05 bins/decade) spectrum for 241Am–Be based-neutron irradiator for SSDL configuration at the selected SDD of 150 cm [Mazrou and Allab, 2012].
Figure 4
The marginal distributions for all the free parameters considered in the parametric spectrum model.
Figure 5
Posterior probabilities of the total fluence rate and the ambient dose equivalent rate derived from CRNA BSS measurements using the Bayesian parameter estimation.
Figure 6
Neutron energy spectra of the CRNA-OB26 neutron calibrator obtained at an irradiation position (SDD=150 cm) using WinBUGS, with a Watt model for fast neutron component, MAXED and GRAVEL with a M.C. Default Spectrum (DS).
Figure 7
Ratios of the calculated readings WinBUGS (Watt model) to the measured readings (BSS), derived from neutron spectrum shown in Figure 6.
Figure 8
The posterior probability for the Beta free parameter considered in Watt fission model.
Figure 9
The posterior probability for the free parameter T considered in the Evaporation model.
Figure 10
The peak shift effect due to the variation of an upper limit of the evaporation temperature free parameter T (from 1.5 MeV to 4.0 MeV) on the fast component of the neutron spectrum.
Figure 11
Neutron energy spectra of the CRNA-OB26 neutron calibrator obtained at an irradiation position (SDD=150 cm) using WinBUGS, with an Evaporation model for fast neutron component, MAXED and GRAVEL with a MC Default Spectrum (DS).
Figure 12
Ratios of the calculated readings WinBUGS (Evap. model) to the measured readings (BSS), derived from neutron spectrum shown in Figure 11.
23
(a)
(b) Figure 1 :
24
Figure 2 :
(MAZROU and BEZOUBIRI)
25
Figure 3 :
(MAZROU and BEZOUBIRI)
26
Figure 4 : (MAZROU and BEZOUBIRI)
27
Figure 5 :
(MAZROU and BEZOUBIRI)
28
Figure 6 :
(MAZROU and BEZOUBIRI)
29
Figure 7 :
(MAZROU and BEZOUBIRI)
30
Figure 8 :
(MAZROU and BEZOUBIRI)
31
Figure 9 :
(MAZROU and BEZOUBIRI)
Figure 10 :
32
(MAZROU and BEZOUBIRI)
Figure 11 :
(MAZROU and BEZOUBIRI)
33
Figure 12 :
(MAZROU and BEZOUBIRI)
34
LIST OF TABLES
Table 1
Comparison of physical characteristics of neutron spectra issued from OB26 irradiator obtained with WinBUGS using Watt and Evaporation models for fast neutron component and, with MAXED and GRAVEL codes.
Table 2
Comparison between fluence rate, ambient dose equivalent rate, mean energy and spectrum averaged fluence-to-ambient dose equivalent conversion factor for OB26 irradiator obtained with MAXED and GRAVEL codes using different default spectra.
Table 1 ḟ (cm-2 s-1)
Codes
E*£5x10-7 5x10-7
̇ *(10)(mSv h-1)
E>10-2
Total
E£5x10-7
5x10-7
E>10-2
Ē h*(10) (MeV) (pSv cm2)
Total
WinBUGS 9.5 ± 0.3 (Watt Model)
1.6 ± 0.1
47.6 ± 0.3 58.8 ± 0.3 0.4 ± 0.0†
0.1 ± 0.0†
64.9 ± 0.4 65.4 ± 0.4
2.8 ± 0.1
308.7 ± 1.1
WinBUGS 9.8 ± 0.2 (Evaporation Model)
1.5 ± 0.1
46.7 ± 0.2 58.0 ± 0.3 0.4 ± 0.0†
0.1 ± 0.0†
64.2 ± 0.4 64.7 ± 0.4
2.6 ± 0.1
309.5 ± 1.0
BSS (MAXED)
9.2
3.0
48.4
60.6
0.4
0.1
64.9
65.4
3.1
300.2
BSS (GRAVEL)
9.1
3.3
48.3
60.7
0.4
0.1
64.1
64.6
2.8
295.7
-
-
-
-
-
-
4.2
391
241 Am-Be [ISO-8529, 2001] *E: Energy is given in MeV. † 0.0 refers to values << 0.001
Table 2 Codes
Default Spectrum E*£5x1 0-7
MAXED
ḟ (cm-2 s-1)
̇ *(10)(mSv h-1)
5x10E>10-2 Total E£5x10 7 -2 -7
Ē h*(10) 5x10E>10-2 Total (MeV (pSv 7
MCNP5
9.2
3.0
48.4
60.6
0.4
0.1
64.9
65.4
3.1
300.2
Watt Model
10.0
1.4
49.8
61.3
0.4
0.1
66.6
67.1
3.3
304.3
Evap. Model
10.1
1.3
49.7
61.1
0.4
0.1
66.2
66.7
3.2
302.8
9.7±0.6
1.9±1
Mean ± Std. Dev.
49.3±0. 61.0 ± 0.4±0.0 8 0.4
0.1±0.0
65.9±0. 66.4 ± 3.2 ± 302.4 9 0.9 0.1 ± 2.1
35
GRAVE L
MCNP5
9.1
3.3
48.3
60.7
0.4
0.1
64.1
64.6
2.8
295.7
Watt Model
10.0
1.5
49.5
61.0
0.4
0.1
66.5
67.0
3.3
304.9
Evap. Model
10.1
1.4
49.4
60.9
0.4
0.1
66.1
66.6
3.2
303.7
9.7±0.6
2.1±1.1
Mean ± Std. Dev.
49.1±0. 60.9 ± 0.4±0.0 7 0.2
0.1±0.0
65.6±1. 66.1 ± 3.1 ± 301.4 3 1.3 0.3 ± 5.0
*E: Energy is given in MeV.
HIGHLIGHTS ·
Neutron spectrum from
241
AmBe irradiator was estimated using a Bayesian
approach. ·
The solution spectrum was derived from BSS counts using WinBUGS Bayesian Software.
·
Parametric model with seven free parameters was assumed for the neutron spectrum.
·
The Bayesian spectrum was compared successfully with the GRAVEL and MAXED spectra.
36