- Email: [email protected]
An unfolding method leading to a positive solution only
Nuclear Instruments and Methods in Physics Research 228 (1984) 129-132 North-Holland, Amsterdam
AN UNFOLDING
METHOD
LEADING
129
TO A POSITIVE
SOLUTION
ONLY
Hiroshi SEKIMOTO
Research Laboratory for Nuclear Reactors, Tokyo Institute of Technology, O-okayama, Meguro-ku, Tokyo, Japan Received 2 December 1983 and in revised form 6 April 1984
A new unfolding method has been developed to minimize an objective function by adjusting the logarithm of the spectrum. This method gives a positive solution over the whole energy region. The solution oscillates much less than conventional solutions using the linear least-squares method, such as FERDOR, even without an oscillation damping term.
1. Introduction F o r n e u t r o n spectrum m e a s u r e m e n t s with liquid scintillators, c o m p u t e r codes such as F E R D O R [1], using the linear least-squares method, are usually employed to unfold a n energy spectrum from a pulse-height distribution. The most i m p o r t a n t p r o b l e m of this m e t h o d is that it allows a negative value as a solution, a l t h o u g h one k n o w s that the spectrum value should not be negative for any energy. This causes large oscillations of the spectrum. In the m e a s u r e m e n t of a spectrum that has a s h a r p d o m i n a n t peak, this peak causes violent oscillations, a n d the part of the spectrum whose value is small becomes unreliable. In the present paper, the logarithm of the spectrum, not the spectrum itself, is adjusted to minimize an objective function. The solution then is positive for any energy, and oscillates m u c h less than in the conventional code.
Pm= ~. am,x,.
W h e n the m e a s u r e d c o u n t of the n th channel is c M, a n d its error is s,, the following function m a y be most naturally chosen as the objective function to be minimized: N J = Y~. c M - c~ )2, (1) n=] where N is the total n u m b e r of channels, a n d c c is the guess for the count of the n th channel, which calculated as
L
(2)
/=1 where A,~, n = 1 . . . . . N, is the response function for the 0 1 6 8 - 9 0 0 2 / 8 4 / $ 0 3 . 0 0 © Elsevier Science Publishers B.V. ( N o r t h - H o l l a n d Physics Publishing Division)
(3)
l=l
The most i m p o r t a n t problem of this m e t h o d is to a d m i t any value over the interval ( - ~ , ~ ) for x t. T h e n x t usually takes negative values for several /. These negative values m a y cause other positive x t to become too large in order to minimize the objective function, a n d x / oscillates violently. Negative values are not theoretically permitted for x t. The violent oscillations cause real spectrum structure obscurities. In F E R D O R , an oscillation d a m p i n g term is added to J given in eq. (1), a n d the following function is employed for the objective function:
J F = n=l ~
2. Objective function
c~ = E A.,x,,
/ t h n e u t r o n energy. The total n u m b e r of response functions is L. The value x / m a y be considered as an integrated spectrum over the l t h energy interval. The final solution is usually o b t a i n e d b y b r o a d e n i n g x~ with a p r o p e r G a u s s i a n distribution, G,a, as L
~ (cy-cf)2+~'~(xlt2'l~l \q']
(4)
where
Cy +Sn
q, = rain,, A,,----~ '
(5)
where the c o n s t a n t ~" in eq. (4) is determined from the b a l a n c e of the first a n d the second term of this equation. The second term works as if the zero is used as a n assumed spectrum, a n d should not be used from a physical point of view [2]. In the present paper, a new m e t h o d is proposed to minimize the objective function defined by eq. (1) subject to the condition x t > 0. Previously the value x t itself is adjusted, but in the present m e t h o d the value Yt defined by x / = e -v~
(6)
130
H. Sekimoto / Unfolding method to obtain a positive solution
is adjusted. Though y I may take any value over ( - oo,oo), x t takes a value only over (0, oo). Since the normal equation for this case is nonlinear, an iterative procedure should be employed. The details of this procedure are presented in the next section.
Marquardt's method [3,4] is employed to minimize the objective function J by adjusting the vector y. The following iteration scheme is derived with this method: y(i+ 1) = y(i) .~_ A y t i ) ,
(7)
where superscript (i) means the ith iteration, and A y (/) is obtained from the following equation: ~tti)A y(i) = #(i).
(8)
The components of the matrix a (i) and the vector fit/) are given by eqs. (9) and (10) respectively. The superscript (i) is omitted in the following for simplification. 10c c acc = Oy k Oyl + ~'Skl' n = l Sn
a k l = ~"
E
±(
f=
~/min Jo N-L
'
M c c. - c ° ,
OYl '
(9)
(10)
where OYt = A ' t x t "
(11)
The variable 8kl in eq. (9) is the Kronecker 8 symbol. This method approaches the gradient method for X ---, oo, and the G a u s s - N e w t o n method for X ~ 0, The parameter X is adjusted for each iteration.
If eq. (2) is an exact model and the error s n is evaluated accurately, the following relation holds from the statistics: min J --- N - L.
(12)
However, it is usually impossible to evaluate s n accurately, and eq. (2) is not an exact model with the errors of An/ and the discrete spectrum approximation. Then rain J usually becomes much larger than N - L when s n is treated to be equal to the statistical error of c M. In the present paper, s, is given by s, =frn,
(13)
where the statistical error is chosen as r, in the present work, but a more proper value for s, should be used for r, if it can be obtained. The error renormalization factor
(15)
By using (cov y ) , the covariance matrices of x and p are given by (cov x ) , t = x , ( c o v y ) k t x / ,
(16)
(COV p ) = G ( c o v x ) G t,
(17)
where the superscript t denotes a transpose operator. The solution x, and therefore p, does not depend on the value of f. The values of a for the corrected f and a 0 for f = 1 satisfy the following relation ~t =f2ot0.
(18)
For the first stage in the computer program, y, x and p with their covariance matrices are calculated for f = 1. The solutions y, x and p for f = 1 are the final solutions for corrected f. Their covariances for corrected f are obtained by modifying the covariances for f = 1 with f given by eq. (14).
5. Demonstrations In this section the spectra unfolded with the U F O / L code, which was programmed using the present method, are compared with the ones with F E R D O R for several problems. The directly compared values are pm+ and p,,~ defined as
P+
=Pro -+ ore,
(19)
where om = ¢ ( c o v P)mer,"
4. Error analysis
(14)
where the subscript 0 means for f = 1. The covariance matrix of y is given by (COV y ) = 0t 1.
3. Solution
n = l $2
f is chosen to hold eq. (12) and is calculated from
(20)
The negative values given by F E R D O R are omitted in the following figures. The unfolded spectra, where the response function A,l. for 14.3 MeV neutron energy is chosen as the measured pulse-height distribution c ~ , are shown in fig. 1. In the U F O / L calculation, though the value of J0 does not converge, the iteration is terminated when J0 becomes less than N - L. If the iteration is continued further, the value of J0 becomes much less than N - L. The exact solution of this problem is such that x t. for the 14.3 MeV neutron energy only takes nonzero value and all the other x t take zero. In F E R D O R , the second term of eq. (4) causes xt. to be smaller, then the other x / takes nonzero value to balance the smaller xt.. In U F O / L , there is no term corresponding to the second term of eq. (4), then any x t other than x t. takes a very small positive value. The error given by U F O / L is much smaller than in F E R D O R .
131
H. Sekimoto / Unfolding method to obtain a positive solution ' ~'
'
'
I
'
'
'
'
I
'
'
'
'
1.0
UFO/L FERDOR
~"
k. .t.,
UFO/L
I ~
.0
<
~
o.1
....f
0
> @
tl It
\k,, ~ :
'X ~"J
/
/
/ 1
/
0
>
C
ob.
®
@
z
Z 0
5
10
Neutron
Energy
15
0
T h e following examples are for actually m e a s u r e d pulse-height distributions. The first is for the n e u t r o n s p r o d u c e d in D - T reactions, which m a y contain some of the b a c k g r o u n d components. The s p e c t r u m unfolded by U F O / L oscillates m u c h less than F E R D O R , especially n e a r the source peak (fig. 2). The error renormalization l
l
l
l
l
,
,
C
UFO/L
i
l I
;
0
,( ,I
11 , ] ,
, , [ 0 , , , I i , , , 5
Neutron
(MeV)
Fig. 1. Comparison of the spectra unfolded with U F O / L and FERDOR. The response function for 14.3 MeV neutron energy is chosen as the measured pulse-height distribution. The upper curves are p + and the lower curves are p - defined by eq. (19).
l
v
o
~fL,
0.01
t;;i
Ikl
•
0
A O
FERDOR
h..
j
F j
~
~_~
//~
Energy
,,
,M
15
10
(MeV)
Fig. 3. Comparison of the spectra unfolded with U F O / L and FERDOR. The pulse-height distribution was measured in a water assembly 3 mean-free-paths apart from the D - T neutron source. The upper curves are p+ and the lower curves are p defined by eq. (19). factor f for U F O / L is 6.3. This m e a n s that the experim e n t contains a considerable a m o u n t of errors in addition to the statistical error. The last example is the m e a s u r e m e n t in a water assembly 3 mean-free-paths apart from the D - T neutron source. The unfolded spectra are shown in fig. 3. U F O / L depresses the oscillation near the source peak in this case also. The peaks near 7 MeV are enhanced. They should be examined in detail to see whether they are real structures. The error renormalization factor is 6.7, which is as large as in the previous case.
1.0 L=
6. D i s c u s s i o n a n d c o n c l u s i o n s
0
> o
0.1
10
C
o t9
z
,,'Lt ,i i
5
Neutron
i
10
Energy
i
I
15
I
(MeV)
Fig. 2. Comparison of the spectra unfolded with U F O / L and FERDOR from the measured pulse-height distribution of the D - T reaction neutrons. The upper curves arep ÷ and the lower curves are p - defined by eq. (19).
A n unfolding program to present a positive-valued spectrum over the whole energy region was developed. This m e t h o d does not need the oscillation d a m p i n g term for the objective function. The spectra unfolded by this m e t h o d show far fewer oscillations than F E R D O R , especially near the sharp d o m i n a n t peak. The error renormalization factors for the second a n d third examples in section 5 are too large. If the covariances of the m e a s u r e d pulse-height distribution a n d the response functions can be o b t a i n e d accurately, better spectra m a y be calculated. In the present paper, the equations are derived u n d e r the condition that the covariance of c ~ has only diagonal elements. Similar equations can be derived for the general form of the covariance also [2].
132
H. Sekimoto / Unfolding method to obtain a positive solution
References [1] H. Kendrick and S.M. Sperling, GA-9882 (1970). [2] H. Sekimoto and N. Yamamuro, Nucl. Sci. Eng. 80 (1982) 101.
[3] P.R. Bevington, Data reduction and error analysis for the physical sciences (McGraw-Hill, New York, 1969) p. 235. [4] R. Fletcher, AERE-R-6799 (1971).
AN UNFOLDING
METHOD
LEADING
129
TO A POSITIVE
SOLUTION
ONLY
Hiroshi SEKIMOTO
Research Laboratory for Nuclear Reactors, Tokyo Institute of Technology, O-okayama, Meguro-ku, Tokyo, Japan Received 2 December 1983 and in revised form 6 April 1984
A new unfolding method has been developed to minimize an objective function by adjusting the logarithm of the spectrum. This method gives a positive solution over the whole energy region. The solution oscillates much less than conventional solutions using the linear least-squares method, such as FERDOR, even without an oscillation damping term.
1. Introduction F o r n e u t r o n spectrum m e a s u r e m e n t s with liquid scintillators, c o m p u t e r codes such as F E R D O R [1], using the linear least-squares method, are usually employed to unfold a n energy spectrum from a pulse-height distribution. The most i m p o r t a n t p r o b l e m of this m e t h o d is that it allows a negative value as a solution, a l t h o u g h one k n o w s that the spectrum value should not be negative for any energy. This causes large oscillations of the spectrum. In the m e a s u r e m e n t of a spectrum that has a s h a r p d o m i n a n t peak, this peak causes violent oscillations, a n d the part of the spectrum whose value is small becomes unreliable. In the present paper, the logarithm of the spectrum, not the spectrum itself, is adjusted to minimize an objective function. The solution then is positive for any energy, and oscillates m u c h less than in the conventional code.
Pm= ~. am,x,.
W h e n the m e a s u r e d c o u n t of the n th channel is c M, a n d its error is s,, the following function m a y be most naturally chosen as the objective function to be minimized: N J = Y~. c M - c~ )2, (1) n=] where N is the total n u m b e r of channels, a n d c c is the guess for the count of the n th channel, which calculated as
L
(2)
/=1 where A,~, n = 1 . . . . . N, is the response function for the 0 1 6 8 - 9 0 0 2 / 8 4 / $ 0 3 . 0 0 © Elsevier Science Publishers B.V. ( N o r t h - H o l l a n d Physics Publishing Division)
(3)
l=l
The most i m p o r t a n t problem of this m e t h o d is to a d m i t any value over the interval ( - ~ , ~ ) for x t. T h e n x t usually takes negative values for several /. These negative values m a y cause other positive x t to become too large in order to minimize the objective function, a n d x / oscillates violently. Negative values are not theoretically permitted for x t. The violent oscillations cause real spectrum structure obscurities. In F E R D O R , an oscillation d a m p i n g term is added to J given in eq. (1), a n d the following function is employed for the objective function:
J F = n=l ~
2. Objective function
c~ = E A.,x,,
/ t h n e u t r o n energy. The total n u m b e r of response functions is L. The value x / m a y be considered as an integrated spectrum over the l t h energy interval. The final solution is usually o b t a i n e d b y b r o a d e n i n g x~ with a p r o p e r G a u s s i a n distribution, G,a, as L
~ (cy-cf)2+~'~(xlt2'l~l \q']
(4)
where
Cy +Sn
q, = rain,, A,,----~ '
(5)
where the c o n s t a n t ~" in eq. (4) is determined from the b a l a n c e of the first a n d the second term of this equation. The second term works as if the zero is used as a n assumed spectrum, a n d should not be used from a physical point of view [2]. In the present paper, a new m e t h o d is proposed to minimize the objective function defined by eq. (1) subject to the condition x t > 0. Previously the value x t itself is adjusted, but in the present m e t h o d the value Yt defined by x / = e -v~
(6)
130
H. Sekimoto / Unfolding method to obtain a positive solution
is adjusted. Though y I may take any value over ( - oo,oo), x t takes a value only over (0, oo). Since the normal equation for this case is nonlinear, an iterative procedure should be employed. The details of this procedure are presented in the next section.
Marquardt's method [3,4] is employed to minimize the objective function J by adjusting the vector y. The following iteration scheme is derived with this method: y(i+ 1) = y(i) .~_ A y t i ) ,
(7)
where superscript (i) means the ith iteration, and A y (/) is obtained from the following equation: ~tti)A y(i) = #(i).
(8)
The components of the matrix a (i) and the vector fit/) are given by eqs. (9) and (10) respectively. The superscript (i) is omitted in the following for simplification. 10c c acc = Oy k Oyl + ~'Skl' n = l Sn
a k l = ~"
E
±(
f=
~/min Jo N-L
'
M c c. - c ° ,
OYl '
(9)
(10)
where OYt = A ' t x t "
(11)
The variable 8kl in eq. (9) is the Kronecker 8 symbol. This method approaches the gradient method for X ---, oo, and the G a u s s - N e w t o n method for X ~ 0, The parameter X is adjusted for each iteration.
If eq. (2) is an exact model and the error s n is evaluated accurately, the following relation holds from the statistics: min J --- N - L.
(12)
However, it is usually impossible to evaluate s n accurately, and eq. (2) is not an exact model with the errors of An/ and the discrete spectrum approximation. Then rain J usually becomes much larger than N - L when s n is treated to be equal to the statistical error of c M. In the present paper, s, is given by s, =frn,
(13)
where the statistical error is chosen as r, in the present work, but a more proper value for s, should be used for r, if it can be obtained. The error renormalization factor
(15)
By using (cov y ) , the covariance matrices of x and p are given by (cov x ) , t = x , ( c o v y ) k t x / ,
(16)
(COV p ) = G ( c o v x ) G t,
(17)
where the superscript t denotes a transpose operator. The solution x, and therefore p, does not depend on the value of f. The values of a for the corrected f and a 0 for f = 1 satisfy the following relation ~t =f2ot0.
(18)
For the first stage in the computer program, y, x and p with their covariance matrices are calculated for f = 1. The solutions y, x and p for f = 1 are the final solutions for corrected f. Their covariances for corrected f are obtained by modifying the covariances for f = 1 with f given by eq. (14).
5. Demonstrations In this section the spectra unfolded with the U F O / L code, which was programmed using the present method, are compared with the ones with F E R D O R for several problems. The directly compared values are pm+ and p,,~ defined as
P+
=Pro -+ ore,
(19)
where om = ¢ ( c o v P)mer,"
4. Error analysis
(14)
where the subscript 0 means for f = 1. The covariance matrix of y is given by (COV y ) = 0t 1.
3. Solution
n = l $2
f is chosen to hold eq. (12) and is calculated from
(20)
The negative values given by F E R D O R are omitted in the following figures. The unfolded spectra, where the response function A,l. for 14.3 MeV neutron energy is chosen as the measured pulse-height distribution c ~ , are shown in fig. 1. In the U F O / L calculation, though the value of J0 does not converge, the iteration is terminated when J0 becomes less than N - L. If the iteration is continued further, the value of J0 becomes much less than N - L. The exact solution of this problem is such that x t. for the 14.3 MeV neutron energy only takes nonzero value and all the other x t take zero. In F E R D O R , the second term of eq. (4) causes xt. to be smaller, then the other x / takes nonzero value to balance the smaller xt.. In U F O / L , there is no term corresponding to the second term of eq. (4), then any x t other than x t. takes a very small positive value. The error given by U F O / L is much smaller than in F E R D O R .
131
H. Sekimoto / Unfolding method to obtain a positive solution ' ~'
'
'
I
'
'
'
'
I
'
'
'
'
1.0
UFO/L FERDOR
~"
k. .t.,
UFO/L
I ~
.0
<
~
o.1
....f
0
> @
tl It
\k,, ~ :
'X ~"J
/
/
/ 1
/
0
>
C
ob.
®
@
z
Z 0
5
10
Neutron
Energy
15
0
T h e following examples are for actually m e a s u r e d pulse-height distributions. The first is for the n e u t r o n s p r o d u c e d in D - T reactions, which m a y contain some of the b a c k g r o u n d components. The s p e c t r u m unfolded by U F O / L oscillates m u c h less than F E R D O R , especially n e a r the source peak (fig. 2). The error renormalization l
l
l
l
l
,
,
C
UFO/L
i
l I
;
0
,( ,I
11 , ] ,
, , [ 0 , , , I i , , , 5
Neutron
(MeV)
Fig. 1. Comparison of the spectra unfolded with U F O / L and FERDOR. The response function for 14.3 MeV neutron energy is chosen as the measured pulse-height distribution. The upper curves are p + and the lower curves are p - defined by eq. (19).
l
v
o
~fL,
0.01
t;;i
Ikl
•
0
A O
FERDOR
h..
j
F j
~
~_~
//~
Energy
,,
,M
15
10
(MeV)
Fig. 3. Comparison of the spectra unfolded with U F O / L and FERDOR. The pulse-height distribution was measured in a water assembly 3 mean-free-paths apart from the D - T neutron source. The upper curves are p+ and the lower curves are p defined by eq. (19). factor f for U F O / L is 6.3. This m e a n s that the experim e n t contains a considerable a m o u n t of errors in addition to the statistical error. The last example is the m e a s u r e m e n t in a water assembly 3 mean-free-paths apart from the D - T neutron source. The unfolded spectra are shown in fig. 3. U F O / L depresses the oscillation near the source peak in this case also. The peaks near 7 MeV are enhanced. They should be examined in detail to see whether they are real structures. The error renormalization factor is 6.7, which is as large as in the previous case.
1.0 L=
6. D i s c u s s i o n a n d c o n c l u s i o n s
0
> o
0.1
10
C
o t9
z
,,'Lt ,i i
5
Neutron
i
10
Energy
i
I
15
I
(MeV)
Fig. 2. Comparison of the spectra unfolded with U F O / L and FERDOR from the measured pulse-height distribution of the D - T reaction neutrons. The upper curves arep ÷ and the lower curves are p - defined by eq. (19).
A n unfolding program to present a positive-valued spectrum over the whole energy region was developed. This m e t h o d does not need the oscillation d a m p i n g term for the objective function. The spectra unfolded by this m e t h o d show far fewer oscillations than F E R D O R , especially near the sharp d o m i n a n t peak. The error renormalization factors for the second a n d third examples in section 5 are too large. If the covariances of the m e a s u r e d pulse-height distribution a n d the response functions can be o b t a i n e d accurately, better spectra m a y be calculated. In the present paper, the equations are derived u n d e r the condition that the covariance of c ~ has only diagonal elements. Similar equations can be derived for the general form of the covariance also [2].
132
H. Sekimoto / Unfolding method to obtain a positive solution
References [1] H. Kendrick and S.M. Sperling, GA-9882 (1970). [2] H. Sekimoto and N. Yamamuro, Nucl. Sci. Eng. 80 (1982) 101.
[3] P.R. Bevington, Data reduction and error analysis for the physical sciences (McGraw-Hill, New York, 1969) p. 235. [4] R. Fletcher, AERE-R-6799 (1971).