- Email: [email protected]
Numerical integration of detector response functions via Monte Carlo simulations
Accepted Manuscript Numerical integration of detector response functions via Monte Carlo simulations K.J. Kelly, J.M. O’Donnell, J.A. Gomez, T.N. Taddeucci, M. Devlin, R.C. Haight, M.C. White, S.M. Mosby, D. Neudecker, M.Q. Buckner, C.Y. Wu, H.Y. Lee
PII: DOI: Reference:
S0168-9002(17)30622-8 http://dx.doi.org/10.1016/j.nima.2017.05.048 NIMA 59886
To appear in:
Nuclear Inst. and Methods in Physics Research, A
Received date : 28 March 2017 Revised date : 19 May 2017 Accepted date : 31 May 2017 Please cite this article as: K.J. Kelly, J.M. O’Donnell, J.A. Gomez, T.N. Taddeucci, M. Devlin, R.C. Haight, M.C. White, S.M. Mosby, D. Neudecker, M.Q. Buckner, C.Y. Wu, H.Y. Lee, Numerical integration of detector response functions via Monte Carlo simulations, Nuclear Inst. and Methods in Physics Research, A (2017), http://dx.doi.org/10.1016/j.nima.2017.05.048 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 proof before it is published in its final 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.
Color Compatible Version Click here to view linked References
Numerical Integration of Detector Response Functions via Monte Carlo Simulations
1
2
3 4
K.J. Kellya,1 , J.M. O’Donnella , J.A. Gomeza , T.N. Taddeuccia , M. Devlina , R.C. Haighta , M.C. Whitea , S.M. Mosbya , D. Neudeckera , M.Q. Bucknerb , C.Y. Wub , H.Y. Leea a
5
b
6
7
Los Alamos National Laboratory, Los Alamos, NM 87545, USA Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
Abstract Calculations of detector response functions are complicated because they include the intricacies of signal creation from the detector itself as well as a complex interplay between the detector, the particle-emitting target, and the entire experimental environment. As such, these functions are typically only accessible through time-consuming Monte Carlo simulations. Furthermore, the output of thousands of Monte Carlo simulations can be necessary in order to extract a physics result from a single experiment. Here we describe a method to obtain a full description of the detector response function using Monte Carlo simulations. We also show that a response function calculated in this way can be used to create Monte Carlo simulation output spectra a factor of ∼1000× faster than running a new Monte Carlo simulation. A detailed discussion of the proper treatment of uncertainties when using this and other similar methods is provided as well. This method is demonstrated and tested using simulated data from the Chi-Nu experiment, which measures prompt fission neutron spectra at the Los Alamos Neutron Science Center.
8
Keywords: Detector Response, mcnp, 6 Li-glass Detector, Neutron Detection
9
1. Introduction
10
The signal rate for a detector in an experiment is a function of the number of particles emitted from the
11
source and the response function of the detector. In a typical experiment, a particle (neutron, γ ray, ...) is
12
emitted from a source or target with an initial energy, E, and the experimental task is to detect that particle
13
and record some properties of it, such as energy, direction, etc. In some γ-ray experiments, the particle is
14
emitted from the target and potentially scatters along the way to the detector. If the γ ray reaches the
15
detector, it does so with an energy, E 0, that may or may not be equal to E. That γ ray then transfers an
16
energy, Et , into the detector which is recorded as a count in a bin of an output spectrum that corresponds
17
to Et . In the case of a γ-ray cascade, multiple γ rays may be emitted nearly simultaneously and coincidence
18
summing can be an issue. Neutron-detection experiments behave similarly with the difference that measured
19
energy, Et , of a neutron is commonly determined via a time-of-flight technique. In either type of experiment,
20
it is likely that Et 6= E or E 0. As such, determining the true distribution of initial energies, E, from a number
21
of particles measured to have a potentially different energy, Et , may be difficult and likely requires knowledge ∗
Corresponding Author Emailsubmitted address: to [email protected] (K.J. Kelly) Preprint Elsevier
May 19, 2017
22
of the detector response function. A detector response function for any kind of experiment must account for all of the above effects. In many experiments, including the experiment discussed in Section 3, experimental efficiencies of individual detectors can be low enough that the chance of multiple particles interacting in a time short enough as to create a coincidence-summing effect is negligible even in the presence of multiple particles being emitted from the target nearly simultaneously. Additionally, the particle multiplicity for the experiment discussed in Section 3 is typically below 3, which further reduces the chance for coincidence summing. Therefore, this potential aspect of detector response is not considered further in this work. As an initial attempt to describe the detector response, it is tempting to describe C(Y (E), Et ), the number of observed counts in a bin corresponding to the measured energy, Et , with an equation of the form C(Y (E), Et ) = ρ(Y (Et ), Et ) Y (Et ),
(1)
23
where Y is any yield relevant to an experiment and Y (E) and Y (Et ) are the value of Y at E and Et ,
24
respectively. This yield could be interpreted as the number of reactions per incident beam particle, the
25
number of particles emitted per unit time, or another pertinent definition. In the formalism of Eq. (1),
26
ρ(Y (Et ), Et ) is assumed to describe the response function, which includes experimental effects such as solid
27
angle. However, Eq. (1) represents a narrow view of detector response; while Eq. (1) may be able to describe
28
C(Y (E), Et ) with the proper ρ(Y (Et ), Et ), ρ depends on the detected energy, Et , as well as the initial energy
29
distribution. Thus, ρ(Y (Et ), Et ) is likely to be unique to each target or sample in an experimental setup, as
30
well as to a single target that emits particles with a different Y (E) under different experimental conditions,
31
and changes with the environment surrounding the detector as well.
32
It is more informative to write C(Y (E), Et ) in the form Z ∞ C(Y (E), Et ) = Y (E) Z0 ∞ × S(E, E 0, Et ) (E 0 ) dE 0dE,
(2)
0
33
34
35
where S is the scattering matrix (not to be confused with the nuclear physics S-matrix) describing the
probability for a particle emitted with an energy, E, to be scattered to an energy, E 0, and determined to have an energy, Et , and (E 0 ) is the detector efficiency corresponding to a particle of energy E 0. We can now define the response function R(E, Et ) as Z ∞ R(E, Et ) = S(E, E 0, Et ) (E 0 ) dE 0,
(3)
0
36
which yields C(Y (E), Et ) =
Z
0
∞
Y (E) R(E, Et ) dE.
(4)
Note that R(E, Et ) depends on E and Et but not on E 0 because the latter has been integrated out. Additionally, although the lower limit of the integral in Eq. (4) is set to 0, R(E, Et ) is essentially null below E ≈ Et . 2
It should also be noted that while Eq. (4) clearly resembles Eq. (1) of Ref. [1] and Eq. (1) of Ref. [2], it will be seen in the subsequent sections that the approach to this issue taken here is very different. Combining Eqs. (1) and (4) we see that ρ(Y (Et ), Et ) =
Z
1 Y (Et )
∞
Et
Y (E) R(E, Et )dE.
(5)
37
This shows that while ρ(Y (Et ), Et ) is related to R(E, Et ), they are fundamentally different quantities and
38
ρ(Y (Et ), Et ) is insufficient to describe the detector response. Even if a form of ρ(Y (Et ), Et ) that correctly
39
describes C(Y (E), Et ) for an energy distribution, Y (E), corresponding to a particular experimental environ-
40
ment is found, that form of ρ(Y (Et ), Et ) is no longer valid if the underlying Y (E) has changed, even if the
41
entire experimental environment remains the same.
42
2. Monte Carlo Simulations
43
In order to obtain a simulated spectrum of C(Et ) for a particular Y (E) it is common to simply run a
44
Monte Carlo simulation of the experimental environment with codes such as mcnp [3, 4] or geant4[5].
45
Running the simulation amounts to a numerical integration of Eq. (4) for the complicated response function
46
and can be time consuming for complex experiments, even with high-performance computers. The problem
47
gets much worse if a result is to be extracted from a data set by iteratively running simulations with various
48
Y (E) distributions until satisfactory agreement with the experimental data is obtained.
R
However, if R(E, Et ) is known, then the integral, i.e. Eq. (4), can be evaluated without running more simulations. Given a single simulation with high statistics obtained with E sampled from a specific initial yield distribution, Yo (E), it is possible to record the values of E and Et associated with each event in matrix form. As will be discussed in Section 3, the choice of what distribution should be used for Yo (E) is somewhat arbitrary and only the sampling rate of energies of interest needs to be considered. Whatever the chosen Yo (E) distribution, that choice is precisely defined and therefore carries no uncertainty in any subsequent calculations. Regardless of the choice of Yo (E), this matrix is exactly the product Yo (E) R(E, Et ), the integration of which over E yields the desired simulation spectrum. Furthermore, since the input Yo (E) distribution is known, it can be factored out of each simulated event, yielding the response function. It then follows that the simulation result for a different input distribution, Y 0 (E), can be obtained by simply integrating Eq. (4), instead of running a new simulation. This procedure yields the new simulation output Z ∞ C 0 (Y 0 (E), Et ) = Y 0 (E) R(E, Et ) dE. (6) Et
49
This integration can be done computationally in less than one second, as opposed to the hours or days that
50
may be required for a single Monte Carlo simulation. As is the case with any output spectrum from a
51
Monte Carlo simulation, the R(E, Et ) calculated in this way will have uncertainties related to the statistics
52
of the matrix, and also to uncertainties in the Monte Carlo simulation itself. For example, uncertainties in
53
the nuclear physics cross sections employed during the simulation and in the placement, size, and shape of 3
54
objects within the simulation can induce systematic uncertainties in R(E, Et ). While these systematic effects
55
are important for every experiment, they will not be discussed here. Even if systematic uncertainties are
56
ignored, the proper propagation of the statistical uncertainties on R(E, Et ) can be quite complicated. This
57
topic is discussed in Section 5.
58
The authors would also like to be clear in the fact that this method itself is not entirely new. Mea-
59
surements of neutron response to flight path environments are commonly carried out at facilities using an
60
incoming neutron beam [6–8], though these applications consider individual simulations when comparing to
61
measured results as opposed to creating a more widely applicable detector response matrix. Lajtai et al. [2]
62
have measured neutrons emitted from neutron-induced fission using a series of monoenergetic simulations
63
to describe their detector response. However, they never calculated a full response matrix. Instead, their
64
simulations were used to estimate their detector efficiency and to correct their data for small differences in
65
the measured time-of-flight neutron energy in an approximate way. The response matrix of fission reactor
66
cores have also been analyzed by running a series of hundreds of monoenergetic simulations to obtain a com-
67
parable description of the experimental environment [9]. Additionally, eigenmodes of reactor cores have been
68
decomposed in Ref. [10] in a similar manner. The present work is distinguished from these and other similar
69
applications primarily in that only a single, well-converged simulation was run to obtain the results shown
70
in subsequent sections, which is a fairly large simplification of the method. Finally, the work of Ref. [11]
71
and its application in Ref. [12] also represent a literature example where a single simulation is created and
72
manipulated to estimate physics parameters of an input distribution. However, the present work is unique
73
compared to those works as well in that a full description of the detector response function is provided here
74
along with crucial details discussed in Section 5 regarding the proper treatment of uncertainties when using
75
simulated spectra generated from a single R(E, Et ).
76
Although the above variations of the technique discussed in this work have been applied in certain areas
77
of physics, they are certainly not commonly employed and so the alternative derivation and description of the
78
method discussed in this work is warranted. Extensions of this method are likely also possible for analysis
79
of data from integral measurements. Furthermore, in Sections 3 and 4 we discuss the first application of
80
this particular method to a modern neutron time-of-flight experiment. Given that this method and similar
81
variants do not seem to be commonly employed in reactor or nuclear physics and are even more rarely, if ever,
82
employed in other areas of physics relying on Monte Carlo simulations, it is our hope that this publication
83
brings awareness to the time-saving potential of this technique.
84
3. Application to Simulated Data from the Chi-Nu Experiment
85
The goal of the Chi-Nu experiment at the Los Alamos Neutron Science Center (LANSCE) is to measure
86
the prompt fission neutron spectrum (pfns) for major actinides at varying incident-neutron energy ranges
87
[13]. To briefly describe the experimental setup, a pulsed, white neutron source is made incident on a parallel-
88
plate avalanche counter (ppac, see Ref. [14] for details) with target material deposited on an array of Ti foils.
4
6
Figure 1: A rendering of the Chi-Nu Li-glass detector array. Typically, the ppac target sits at the center of the array. Note R 6
that although only the Li-glass detector array is shown here, the mcnp simulation of the Chi-Nu experiment includes the floor, signal cables, and every other object included in the experimental environment.
89
Prompt fission neutrons are emitted from the target and detected with either a 6 Li-glass or a liquid scintillator
90
detector array [15]. 6 Li-glass and liquid scintillator arrays are intended to measure neutrons with E < 2 MeV
91
and E > 800 keV, respectively. Here we will focus only on the 6 Li-glass array, a rendering of which is shown
92
in Fig. 1, though treatment of the liquid array is analogous. In practice, data from 6 Li-glass detectors are
93
reliable up to ∼1 MeV because the 6 Li(n,α)t reaction cross section on which the 6 Li-glass detector efficiency
94
is based is only a standard up to 1 MeV [16]. However, the uncertainty on the 6 Li(n,α)t cross section reported
95
in ENDF/B-VII.1 [17] is below 2% up to approximately 2.5 MeV and so the simulated 6 Li-glass detector
96
spectra shown in this work are generally reported up to 2.5 MeV. Additionally, it is important to understand
97
that the pfns used in these calculations goes up much higher than this approximate upper energy limit
98
because neutrons with initial energies E > 2.5 MeV still contribute to the counts observed at Et < 2.5 MeV.
99
The experimental environment was modeled with mcnpx-PoliMi[3]. Neutron scattering, which results in
100
measured energies, Et , that are less than the true initial energies, E, adversely affects the Chi-Nu experiment
101
as well as other similar experiments. Scattering causes an excess in counts at low neutron energies and a
102
deficiency at high energies, which effectively distorts the measured detector efficiency and can distort the
103
extracted physics results as well. Significant efforts have been made to reduce this effect as much as possible
104
for Chi-Nu. During the course of characterizing the scattering effects of the Chi-Nu experiment, data from
5
105
older experiments were also discovered to be biased by uncorrected scattering as well [18]. Despite this
106
effort, data from the Chi-Nu experiment still contain a large amount of neutron scattering. This means the
107
response function must be accurately known to extract as much information out of the data as possible. The
108
application of the method described in the previous sections to data analysis at the Chi-Nu experiment has
109
significantly impacted how data are analyzed within our collaboration by reducing the amount of time it
110
takes to produce a single simulation spectrum to less than one second and reducing the weeks or months of
111
simulation CPU time required for a full data analysis to minutes.
112
In an experiment like Chi-Nu, neutrons are typically emitted from a target with a roughly Maxwellian
113
energy distribution. An example Maxwellian distribution with kT = 1.424 MeV is shown in Fig. 2a as the
114
solid line compared to literature evaluations of the neutron spectra from
115
as the dashed, dotted, and dash-dotted lines, respectively, and the ratios of these distributions to the kT =
116
1.424 MeV Maxwellian are shown in Fig. 2b to emphasize the differences between these spectra. In general,
117
the literature pfns distributions are within 5–10% of the Maxwellian spectrum. Furthermore, the initial
118
energy distribution observed during an experiment is not only a function of the target, but also of the energy
119
of the incoming neutrons that induce fission within the target. A pfns corresponding to Maxwellian with kT
120
= 1.3 MeV (as opposed to the 1.424 MeV Maxwellian shown in Fig. 2) was used as Yo (E) during integration
121
of Eq. (3) to obtain R(E, Et ) for the Chi-Nu experiment. The choice of this particular pfns was somewhat
122
arbitrary. It could have been chosen to be any distribution and the choice does not carry any additional
123
uncertainty with it since the chosen pfns is precisely defined. In principle, virtually the same R(E, Et ) can
124
be obtained regardless of the chosen Yo (E). The choice is motivated only by the desire to sample the initial
125
energies in accordance with how frequently they are encountered in an experiment, i.e. a pfns in the case
126
of Chi-Nu. A patch was made to mcnpx-PoliMi by M. Marcath [19] in order to obtain the initial neutron
127
energies for each simulated event so that Yo (E) could be factored out of Eq. (4). Following the steps described
128
in Section 2, the Maxwellian pfns used during the mcnp simulation was factored out of each event, yielding
129
R(E, Et ) for the 6 Li-glass array.
252
Cf,
235
U, and
239
Pu [17] shown
130
The resulting response function is shown in Fig. 3. A few comments are in order regarding the structure
131
of this response function. First, the dark band along E = Et corresponds to neutrons that did not scatter
132
along the path to the detector. This band is broadened according to the known detector resolution and light
133
curves, which allows for some counts to be accumulated with E < Et . The next most prominent feature
134
is the broad band of counts around E = 240 keV. These counts correspond to the 6 Li(n, α)t resonance at
135
approximately the same energy, with the lower Et counts in this band arising from neutrons emitted with the
136
resonance energy that then scattered one or more times along the path to the detector. The same statement
137
is true of all counts observed to be measured with energies Et < E, i.e. these are neutrons that were affected
138
by multiple scattering before detection. It is remarkable that this level of multiple scattering can still be seen
139
in data from the Chi-Nu experiment given that the experiment and the experimental environment were all
140
designed with an effort towards reducing the amount of multiple scattering in the data. Finally, it should be
6
0.3
PFNS (1/MeV)
PFNS Ratio to Maxwellian (kT = 1.424 MeV)
Maxwellian (kT = 1.424 MeV) 252 Cf (ENDF/B-VII.1) 235 U (ENDF/B-VII.1) 239 Pu (ENDF/B-VII.1)
0.35
0.25
0.2
0.15
0.1
0.05 10
−2
10
−1
Maxwellian (kT = 1.424 MeV) 252 Cf (ENDF/B-VII.1) 235 U (ENDF/B-VII.1) 239 Pu (ENDF/B-VII.1)
1.06 1.04 1.02 1 0.98 0.96 0.94 0.92 0.9 10
1
Initial Outgoing Neutron Energy (MeV)
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
(a) Examples of pfns distributions from commonly-studied actinides compared to a Maxwellian distribution with kT = 1.424 MeV.
(b) The same pfns distributions as shown in panel (a), but shown as a ratio to a Maxwellian distribution with kT = 1.424 MeV
Figure 2: Panels (a) showa a comparison of common pfns distributions and panel (b) shows the same distributions as a ratio to a Maxwellian of kT = 1.424 MeV. The solid, dashed, dotted, and dash-dotted lines are a Maxwellian with kT = 1.424 MeV, and 252 235 239 235 239 ENDF/B-VII.1 pfns distributions for Cf, U, and Pu, respectively. The U and Pu pfns distributions correspond to an incoming neutron energy of 1.0 MeV. These distributions are red, black, blue, and green online, respectively.
141
emphasized that R(E, Et ) depends on the entire environment surrounding the experimental area, in addition
142
to the detectors themselves and even to the target material under study. While it is likely true that, for
143
Chi-Nu, a change in target material from, for example,
144
the detector response matrix, the change is not zero. Therefore, even identical ppac’s with different target
145
materials require separately calculated response matrices. This environmental dependence can be seen in
146
Fig. 3. The thin bands seen at E ≈ 3.0, 2.5, 1.2, and 0.5 MeV and other energies correspond to resonances in
147
aluminum, oxygen, etc., as well as to objects in the experimental environment off which neutrons are likely to
148
scatter, thereby creating bands at specific times of flight [20]. Any change in the experimental environment
149
surrounding the Chi-Nu detector array will induce a change in R(E, Et ).
235
U to
239
Pu does not cause a significant change in
150
Given this R(E, Et ) obtained from a single, well-converged simulation, new Monte Carlo spectra were
151
created via Eq. (6) using various pfns distributions, each of which required time on the order of milliseconds
152
to create, and tests were carried out to confirm the accuracy of the new simulation output. The results of
153
these tests are discussed in the following section.
154
4. Results
155
Figure 4a shows some example Maxwellian pfns distributions at temperatures ranging from kT = 1.1-
156
1.4 MeV and Fig. 4b shows the corresponding spectra created using the R(E, Et ) shown in Eq. (6) and
157
Fig. 3. Note that the large peak at ∼240 keV present in each spectrum shown in Fig. 4 is caused by the
158
resonance in the 6 Li(n,α)t reaction at the same energy. These results demonstrate our ability to produce
159
simulation output for any desired input pfns. As a test of these and other new simulation spectra, it is
160
natural to simply compare simulation output created via the method described in this work to that of an
7
Initial Outgoing Neutron Energy, E (MeV)
10
103 1 102
10−1 10
10−2
10−1
1
10
1
Outgoing Neutron Energy from Time of Flight, E t (MeV)
6
Figure 3: The detector response function, R(E, Et ), for the Chi-Nu experimental setup using a Li-glass detector array. This 9 response function has not been normalized in any way and was calculated with a total of 1.98×10 neutrons emitted evenly over the ten targets within a single Chi-Nu ppac. The dashed line along E = Et shown to guide the eye. See the text for a discussion of the main features.
235
161
independently-calculated simulation. An independent simulation of a literature
U pfns [17] was carried
162
out for this comparison. The resulting simulation spectra from mcnp and from the convolution of the
163
pfns with R(E, Et ) are shown in Fig. 5a and the ratio of these spectra is shown in Fig. 5b. Both spectra
164
appear to be in excellent agreement. Although it may be noted that the statistical fluctuations become quite
165
large at lower time-of-flight energies where fewer counts are observed, the statistical uncertainty on these bins
166
is approximately 5%. This uncertainty aligns well with the variation about unity observed in Fig. 5b. It is
167
also worth noting that R(E, Et ) and the simulation shown in Fig. 5a were calculated using different random
168
number seeds in mcnpx-PoliMi.
235
U
169
A further test of the accuracy of the simulation output produced via Eq. (6) is to fit these output spectra
170
and use R(E, Et ) to extract the pfns used to create the calculated spectrum. Fits described here utilize
171
the TFractionFitter [21] class of root [22]. Detailed descriptions of this fitting routine have been provided
172
in previous works (see Refs. [23–26]), and so the description here will be kept brief. In short, a spectrum
173
is produced from R(E, Et ) for narrow ranges of initial neutron energy, E (i.e. an Et spectrum is produced
174
for each slice of the R(E, Et ) shown in Fig. 3 on a 10 bins per decade basis, in this case). Collectively,
175
these spectra span all E, are referred to as templates, and describe the portion of the detector response
176
corresponding to that particular E range. The goal is to determine the pfns used during creation of the
177
simulation output, which amounts to determining the relative contribution of each template spectrum to the
178
simulated data. The TFractionFitter class accomplishes this goal via a log-likelihood maximization routine
179
considering statistical variations in both the template spectra and in the data (real or simulated) being fit.
180
It is important to note that these fits do not assume anything about the shape of the pfns used to create the
181
simulation output. In this way these fits are a stringent test of the simulation output calculated via Eq. (6).
182
The resulting fit to simulation output created via a Maxwellian pfns of kT = 1.1 MeV is shown in Fig. 6a
183
with the input pfns as the solid line and the pfns from the fit as the diamonds. The same pfns and fit result
8
20000
Maxwellian kT = 1.1 MeV
18000
Maxwellian kT = 1.2 MeV
Maxwellian kT = 1.1 MeV Maxwellian kT = 1.2 MeV
0.4
Maxwellian kT = 1.3 MeV
16000
Maxwellian kT = 1.4 MeV
14000
Maxwellian kT = 1.3 MeV Maxwellian kT = 1.4 MeV
Counts
PFNS (1/MeV)
0.5
0.3
12000 10000 8000
0.2
6000 4000
0.1
2000
10
−1
10
1
−1
1
Outgoing Neutron Energy from Time of Flight (MeV)
Initial Outgoing Neutron Energy (MeV)
(b) A series of histograms created using the R(E, Et ) shown in Fig. 3 with the series of Maxwellian distributions on panel (a). The full circle, open diamond, closed diamond, and open circle data points correspond to temperatures of kT = 1.1, 1.2, 1.3, and 1.4 MeV (black, red, green, and blue online, respectively).
(a) A series of Maxwellian pfns distributions. The solid, dashed, dotted, and dash-dotted distributions correspond to temperatures of kT = 1.1, 1.2, 1.3, and 1.4 MeV (black, red, green, and blue, online, respectively).
Figure 4: The series of Maxwellian pfns distributions shown in panel (a) were used with the R(E, Et ) shown in Fig. 3 via Eq. (6) to create the histograms shown in panel (b). No normalization was applied to the simulated spectra because they were all created with the same R(E, Et ).
1.2
12000
Independent MCNP Simulation
1.15
Spectrum Calculated via Eq. (6) 6
Ratio of Simulation Spectra
14000
Li(n ,α )t Cross Section Trend
Counts
10000 8000 6000 4000 2000
1.1 1.05 1 0.95 0.9 0.85
10−1
0.8 1
Outgoing Neutron Energy from Time of Flight (MeV)
10−1
1
Outgoing Neutron Energy from Time of Flight (MeV)
(a) A comparison of an independent mcnp simulation to the simulation output calculated via the method described in this work.
(b) The ratio of spectra shown in panel (a). A dashed line is shown at unity for reference.
Figure 5: Panel (a) shows a comparison of an independent mcnp simulation spectrum obtained using an ENDF/B-VII.1 [17] 235 U pfns corresponding to an incoming neutron energy of 1.0 MeV shown as the solid line (black online, calculated with 8 9.9×10 neutrons emitted evenly over the ten targets within a single Chi-Nu ppac) to a spectrum calculated with Eq. (6) using 235 the same U input pfns convolved with R(E, Et ) shown as the open circle data points (red online). Excellent visual agreement 6 across the entire data spectrum. The Li(n,α)t cross section trend from ENDF/B-VII.1 [17] is shown as the dashed (blue online) 6 line as a reference for the structure seen at ∼240 keV in the Li-glass spectra. The ratio of the spectra shown in panel (a) is displayed in panel (b). See the text for a discussion of this panel.
9
Simulated Maxwellian (kT=1.10 MeV) PFNS from Fit
0.4
PFNS (1/MeV)
PFNS Relative To Maxwellian (kT=1.3 MeV)
0.45
0.35 0.3 0.25 0.2 0.15 0.1
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
1.6
Simulated Maxwellian (kT=1.10 MeV) PFNS from Fit
1.4
1.2
1
0.8
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
(a) Results of a fit to simulated data created via Eq. (6) using a Maxwellian of kT = 1.1 MeV.
(b) Results of a fit to simulated data created via Eq. (6) using a Maxwellian of kT = 1.1 MeV shown as a ratio to a Maxwellian of kT = 1.3 MeV.
Figure 6: A Maxwellian pfns with kT = 1.1 MeV is shown in panel (a) as the solid line (blue online). A simulation spectrum was created via Eq. (6) using this pfns and a fit was carried out to extract the pfns from the simulation data. The result of the fit is shown as the diamond data points (red online) in the same panel. The simulated Maxwellian pfns and fit result are shown in panel (b) as a ratio to a Maxwellian of kT = 1.3 MeV.
184
are shown in Fig. 6b as a ratio to a kT = 1.3 MeV Maxwellian to emphasize the details of this fit. The pfns
185
from the fit is nearly identical to the simulated Maxwellian distribution. An analogous test was also carried
186
out with a drastically different input pfns, a Maxwellian of kT = 1.65 MeV. Figure 7a shows the kT = 1.65
187
MeV Maxwellian distribution as the solid line with the extracted pfns from the fit as the diamonds and the
188
result is again shown as a ratio to a kT = 1.3 MeV Maxwellian in Fig. 7b. The kT = 1.1 MeV Maxwellian
189
from Fig. 6 is reproduced as the dashed line in Fig. 7 as a reference for the large change from one pfns to
190
the other. Again, a nearly identical reproduction of the simulated pfns was obtained from the fit.
191
Given that the fit results discussed in Figs. 6 and 7 as well as the simulation output data themselves
192
were obtained with the same R(E, Et ), it is not entirely surprising that a nearly exact reproduction of the
193
input pfns can be extracted. So, as a final test of this method, simulation data calculated directly with
194
mcnpx-PoliMi using a literature
195
and with different random seeds were fit as well. This input pfns is shown as the solid line and the fit pfns
196
is shown as the diamonds in Fig. 8a. The same results are shown in Fig. 8b relative to a Maxwellian pfns
197
of kT = 1.3 MeV to emphasize the details of Fig. 8a. The Maxwellian pfns distributions that were fit in
198
Figs. 6a and 7a are reproduced in Fig. 8 as the dashed and dash-dotted lines for reference. Clearly, there is
199
more variation in the resulting fit pfns of Fig. 8a than in Figs. 6 or 7.
200
235
U pfns [17] corresponding to an incoming neutron energy of 1.0 MeV
Although an increase in simulation time for both R(E, Et ) and the
235
U mcnp simulation would likely
201
improve the fit results shown in Figs. 8a and 8b because of the resulting smaller statistical uncertainties in
202
percent of the total counts in each bin, the variations in the extracted pfns are primarily caused by the
203
fitting routine itself. For the purposes of this work, the TFractionFitter class is an unfolding method [13]
204
in the sense that the input energy distribution, i.e. the pfns, is extracted without assuming any shape of
10
PFNS Relative To Maxwellian (kT=1.3 MeV)
Reference Maxwellian (kT=1.10 MeV)
0.45
Simulated Maxwellian (kT=1.65 MeV) PFNS from Fit
PFNS (1/MeV)
0.4 0.35 0.3 0.25 0.2 0.15 0.1
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
1.6
Reference Maxwellian (kT=1.10 MeV) Simulated Maxwellian (kT=1.65 MeV) PFNS from Fit
1.4
1.2
1
0.8
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
(a) Results of a fit to simulated data created via Eq. (6) using a Maxwellian of kT = 1.65 MeV. The kT = 1.1 MeV Maxwellian from Fig. 6a is shown for reference.
(b) Results of a fit to simulated data created via Eq. (6) using a Maxwellian of kT = 1.1 MeV shown as a ratio to a Maxwellian of kT = 1.3 MeV. The kT = 1.1 MeV Maxwellian from Fig. 6b is shown for reference.
Figure 7: A Maxwellian pfns with kT = 1.65 MeV is shown in panel (a) as the solid line (green online). A simulation spectrum was created via Eq. (6) using this pfns and a fit was carried out to extract the pfns from the simulation data. The result of the fit is shown as the diamond data points (red online) in the same panel. The simulated Maxwellian pfns and fit result are shown in panel (b) as a ratio to a Maxwellian of kT = 1.3 MeV. The kT = 1.1 MeV Maxwellian from Figs. 6a and 6b is reproduced here as the dashed line (blue online) for reference.
205
that distribution and using only knowledge of the detector response function. Unfolding methods have the
206
inherent property that statistical noise, however small it may appear in the measured or simulated data, is
207
amplified into an oscillatory variation in the resulting fit. This is most easily seen below 70 keV in Figs. 8a
208
and 8b. In the fitting performed to obtain the results shown in Figs. 8a and 8b, the small changes between
209
the 235 U simulated spectrum and R(E, Et ) caused by the different choice of random number seed are enough
210
to cause these strong variations in the fit result. The magnitude and character of these variations may be
211
different for other unfolding techniques, such as single value decomposition or Bayesian unfolding, but the
212
variations themselves will likely be present in any unfolding result. Furthermore, these variations can be
213
reduced by enforcing a smoothing condition on the unfolding, but this enforcement is somewhat artificial and
214
would likely make the results discussed in this work less clear. As such, no such smoothing was applied to
215
any fit described in this work. Despite these variations in the fits, it is still clear that the correct input pfns
216
has been extracted from this final test, thereby verifying that the correct mcnp spectrum was produced by
217
the method discussed in this work.
218
5. Notes on Uncertainties Associated with the Monte Carlo Simulation Output
219
When an independent simulation is run, the output spectrum can be assumed to have a Poisson distributed
220
statistical uncertainty, and so the statistical uncertainty on each bin can generally be taken to simply be the
221
square root of the number of counts. Given that the output of an independent simulation spectrum can be
222
nearly exactly reproduced using the R(E, Et ) shown in Fig. 3 and Eq. (6), it may seem natural to assume
223
that the statistical uncertainty on each bin of the spectrum calculated via Eq. (6) is also the square root of the 11
PFNS Relative To Maxwellian (kT=1.3 MeV)
Reference Maxwellian (kT=1.10 MeV)
0.45
Reference Maxwellian (kT=1.65 MeV) 235
0.4
U PFNS (ENDF/B-VII.1)
PFNS (1/MeV)
PFNS from Fit
0.35 0.3 0.25 0.2 0.15 0.1
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
Reference Maxwellian (kT=1.10 MeV)
1.6
Reference Maxwellian (kT=1.65 MeV) 235
U PFNS (ENDF/B-VII.1)
PFNS from Fit
1.4
1.2
1
0.8
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
(a) Results of a fit to simulated data created via Eq. (6) 235 using a U pfns from ENDF/B-VII.1 [17].
(b) Results of a fit to simulated data created via Eq. (6) 235 using a U pfns from ENDF/B-VII.1 [17]. shown as a ratio to a Maxwellian with kT = 1.3 MeV.
235
Figure 8: (a) A literature U pfns is shown as the solid line (black online) and the extracted pfns from the mcnp spectrum using different initial random number seeds than in Fig. 7a is shown as the diamonds (red online). The Maxwellian distributions used in Figs. 6 and 7 are shown as the dashed and dash-dotted lines (blue and green online) for reference, respectively. (b) A 235 reproduction of Fig. 8a shown relative to a Maxwellian of kT = 1.3 MeV in order to emphasize variations of the U pfns from the fit.
q q C(Y 0 (E), Et ). However, the simple C(Y 0 (E), Et ) statistical uncertainty
224
number of counts in that bin,
225
assumption does not apply to spectra calculated via Eq. (6) because all simulation output spectra calculated
226
using R(E, Et ) are correlated through the use of a single simulation. It is important to understand that only
229
the original, independently-calculated simulation spectrum used to obtain R(E, Et ) maintains the simple q C(Y 0 (E), Et ) uncertainty relationship.
230
Chi-Nu experiment. This pfns has no uncertainty associated with it because it is a somewhat arbitrary,
231
but precisely-defined choice for a probability density function from which initial energies are sampled and
232
that choice is motivated only by the requirement that initial energies are sampled in accordance with their
233
importance in the approximate expectation from experimental data, i.e. a pfns distribution. So, if we ignore
234
systematic uncertainties in the Monte Carlo calculation, the only source of uncertainty in spectra calculated
235
via Eq. (6) is the statistical uncertainty number of counts in each (E,Et ) bin of the original simulation
236
calculated with Yo (E), NE,Et . The associated covariance matrix for NE,Et is described by
227
228
As was mentioned earlier, an arbitrary Maxwellian pfns was used as Yo (E) to calculate R(E, Et ) for the
NEi ,Et
= σ(NE ,E )σ(N 0 0 Cov Ej ,Et )δEi Ej δEt Et , i t NEj ,Et0
(7)
237
where the δ’s in Eq. (7) are Kronecker deltas and σ 2 (NEi ,Et ) is the variance of NEi ,Et . Calculation of a
238
new simulation spectrum corresponding to any input pfns distribution, Yα (E), then amounts to scaling the
239
contents of each bin in R(E, Et ) by the ratio Yα (E)/Yo (E) and summing over all initial energies, E, that
12
240
contribute to each measured energy bin, Et , i.e. C(Yα (E), Et ) = =
Z
∞
0 n X i=1
Yα (E) R(E, Et ) dE
Yα (Ei ) N . Yo (Ei ) Ei ,Et
(8)
241
Using Eqs. (7) and (8), it can be shown that the covariance between C(Yα (E), Et ), C(Yβ (E), Et ), ...,
242
C(Yζ (E), Et ), i.e. the number of counts in the bin of a Monte Carlo simulation output spectrum corre-
243
sponding to Et calculated with Eq. (6) using various pfns distributions, Yα (E), Yβ (E), ..., Yζ (E), in place
244
of Y 0 (E), is given by the matrix σ2 C(Yα (E), Et ) α C(Yβ (E), Et ) σβα = Cov .. .. . . σζα
C(Yζ (E), Et )
with
σαβ =
n X Yα (Ei )Yβ (Ei ) i=1
Yo2 (Ei )
σαβ σβ2
··· ..
σζβ
.
···
σαζ
σβζ .. , . σζ2
(9)
σ 2 (NEi ,Et ).
(10)
245
The remaining elements of the above covariance matrix can be easily obtained by replacing the indices α and
246
β with the appropriate index in Eq. (10) keeping in mind that σαα = σα2 . Note that all pfns distributions
247
in the above equations are evaluated at the initial neutron energy, E, and not at Et . Also note that this
248
covariance matrix could be simplified by choosing a flat distribution instead of the chosen Yo (E). This may
249
be convenient for certain applications, but if the expected initial distribution is far from uniform, then this
250
simplification would be made at the cost of running the initial simulation for a much longer time in order to
251
obtain sufficient statistics at the important energies in R(E, Et ), as was the case with the Chi-Nu experiment
252
discussed in Sections 3 and 4.
253
The above covariance matrix assumes that each Et bin is independent of each other bin, Et0 . However, if
254
the desired new input distribution, Yµ (E), is defined by some function and the parameters describing that
255
function have uncertainties associated with them, then there may well be a non-zero covariance between each
256
value of Et . This is frequently the case for pfns distributions described by a functional form such as, for
257
example, a Maxwellian or a Watt distribution. In these distributions, changing a single parameter can change
258
the entire pfns across all initial energies and, in turn, across all measured energies, Et . For simplicity, we
259
will consider only two E bins, E1 and E2 , and two Et bins, Et and Et0 . If the covariance matrix for the
260
function describing the initial energy distribution Yµ (E) is given by
261
262
Yµ (E1 ) σµ2 (E1 ) σµ (E1 , E2 ) σ2 = = µ1 Cov Yµ (E2 ) σµ (E2 , E1 ) σµ2 (E2 ) σµ21
σµ12 2 σµ2
,
(11)
with all elements greater than zero, then the covariance between bins Et and Et0 of the simulation output 13
263
spectrum is given by C(Yµ (E), Et ) σ 2 (Et ) σ(Et , Et0 ) = . Cov C(Yµ (E), Et0 ) σ(Et0 , Et ) σ 2 (Et0 )
264
σ 2 (Et ) σ(Et , Et0 )
2 2 2 X X N N Y (E ) Ei ,Et Ej ,Et µ i σ 2 (NE ,E )+ = σ i t Y (E ) Y (Ei )Yo (Ej ) µij o i i=1 j=1 o
= σ(Et0 , Et ) =
2 X 2 N X Ei ,Et NEj ,Et0 i=1 j=1
Yo (Ei )Yo (Ej )
σµij
(12)
(13)
(14)
265
Note that σ 2 (Et0 ) can be obtained by replacing Et with Et0 in Eq. (13). The first term in the diagonal com-
266
ponents of this matrix corresponds to the uncertainty associated with NEi ,Et and NEi ,Et0 , but the remaining
267
diagonal terms and the off-diagonal terms arise from the covariance between the simulation output bins.
268
These additional components impact the assumed uncertainty in the resulting simulation output spectrum
269
and must be propagated through to the final covariance matrix.
270
It should be noted that although uncertainty propagation with simulation output spectra calculated via
271
Eq. (6) is significantly more complicated than a simple Gaussian or Poisson statistical approach, a much longer
272
simulation run time can be, but is not required to be, chosen for the R(E, Et ) calculation because there is
273
only one simulation to be run. This guarantees that the statistical uncertainty on the counts in each bin of
274
an output spectrum calculated via Eq. (6) is lower than would be obtained for each of potentially hundreds or
275
thousands of independent Monte Carlo simulations. Additionally, although the resulting covariance matrices
276
are not shown here, the uncertainties obtained with this method are even further reduced when considering the
277
multiplication or division of multiple simulation output spectra. In these cases, the covariance of the result
278
can be significantly lower than that obtained by multiplying or dividing simulation output spectra from
279
multiple independent calculations. Clearly, the covariances that must be considered in order to properly
280
handle the uncertainties on the simulation output produced via the method described in this work can be
281
quite complicated. Nonetheless, analytic expressions for each element of the necessary covariance matrices
282
can be derived, and thus the problem is tractable.
283
6. Conclusions and Future Work
284
The results shown in Section 4 demonstrate that the output of a Monte Carlo simulation for any input
285
distribution can be obtained by calculating the detector response function just once with a single simulation
286
and incorporating the desired input energy distribution in a subsequent integration of the detector response.
287
Furthermore, these output spectra can be calculated a factor of ∼1000× faster than running a new simulation.
288
The amount of time saved with this method will have a profound impact on any experiment that relies on
289
Monte Carlo simulations in which the input distribution is varied, but the experimental geometry remains
290
the same. While the uncertainties on the simulation output spectra created with this method are much
14
291
more complicated than a simple Gaussian or Poisson statistics uncertainty on each bin, the uncertainty on
292
each bin can actually be significantly smaller than would be obtained by running hundreds or thousands of
293
independent simulations. Despite the fact that variants of this method appear to be used primarily, though
294
not commonly, in reactor physics at the present time, its application can significantly impact research in
295
nuclear physics, nuclear astrophysics, and many other fields.
296
7. Acnknowledgments
297
298
The authors would like to thank M.E. Rising and F.B. Brown for fruitful discussions, M. Marcath for assistance with mcnp-related issues, and also P. Talou for helpful comments on this manuscript.
299
[1] T. N. Taddeucci, R. C. Haight, H. Y. Lee, D. Neudecker, J. M. O’Donnell, et al., Multiple-scattering
300
Corrections to Measurements of the Prompt Fission Neutron Spectrum, Nucl. Data Sheets 1232 (2015)
301
135.
302
303
304
305
306
307
308
309
310
311
[2] A. Lajtai, P. P. Dyachenko, V. N. Kononov, E. A. Seregina, Low-Energy Neutron Spectrometer and its Application for
252
Cf Neutron Spectrum Measurements, Nucl. Instrum. and Methods A 293 (1990) 555.
[3] S. A. Pozzi, E. Padovani, M. Marsequerra, A Monte-Carlo Code for Correlation Measurements, Nucl. Instrum. and Methods A 513 (2003) 550. [4] S. A. Pozzi, S. D. Clarke, W. J. Walsh, E. C. Miller, J. L. Dolan, et al., MCNPX-PoliMi for Nuclear Nonproliferation Applications, Nucl. Instrum. and Methods A 694 (2012) 119. [5] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo, et al., Geant4-A Simulation Toolkit, Nucl. Instrum. and Methods A 506 (2003) 250. [6] C. Coceva, R. Simonini, Calculation of the ORELA Neutron Moderator Spectrum and Resolution Function, Nucl. Instrum. and Methods A 211 (1983) 459.
312
[7] M. Mocko, G. Muhrer, T. Ino, M. Ooi, L. L. Daemen, et al., Monte Carlo study of the neutron time
313
emission spectra at the Manuel Lujan Jr. Neutron Scattering Center, Nucl. Instrum. and Methods A
314
594 (2008) 373.
315
316
317
318
319
320
[8] C. Guerrero, A. Tsinganis, E. Berthoumieux, M. Barbagallo, F. Belloni, et al., Performance of the Neutron Time-of-Flight Facility n TOF at CERN, Eur. Phys. J. A 49 (2013) 27. [9] R. G. Gamino, F. B. Brown, M. R. Mendelson, A Monte Carlo Green’s Function Method for ThreeDimensional Neutron Transport, Trans. Am. Nucl. Soc. 65 (1992) 237. [10] F. B. Brown, S. E. Carney, B. C. Kiedrowski, W. R. Martin, Fission Matrix Capability for MCNP Monte Carlo, Los Alamos National Laboratory Report LA-UR-13-26962.
15
321
[11] D. M. Schmidt, R. J. Morrison, M. S. Witherell, A General Method of Estimating Physics Parameters
322
from a Distribution with Acceptance and Smearing Effects, Nucl. Instrum. and Methods A 328 (1993)
323
547.
324
325
[12] G. Christian, N. Frank, S. Ash, T. Baumann, P. A. DeYoung, et al., Spectroscopy of Neutron-Unbound 27,28
F, Phys. Rev. C 85 (2012) 034327.
326
[13] R. C. Haight, C. Y. Wu, H. Y. Lee, T. N. Taddeucci, B. A. Perdue, et al., The LANL/LLNL Prompt
327
Fission Neutron Spectrum Program at LANSCE and Approach to Uncertainties, Nucl. Data Sheets 123
328
(2015) 130.
329
[14] C. Y. Wu, R. A. Henderson, R. C. Haight, H. Y. Lee, T. N. Taddeucci, et al., A Multiple Parallel-plate
330
Avalanche Counter for Fission-Fragment Detection, Nucl. Instrum. and Methods A 794 (2015) 76.
331
[15] H. Y. Lee, T. N. Taddeucci, R. C. Haight, T. A. Bredeweg, A. Chyzh, et al., Li-glass Detector Response 252
332
Study with a
333
703 (2013) 213.
Cf Source for Low-Energy Prompt Fission Neutrons, Nucl. Instrum. and Methods A
334
[16] https://www-nds.iaea.org/standards/.
335
[17] M. B. Chadwick, M. Herman, P. Obl˘ ozinsk´ y, M. E. Dunn, Y. Danon, et al., ENDF/B-VII.1 Nuclear
336
Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields, and Decay Data,
337
Nucl. Data Sheets 112 (2011) 2887.
338
[18] D. Neudecker, T. N. Taddeucci, R. C. Haight, H. Y. Lee, M. C. White, et al., The Need for Precise and
339
Well-documented Experimental Data on Prompt Fission Neutron Spectra from Neutron-induced Fission
340
of
239
Pu, Nucl. Data Sheets 131 (2016) 289.
341
[19] M. Marcath, mcnp Patch to Provide Initial Neutron Energies, private communication (October 2016).
342
[20] https://www-nds.iaea.org/exfor/endf.htm.
343
[21] R. Barlow, C. Beeston, Fitting using Finite Monte Carlo Samples, Comp. Phys. Commun. 77 (1993)
344
345
346
347
348
349
350
219. [22] R. Brun, F. Rademakers, ROOT – An Object Oriented Data Analysis Framework, Nucl. Instrum. Methods A 389 (1997) 81. [23] M. Q. Buckner, C. Iliadis, K. J. Kelly, L. N. Downen, A. E. Champagne, et al., Thermonuclear reaction rate of
18
O(p,γ)19 F, Phys. Rev. C 91 (2015) 015812.
[24] K. J. Kelly, A. E. Champagne, R. Longland, M. Q. Buckner, New Recommended ωγ for the Ecm = 458 r keV Resonance in
22
Ne(p,γ)23 Na, Phys. Rev. C 92 (2015) 035805.
16
351
352
353
354
[25] J. Dermigny, C. Iliadis, M. Q. Buckner, K. J. Kelly, γ-ray Spectroscopy using a Binned Likelihood Approach, Nucl. Instrum. and Methods A 830 (2016) 427. [26] K. J. Kelly, A. E. Champagne, L. N. Downen, J. R. Dermigny, S. Hunt, et al., New Measurements of Low-Energy Resonances in the
22
Ne(p,γ)23 Na Reaction, Phys. Rev. C 95 (2017) 015806.
17
PII: DOI: Reference:
S0168-9002(17)30622-8 http://dx.doi.org/10.1016/j.nima.2017.05.048 NIMA 59886
To appear in:
Nuclear Inst. and Methods in Physics Research, A
Received date : 28 March 2017 Revised date : 19 May 2017 Accepted date : 31 May 2017 Please cite this article as: K.J. Kelly, J.M. O’Donnell, J.A. Gomez, T.N. Taddeucci, M. Devlin, R.C. Haight, M.C. White, S.M. Mosby, D. Neudecker, M.Q. Buckner, C.Y. Wu, H.Y. Lee, Numerical integration of detector response functions via Monte Carlo simulations, Nuclear Inst. and Methods in Physics Research, A (2017), http://dx.doi.org/10.1016/j.nima.2017.05.048 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 proof before it is published in its final 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.
Color Compatible Version Click here to view linked References
Numerical Integration of Detector Response Functions via Monte Carlo Simulations
1
2
3 4
K.J. Kellya,1 , J.M. O’Donnella , J.A. Gomeza , T.N. Taddeuccia , M. Devlina , R.C. Haighta , M.C. Whitea , S.M. Mosbya , D. Neudeckera , M.Q. Bucknerb , C.Y. Wub , H.Y. Leea a
5
b
6
7
Los Alamos National Laboratory, Los Alamos, NM 87545, USA Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
Abstract Calculations of detector response functions are complicated because they include the intricacies of signal creation from the detector itself as well as a complex interplay between the detector, the particle-emitting target, and the entire experimental environment. As such, these functions are typically only accessible through time-consuming Monte Carlo simulations. Furthermore, the output of thousands of Monte Carlo simulations can be necessary in order to extract a physics result from a single experiment. Here we describe a method to obtain a full description of the detector response function using Monte Carlo simulations. We also show that a response function calculated in this way can be used to create Monte Carlo simulation output spectra a factor of ∼1000× faster than running a new Monte Carlo simulation. A detailed discussion of the proper treatment of uncertainties when using this and other similar methods is provided as well. This method is demonstrated and tested using simulated data from the Chi-Nu experiment, which measures prompt fission neutron spectra at the Los Alamos Neutron Science Center.
8
Keywords: Detector Response, mcnp, 6 Li-glass Detector, Neutron Detection
9
1. Introduction
10
The signal rate for a detector in an experiment is a function of the number of particles emitted from the
11
source and the response function of the detector. In a typical experiment, a particle (neutron, γ ray, ...) is
12
emitted from a source or target with an initial energy, E, and the experimental task is to detect that particle
13
and record some properties of it, such as energy, direction, etc. In some γ-ray experiments, the particle is
14
emitted from the target and potentially scatters along the way to the detector. If the γ ray reaches the
15
detector, it does so with an energy, E 0, that may or may not be equal to E. That γ ray then transfers an
16
energy, Et , into the detector which is recorded as a count in a bin of an output spectrum that corresponds
17
to Et . In the case of a γ-ray cascade, multiple γ rays may be emitted nearly simultaneously and coincidence
18
summing can be an issue. Neutron-detection experiments behave similarly with the difference that measured
19
energy, Et , of a neutron is commonly determined via a time-of-flight technique. In either type of experiment,
20
it is likely that Et 6= E or E 0. As such, determining the true distribution of initial energies, E, from a number
21
of particles measured to have a potentially different energy, Et , may be difficult and likely requires knowledge ∗
Corresponding Author Emailsubmitted address: to [email protected] (K.J. Kelly) Preprint Elsevier
May 19, 2017
22
of the detector response function. A detector response function for any kind of experiment must account for all of the above effects. In many experiments, including the experiment discussed in Section 3, experimental efficiencies of individual detectors can be low enough that the chance of multiple particles interacting in a time short enough as to create a coincidence-summing effect is negligible even in the presence of multiple particles being emitted from the target nearly simultaneously. Additionally, the particle multiplicity for the experiment discussed in Section 3 is typically below 3, which further reduces the chance for coincidence summing. Therefore, this potential aspect of detector response is not considered further in this work. As an initial attempt to describe the detector response, it is tempting to describe C(Y (E), Et ), the number of observed counts in a bin corresponding to the measured energy, Et , with an equation of the form C(Y (E), Et ) = ρ(Y (Et ), Et ) Y (Et ),
(1)
23
where Y is any yield relevant to an experiment and Y (E) and Y (Et ) are the value of Y at E and Et ,
24
respectively. This yield could be interpreted as the number of reactions per incident beam particle, the
25
number of particles emitted per unit time, or another pertinent definition. In the formalism of Eq. (1),
26
ρ(Y (Et ), Et ) is assumed to describe the response function, which includes experimental effects such as solid
27
angle. However, Eq. (1) represents a narrow view of detector response; while Eq. (1) may be able to describe
28
C(Y (E), Et ) with the proper ρ(Y (Et ), Et ), ρ depends on the detected energy, Et , as well as the initial energy
29
distribution. Thus, ρ(Y (Et ), Et ) is likely to be unique to each target or sample in an experimental setup, as
30
well as to a single target that emits particles with a different Y (E) under different experimental conditions,
31
and changes with the environment surrounding the detector as well.
32
It is more informative to write C(Y (E), Et ) in the form Z ∞ C(Y (E), Et ) = Y (E) Z0 ∞ × S(E, E 0, Et ) (E 0 ) dE 0dE,
(2)
0
33
34
35
where S is the scattering matrix (not to be confused with the nuclear physics S-matrix) describing the
probability for a particle emitted with an energy, E, to be scattered to an energy, E 0, and determined to have an energy, Et , and (E 0 ) is the detector efficiency corresponding to a particle of energy E 0. We can now define the response function R(E, Et ) as Z ∞ R(E, Et ) = S(E, E 0, Et ) (E 0 ) dE 0,
(3)
0
36
which yields C(Y (E), Et ) =
Z
0
∞
Y (E) R(E, Et ) dE.
(4)
Note that R(E, Et ) depends on E and Et but not on E 0 because the latter has been integrated out. Additionally, although the lower limit of the integral in Eq. (4) is set to 0, R(E, Et ) is essentially null below E ≈ Et . 2
It should also be noted that while Eq. (4) clearly resembles Eq. (1) of Ref. [1] and Eq. (1) of Ref. [2], it will be seen in the subsequent sections that the approach to this issue taken here is very different. Combining Eqs. (1) and (4) we see that ρ(Y (Et ), Et ) =
Z
1 Y (Et )
∞
Et
Y (E) R(E, Et )dE.
(5)
37
This shows that while ρ(Y (Et ), Et ) is related to R(E, Et ), they are fundamentally different quantities and
38
ρ(Y (Et ), Et ) is insufficient to describe the detector response. Even if a form of ρ(Y (Et ), Et ) that correctly
39
describes C(Y (E), Et ) for an energy distribution, Y (E), corresponding to a particular experimental environ-
40
ment is found, that form of ρ(Y (Et ), Et ) is no longer valid if the underlying Y (E) has changed, even if the
41
entire experimental environment remains the same.
42
2. Monte Carlo Simulations
43
In order to obtain a simulated spectrum of C(Et ) for a particular Y (E) it is common to simply run a
44
Monte Carlo simulation of the experimental environment with codes such as mcnp [3, 4] or geant4[5].
45
Running the simulation amounts to a numerical integration of Eq. (4) for the complicated response function
46
and can be time consuming for complex experiments, even with high-performance computers. The problem
47
gets much worse if a result is to be extracted from a data set by iteratively running simulations with various
48
Y (E) distributions until satisfactory agreement with the experimental data is obtained.
R
However, if R(E, Et ) is known, then the integral, i.e. Eq. (4), can be evaluated without running more simulations. Given a single simulation with high statistics obtained with E sampled from a specific initial yield distribution, Yo (E), it is possible to record the values of E and Et associated with each event in matrix form. As will be discussed in Section 3, the choice of what distribution should be used for Yo (E) is somewhat arbitrary and only the sampling rate of energies of interest needs to be considered. Whatever the chosen Yo (E) distribution, that choice is precisely defined and therefore carries no uncertainty in any subsequent calculations. Regardless of the choice of Yo (E), this matrix is exactly the product Yo (E) R(E, Et ), the integration of which over E yields the desired simulation spectrum. Furthermore, since the input Yo (E) distribution is known, it can be factored out of each simulated event, yielding the response function. It then follows that the simulation result for a different input distribution, Y 0 (E), can be obtained by simply integrating Eq. (4), instead of running a new simulation. This procedure yields the new simulation output Z ∞ C 0 (Y 0 (E), Et ) = Y 0 (E) R(E, Et ) dE. (6) Et
49
This integration can be done computationally in less than one second, as opposed to the hours or days that
50
may be required for a single Monte Carlo simulation. As is the case with any output spectrum from a
51
Monte Carlo simulation, the R(E, Et ) calculated in this way will have uncertainties related to the statistics
52
of the matrix, and also to uncertainties in the Monte Carlo simulation itself. For example, uncertainties in
53
the nuclear physics cross sections employed during the simulation and in the placement, size, and shape of 3
54
objects within the simulation can induce systematic uncertainties in R(E, Et ). While these systematic effects
55
are important for every experiment, they will not be discussed here. Even if systematic uncertainties are
56
ignored, the proper propagation of the statistical uncertainties on R(E, Et ) can be quite complicated. This
57
topic is discussed in Section 5.
58
The authors would also like to be clear in the fact that this method itself is not entirely new. Mea-
59
surements of neutron response to flight path environments are commonly carried out at facilities using an
60
incoming neutron beam [6–8], though these applications consider individual simulations when comparing to
61
measured results as opposed to creating a more widely applicable detector response matrix. Lajtai et al. [2]
62
have measured neutrons emitted from neutron-induced fission using a series of monoenergetic simulations
63
to describe their detector response. However, they never calculated a full response matrix. Instead, their
64
simulations were used to estimate their detector efficiency and to correct their data for small differences in
65
the measured time-of-flight neutron energy in an approximate way. The response matrix of fission reactor
66
cores have also been analyzed by running a series of hundreds of monoenergetic simulations to obtain a com-
67
parable description of the experimental environment [9]. Additionally, eigenmodes of reactor cores have been
68
decomposed in Ref. [10] in a similar manner. The present work is distinguished from these and other similar
69
applications primarily in that only a single, well-converged simulation was run to obtain the results shown
70
in subsequent sections, which is a fairly large simplification of the method. Finally, the work of Ref. [11]
71
and its application in Ref. [12] also represent a literature example where a single simulation is created and
72
manipulated to estimate physics parameters of an input distribution. However, the present work is unique
73
compared to those works as well in that a full description of the detector response function is provided here
74
along with crucial details discussed in Section 5 regarding the proper treatment of uncertainties when using
75
simulated spectra generated from a single R(E, Et ).
76
Although the above variations of the technique discussed in this work have been applied in certain areas
77
of physics, they are certainly not commonly employed and so the alternative derivation and description of the
78
method discussed in this work is warranted. Extensions of this method are likely also possible for analysis
79
of data from integral measurements. Furthermore, in Sections 3 and 4 we discuss the first application of
80
this particular method to a modern neutron time-of-flight experiment. Given that this method and similar
81
variants do not seem to be commonly employed in reactor or nuclear physics and are even more rarely, if ever,
82
employed in other areas of physics relying on Monte Carlo simulations, it is our hope that this publication
83
brings awareness to the time-saving potential of this technique.
84
3. Application to Simulated Data from the Chi-Nu Experiment
85
The goal of the Chi-Nu experiment at the Los Alamos Neutron Science Center (LANSCE) is to measure
86
the prompt fission neutron spectrum (pfns) for major actinides at varying incident-neutron energy ranges
87
[13]. To briefly describe the experimental setup, a pulsed, white neutron source is made incident on a parallel-
88
plate avalanche counter (ppac, see Ref. [14] for details) with target material deposited on an array of Ti foils.
4
6
Figure 1: A rendering of the Chi-Nu Li-glass detector array. Typically, the ppac target sits at the center of the array. Note R 6
that although only the Li-glass detector array is shown here, the mcnp simulation of the Chi-Nu experiment includes the floor, signal cables, and every other object included in the experimental environment.
89
Prompt fission neutrons are emitted from the target and detected with either a 6 Li-glass or a liquid scintillator
90
detector array [15]. 6 Li-glass and liquid scintillator arrays are intended to measure neutrons with E < 2 MeV
91
and E > 800 keV, respectively. Here we will focus only on the 6 Li-glass array, a rendering of which is shown
92
in Fig. 1, though treatment of the liquid array is analogous. In practice, data from 6 Li-glass detectors are
93
reliable up to ∼1 MeV because the 6 Li(n,α)t reaction cross section on which the 6 Li-glass detector efficiency
94
is based is only a standard up to 1 MeV [16]. However, the uncertainty on the 6 Li(n,α)t cross section reported
95
in ENDF/B-VII.1 [17] is below 2% up to approximately 2.5 MeV and so the simulated 6 Li-glass detector
96
spectra shown in this work are generally reported up to 2.5 MeV. Additionally, it is important to understand
97
that the pfns used in these calculations goes up much higher than this approximate upper energy limit
98
because neutrons with initial energies E > 2.5 MeV still contribute to the counts observed at Et < 2.5 MeV.
99
The experimental environment was modeled with mcnpx-PoliMi[3]. Neutron scattering, which results in
100
measured energies, Et , that are less than the true initial energies, E, adversely affects the Chi-Nu experiment
101
as well as other similar experiments. Scattering causes an excess in counts at low neutron energies and a
102
deficiency at high energies, which effectively distorts the measured detector efficiency and can distort the
103
extracted physics results as well. Significant efforts have been made to reduce this effect as much as possible
104
for Chi-Nu. During the course of characterizing the scattering effects of the Chi-Nu experiment, data from
5
105
older experiments were also discovered to be biased by uncorrected scattering as well [18]. Despite this
106
effort, data from the Chi-Nu experiment still contain a large amount of neutron scattering. This means the
107
response function must be accurately known to extract as much information out of the data as possible. The
108
application of the method described in the previous sections to data analysis at the Chi-Nu experiment has
109
significantly impacted how data are analyzed within our collaboration by reducing the amount of time it
110
takes to produce a single simulation spectrum to less than one second and reducing the weeks or months of
111
simulation CPU time required for a full data analysis to minutes.
112
In an experiment like Chi-Nu, neutrons are typically emitted from a target with a roughly Maxwellian
113
energy distribution. An example Maxwellian distribution with kT = 1.424 MeV is shown in Fig. 2a as the
114
solid line compared to literature evaluations of the neutron spectra from
115
as the dashed, dotted, and dash-dotted lines, respectively, and the ratios of these distributions to the kT =
116
1.424 MeV Maxwellian are shown in Fig. 2b to emphasize the differences between these spectra. In general,
117
the literature pfns distributions are within 5–10% of the Maxwellian spectrum. Furthermore, the initial
118
energy distribution observed during an experiment is not only a function of the target, but also of the energy
119
of the incoming neutrons that induce fission within the target. A pfns corresponding to Maxwellian with kT
120
= 1.3 MeV (as opposed to the 1.424 MeV Maxwellian shown in Fig. 2) was used as Yo (E) during integration
121
of Eq. (3) to obtain R(E, Et ) for the Chi-Nu experiment. The choice of this particular pfns was somewhat
122
arbitrary. It could have been chosen to be any distribution and the choice does not carry any additional
123
uncertainty with it since the chosen pfns is precisely defined. In principle, virtually the same R(E, Et ) can
124
be obtained regardless of the chosen Yo (E). The choice is motivated only by the desire to sample the initial
125
energies in accordance with how frequently they are encountered in an experiment, i.e. a pfns in the case
126
of Chi-Nu. A patch was made to mcnpx-PoliMi by M. Marcath [19] in order to obtain the initial neutron
127
energies for each simulated event so that Yo (E) could be factored out of Eq. (4). Following the steps described
128
in Section 2, the Maxwellian pfns used during the mcnp simulation was factored out of each event, yielding
129
R(E, Et ) for the 6 Li-glass array.
252
Cf,
235
U, and
239
Pu [17] shown
130
The resulting response function is shown in Fig. 3. A few comments are in order regarding the structure
131
of this response function. First, the dark band along E = Et corresponds to neutrons that did not scatter
132
along the path to the detector. This band is broadened according to the known detector resolution and light
133
curves, which allows for some counts to be accumulated with E < Et . The next most prominent feature
134
is the broad band of counts around E = 240 keV. These counts correspond to the 6 Li(n, α)t resonance at
135
approximately the same energy, with the lower Et counts in this band arising from neutrons emitted with the
136
resonance energy that then scattered one or more times along the path to the detector. The same statement
137
is true of all counts observed to be measured with energies Et < E, i.e. these are neutrons that were affected
138
by multiple scattering before detection. It is remarkable that this level of multiple scattering can still be seen
139
in data from the Chi-Nu experiment given that the experiment and the experimental environment were all
140
designed with an effort towards reducing the amount of multiple scattering in the data. Finally, it should be
6
0.3
PFNS (1/MeV)
PFNS Ratio to Maxwellian (kT = 1.424 MeV)
Maxwellian (kT = 1.424 MeV) 252 Cf (ENDF/B-VII.1) 235 U (ENDF/B-VII.1) 239 Pu (ENDF/B-VII.1)
0.35
0.25
0.2
0.15
0.1
0.05 10
−2
10
−1
Maxwellian (kT = 1.424 MeV) 252 Cf (ENDF/B-VII.1) 235 U (ENDF/B-VII.1) 239 Pu (ENDF/B-VII.1)
1.06 1.04 1.02 1 0.98 0.96 0.94 0.92 0.9 10
1
Initial Outgoing Neutron Energy (MeV)
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
(a) Examples of pfns distributions from commonly-studied actinides compared to a Maxwellian distribution with kT = 1.424 MeV.
(b) The same pfns distributions as shown in panel (a), but shown as a ratio to a Maxwellian distribution with kT = 1.424 MeV
Figure 2: Panels (a) showa a comparison of common pfns distributions and panel (b) shows the same distributions as a ratio to a Maxwellian of kT = 1.424 MeV. The solid, dashed, dotted, and dash-dotted lines are a Maxwellian with kT = 1.424 MeV, and 252 235 239 235 239 ENDF/B-VII.1 pfns distributions for Cf, U, and Pu, respectively. The U and Pu pfns distributions correspond to an incoming neutron energy of 1.0 MeV. These distributions are red, black, blue, and green online, respectively.
141
emphasized that R(E, Et ) depends on the entire environment surrounding the experimental area, in addition
142
to the detectors themselves and even to the target material under study. While it is likely true that, for
143
Chi-Nu, a change in target material from, for example,
144
the detector response matrix, the change is not zero. Therefore, even identical ppac’s with different target
145
materials require separately calculated response matrices. This environmental dependence can be seen in
146
Fig. 3. The thin bands seen at E ≈ 3.0, 2.5, 1.2, and 0.5 MeV and other energies correspond to resonances in
147
aluminum, oxygen, etc., as well as to objects in the experimental environment off which neutrons are likely to
148
scatter, thereby creating bands at specific times of flight [20]. Any change in the experimental environment
149
surrounding the Chi-Nu detector array will induce a change in R(E, Et ).
235
U to
239
Pu does not cause a significant change in
150
Given this R(E, Et ) obtained from a single, well-converged simulation, new Monte Carlo spectra were
151
created via Eq. (6) using various pfns distributions, each of which required time on the order of milliseconds
152
to create, and tests were carried out to confirm the accuracy of the new simulation output. The results of
153
these tests are discussed in the following section.
154
4. Results
155
Figure 4a shows some example Maxwellian pfns distributions at temperatures ranging from kT = 1.1-
156
1.4 MeV and Fig. 4b shows the corresponding spectra created using the R(E, Et ) shown in Eq. (6) and
157
Fig. 3. Note that the large peak at ∼240 keV present in each spectrum shown in Fig. 4 is caused by the
158
resonance in the 6 Li(n,α)t reaction at the same energy. These results demonstrate our ability to produce
159
simulation output for any desired input pfns. As a test of these and other new simulation spectra, it is
160
natural to simply compare simulation output created via the method described in this work to that of an
7
Initial Outgoing Neutron Energy, E (MeV)
10
103 1 102
10−1 10
10−2
10−1
1
10
1
Outgoing Neutron Energy from Time of Flight, E t (MeV)
6
Figure 3: The detector response function, R(E, Et ), for the Chi-Nu experimental setup using a Li-glass detector array. This 9 response function has not been normalized in any way and was calculated with a total of 1.98×10 neutrons emitted evenly over the ten targets within a single Chi-Nu ppac. The dashed line along E = Et shown to guide the eye. See the text for a discussion of the main features.
235
161
independently-calculated simulation. An independent simulation of a literature
U pfns [17] was carried
162
out for this comparison. The resulting simulation spectra from mcnp and from the convolution of the
163
pfns with R(E, Et ) are shown in Fig. 5a and the ratio of these spectra is shown in Fig. 5b. Both spectra
164
appear to be in excellent agreement. Although it may be noted that the statistical fluctuations become quite
165
large at lower time-of-flight energies where fewer counts are observed, the statistical uncertainty on these bins
166
is approximately 5%. This uncertainty aligns well with the variation about unity observed in Fig. 5b. It is
167
also worth noting that R(E, Et ) and the simulation shown in Fig. 5a were calculated using different random
168
number seeds in mcnpx-PoliMi.
235
U
169
A further test of the accuracy of the simulation output produced via Eq. (6) is to fit these output spectra
170
and use R(E, Et ) to extract the pfns used to create the calculated spectrum. Fits described here utilize
171
the TFractionFitter [21] class of root [22]. Detailed descriptions of this fitting routine have been provided
172
in previous works (see Refs. [23–26]), and so the description here will be kept brief. In short, a spectrum
173
is produced from R(E, Et ) for narrow ranges of initial neutron energy, E (i.e. an Et spectrum is produced
174
for each slice of the R(E, Et ) shown in Fig. 3 on a 10 bins per decade basis, in this case). Collectively,
175
these spectra span all E, are referred to as templates, and describe the portion of the detector response
176
corresponding to that particular E range. The goal is to determine the pfns used during creation of the
177
simulation output, which amounts to determining the relative contribution of each template spectrum to the
178
simulated data. The TFractionFitter class accomplishes this goal via a log-likelihood maximization routine
179
considering statistical variations in both the template spectra and in the data (real or simulated) being fit.
180
It is important to note that these fits do not assume anything about the shape of the pfns used to create the
181
simulation output. In this way these fits are a stringent test of the simulation output calculated via Eq. (6).
182
The resulting fit to simulation output created via a Maxwellian pfns of kT = 1.1 MeV is shown in Fig. 6a
183
with the input pfns as the solid line and the pfns from the fit as the diamonds. The same pfns and fit result
8
20000
Maxwellian kT = 1.1 MeV
18000
Maxwellian kT = 1.2 MeV
Maxwellian kT = 1.1 MeV Maxwellian kT = 1.2 MeV
0.4
Maxwellian kT = 1.3 MeV
16000
Maxwellian kT = 1.4 MeV
14000
Maxwellian kT = 1.3 MeV Maxwellian kT = 1.4 MeV
Counts
PFNS (1/MeV)
0.5
0.3
12000 10000 8000
0.2
6000 4000
0.1
2000
10
−1
10
1
−1
1
Outgoing Neutron Energy from Time of Flight (MeV)
Initial Outgoing Neutron Energy (MeV)
(b) A series of histograms created using the R(E, Et ) shown in Fig. 3 with the series of Maxwellian distributions on panel (a). The full circle, open diamond, closed diamond, and open circle data points correspond to temperatures of kT = 1.1, 1.2, 1.3, and 1.4 MeV (black, red, green, and blue online, respectively).
(a) A series of Maxwellian pfns distributions. The solid, dashed, dotted, and dash-dotted distributions correspond to temperatures of kT = 1.1, 1.2, 1.3, and 1.4 MeV (black, red, green, and blue, online, respectively).
Figure 4: The series of Maxwellian pfns distributions shown in panel (a) were used with the R(E, Et ) shown in Fig. 3 via Eq. (6) to create the histograms shown in panel (b). No normalization was applied to the simulated spectra because they were all created with the same R(E, Et ).
1.2
12000
Independent MCNP Simulation
1.15
Spectrum Calculated via Eq. (6) 6
Ratio of Simulation Spectra
14000
Li(n ,α )t Cross Section Trend
Counts
10000 8000 6000 4000 2000
1.1 1.05 1 0.95 0.9 0.85
10−1
0.8 1
Outgoing Neutron Energy from Time of Flight (MeV)
10−1
1
Outgoing Neutron Energy from Time of Flight (MeV)
(a) A comparison of an independent mcnp simulation to the simulation output calculated via the method described in this work.
(b) The ratio of spectra shown in panel (a). A dashed line is shown at unity for reference.
Figure 5: Panel (a) shows a comparison of an independent mcnp simulation spectrum obtained using an ENDF/B-VII.1 [17] 235 U pfns corresponding to an incoming neutron energy of 1.0 MeV shown as the solid line (black online, calculated with 8 9.9×10 neutrons emitted evenly over the ten targets within a single Chi-Nu ppac) to a spectrum calculated with Eq. (6) using 235 the same U input pfns convolved with R(E, Et ) shown as the open circle data points (red online). Excellent visual agreement 6 across the entire data spectrum. The Li(n,α)t cross section trend from ENDF/B-VII.1 [17] is shown as the dashed (blue online) 6 line as a reference for the structure seen at ∼240 keV in the Li-glass spectra. The ratio of the spectra shown in panel (a) is displayed in panel (b). See the text for a discussion of this panel.
9
Simulated Maxwellian (kT=1.10 MeV) PFNS from Fit
0.4
PFNS (1/MeV)
PFNS Relative To Maxwellian (kT=1.3 MeV)
0.45
0.35 0.3 0.25 0.2 0.15 0.1
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
1.6
Simulated Maxwellian (kT=1.10 MeV) PFNS from Fit
1.4
1.2
1
0.8
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
(a) Results of a fit to simulated data created via Eq. (6) using a Maxwellian of kT = 1.1 MeV.
(b) Results of a fit to simulated data created via Eq. (6) using a Maxwellian of kT = 1.1 MeV shown as a ratio to a Maxwellian of kT = 1.3 MeV.
Figure 6: A Maxwellian pfns with kT = 1.1 MeV is shown in panel (a) as the solid line (blue online). A simulation spectrum was created via Eq. (6) using this pfns and a fit was carried out to extract the pfns from the simulation data. The result of the fit is shown as the diamond data points (red online) in the same panel. The simulated Maxwellian pfns and fit result are shown in panel (b) as a ratio to a Maxwellian of kT = 1.3 MeV.
184
are shown in Fig. 6b as a ratio to a kT = 1.3 MeV Maxwellian to emphasize the details of this fit. The pfns
185
from the fit is nearly identical to the simulated Maxwellian distribution. An analogous test was also carried
186
out with a drastically different input pfns, a Maxwellian of kT = 1.65 MeV. Figure 7a shows the kT = 1.65
187
MeV Maxwellian distribution as the solid line with the extracted pfns from the fit as the diamonds and the
188
result is again shown as a ratio to a kT = 1.3 MeV Maxwellian in Fig. 7b. The kT = 1.1 MeV Maxwellian
189
from Fig. 6 is reproduced as the dashed line in Fig. 7 as a reference for the large change from one pfns to
190
the other. Again, a nearly identical reproduction of the simulated pfns was obtained from the fit.
191
Given that the fit results discussed in Figs. 6 and 7 as well as the simulation output data themselves
192
were obtained with the same R(E, Et ), it is not entirely surprising that a nearly exact reproduction of the
193
input pfns can be extracted. So, as a final test of this method, simulation data calculated directly with
194
mcnpx-PoliMi using a literature
195
and with different random seeds were fit as well. This input pfns is shown as the solid line and the fit pfns
196
is shown as the diamonds in Fig. 8a. The same results are shown in Fig. 8b relative to a Maxwellian pfns
197
of kT = 1.3 MeV to emphasize the details of Fig. 8a. The Maxwellian pfns distributions that were fit in
198
Figs. 6a and 7a are reproduced in Fig. 8 as the dashed and dash-dotted lines for reference. Clearly, there is
199
more variation in the resulting fit pfns of Fig. 8a than in Figs. 6 or 7.
200
235
U pfns [17] corresponding to an incoming neutron energy of 1.0 MeV
Although an increase in simulation time for both R(E, Et ) and the
235
U mcnp simulation would likely
201
improve the fit results shown in Figs. 8a and 8b because of the resulting smaller statistical uncertainties in
202
percent of the total counts in each bin, the variations in the extracted pfns are primarily caused by the
203
fitting routine itself. For the purposes of this work, the TFractionFitter class is an unfolding method [13]
204
in the sense that the input energy distribution, i.e. the pfns, is extracted without assuming any shape of
10
PFNS Relative To Maxwellian (kT=1.3 MeV)
Reference Maxwellian (kT=1.10 MeV)
0.45
Simulated Maxwellian (kT=1.65 MeV) PFNS from Fit
PFNS (1/MeV)
0.4 0.35 0.3 0.25 0.2 0.15 0.1
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
1.6
Reference Maxwellian (kT=1.10 MeV) Simulated Maxwellian (kT=1.65 MeV) PFNS from Fit
1.4
1.2
1
0.8
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
(a) Results of a fit to simulated data created via Eq. (6) using a Maxwellian of kT = 1.65 MeV. The kT = 1.1 MeV Maxwellian from Fig. 6a is shown for reference.
(b) Results of a fit to simulated data created via Eq. (6) using a Maxwellian of kT = 1.1 MeV shown as a ratio to a Maxwellian of kT = 1.3 MeV. The kT = 1.1 MeV Maxwellian from Fig. 6b is shown for reference.
Figure 7: A Maxwellian pfns with kT = 1.65 MeV is shown in panel (a) as the solid line (green online). A simulation spectrum was created via Eq. (6) using this pfns and a fit was carried out to extract the pfns from the simulation data. The result of the fit is shown as the diamond data points (red online) in the same panel. The simulated Maxwellian pfns and fit result are shown in panel (b) as a ratio to a Maxwellian of kT = 1.3 MeV. The kT = 1.1 MeV Maxwellian from Figs. 6a and 6b is reproduced here as the dashed line (blue online) for reference.
205
that distribution and using only knowledge of the detector response function. Unfolding methods have the
206
inherent property that statistical noise, however small it may appear in the measured or simulated data, is
207
amplified into an oscillatory variation in the resulting fit. This is most easily seen below 70 keV in Figs. 8a
208
and 8b. In the fitting performed to obtain the results shown in Figs. 8a and 8b, the small changes between
209
the 235 U simulated spectrum and R(E, Et ) caused by the different choice of random number seed are enough
210
to cause these strong variations in the fit result. The magnitude and character of these variations may be
211
different for other unfolding techniques, such as single value decomposition or Bayesian unfolding, but the
212
variations themselves will likely be present in any unfolding result. Furthermore, these variations can be
213
reduced by enforcing a smoothing condition on the unfolding, but this enforcement is somewhat artificial and
214
would likely make the results discussed in this work less clear. As such, no such smoothing was applied to
215
any fit described in this work. Despite these variations in the fits, it is still clear that the correct input pfns
216
has been extracted from this final test, thereby verifying that the correct mcnp spectrum was produced by
217
the method discussed in this work.
218
5. Notes on Uncertainties Associated with the Monte Carlo Simulation Output
219
When an independent simulation is run, the output spectrum can be assumed to have a Poisson distributed
220
statistical uncertainty, and so the statistical uncertainty on each bin can generally be taken to simply be the
221
square root of the number of counts. Given that the output of an independent simulation spectrum can be
222
nearly exactly reproduced using the R(E, Et ) shown in Fig. 3 and Eq. (6), it may seem natural to assume
223
that the statistical uncertainty on each bin of the spectrum calculated via Eq. (6) is also the square root of the 11
PFNS Relative To Maxwellian (kT=1.3 MeV)
Reference Maxwellian (kT=1.10 MeV)
0.45
Reference Maxwellian (kT=1.65 MeV) 235
0.4
U PFNS (ENDF/B-VII.1)
PFNS (1/MeV)
PFNS from Fit
0.35 0.3 0.25 0.2 0.15 0.1
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
Reference Maxwellian (kT=1.10 MeV)
1.6
Reference Maxwellian (kT=1.65 MeV) 235
U PFNS (ENDF/B-VII.1)
PFNS from Fit
1.4
1.2
1
0.8
10
−2
10
−1
1
Initial Outgoing Neutron Energy (MeV)
(a) Results of a fit to simulated data created via Eq. (6) 235 using a U pfns from ENDF/B-VII.1 [17].
(b) Results of a fit to simulated data created via Eq. (6) 235 using a U pfns from ENDF/B-VII.1 [17]. shown as a ratio to a Maxwellian with kT = 1.3 MeV.
235
Figure 8: (a) A literature U pfns is shown as the solid line (black online) and the extracted pfns from the mcnp spectrum using different initial random number seeds than in Fig. 7a is shown as the diamonds (red online). The Maxwellian distributions used in Figs. 6 and 7 are shown as the dashed and dash-dotted lines (blue and green online) for reference, respectively. (b) A 235 reproduction of Fig. 8a shown relative to a Maxwellian of kT = 1.3 MeV in order to emphasize variations of the U pfns from the fit.
q q C(Y 0 (E), Et ). However, the simple C(Y 0 (E), Et ) statistical uncertainty
224
number of counts in that bin,
225
assumption does not apply to spectra calculated via Eq. (6) because all simulation output spectra calculated
226
using R(E, Et ) are correlated through the use of a single simulation. It is important to understand that only
229
the original, independently-calculated simulation spectrum used to obtain R(E, Et ) maintains the simple q C(Y 0 (E), Et ) uncertainty relationship.
230
Chi-Nu experiment. This pfns has no uncertainty associated with it because it is a somewhat arbitrary,
231
but precisely-defined choice for a probability density function from which initial energies are sampled and
232
that choice is motivated only by the requirement that initial energies are sampled in accordance with their
233
importance in the approximate expectation from experimental data, i.e. a pfns distribution. So, if we ignore
234
systematic uncertainties in the Monte Carlo calculation, the only source of uncertainty in spectra calculated
235
via Eq. (6) is the statistical uncertainty number of counts in each (E,Et ) bin of the original simulation
236
calculated with Yo (E), NE,Et . The associated covariance matrix for NE,Et is described by
227
228
As was mentioned earlier, an arbitrary Maxwellian pfns was used as Yo (E) to calculate R(E, Et ) for the
NEi ,Et
= σ(NE ,E )σ(N 0 0 Cov Ej ,Et )δEi Ej δEt Et , i t NEj ,Et0
(7)
237
where the δ’s in Eq. (7) are Kronecker deltas and σ 2 (NEi ,Et ) is the variance of NEi ,Et . Calculation of a
238
new simulation spectrum corresponding to any input pfns distribution, Yα (E), then amounts to scaling the
239
contents of each bin in R(E, Et ) by the ratio Yα (E)/Yo (E) and summing over all initial energies, E, that
12
240
contribute to each measured energy bin, Et , i.e. C(Yα (E), Et ) = =
Z
∞
0 n X i=1
Yα (E) R(E, Et ) dE
Yα (Ei ) N . Yo (Ei ) Ei ,Et
(8)
241
Using Eqs. (7) and (8), it can be shown that the covariance between C(Yα (E), Et ), C(Yβ (E), Et ), ...,
242
C(Yζ (E), Et ), i.e. the number of counts in the bin of a Monte Carlo simulation output spectrum corre-
243
sponding to Et calculated with Eq. (6) using various pfns distributions, Yα (E), Yβ (E), ..., Yζ (E), in place
244
of Y 0 (E), is given by the matrix σ2 C(Yα (E), Et ) α C(Yβ (E), Et ) σβα = Cov .. .. . . σζα
C(Yζ (E), Et )
with
σαβ =
n X Yα (Ei )Yβ (Ei ) i=1
Yo2 (Ei )
σαβ σβ2
··· ..
σζβ
.
···
σαζ
σβζ .. , . σζ2
(9)
σ 2 (NEi ,Et ).
(10)
245
The remaining elements of the above covariance matrix can be easily obtained by replacing the indices α and
246
β with the appropriate index in Eq. (10) keeping in mind that σαα = σα2 . Note that all pfns distributions
247
in the above equations are evaluated at the initial neutron energy, E, and not at Et . Also note that this
248
covariance matrix could be simplified by choosing a flat distribution instead of the chosen Yo (E). This may
249
be convenient for certain applications, but if the expected initial distribution is far from uniform, then this
250
simplification would be made at the cost of running the initial simulation for a much longer time in order to
251
obtain sufficient statistics at the important energies in R(E, Et ), as was the case with the Chi-Nu experiment
252
discussed in Sections 3 and 4.
253
The above covariance matrix assumes that each Et bin is independent of each other bin, Et0 . However, if
254
the desired new input distribution, Yµ (E), is defined by some function and the parameters describing that
255
function have uncertainties associated with them, then there may well be a non-zero covariance between each
256
value of Et . This is frequently the case for pfns distributions described by a functional form such as, for
257
example, a Maxwellian or a Watt distribution. In these distributions, changing a single parameter can change
258
the entire pfns across all initial energies and, in turn, across all measured energies, Et . For simplicity, we
259
will consider only two E bins, E1 and E2 , and two Et bins, Et and Et0 . If the covariance matrix for the
260
function describing the initial energy distribution Yµ (E) is given by
261
262
Yµ (E1 ) σµ2 (E1 ) σµ (E1 , E2 ) σ2 = = µ1 Cov Yµ (E2 ) σµ (E2 , E1 ) σµ2 (E2 ) σµ21
σµ12 2 σµ2
,
(11)
with all elements greater than zero, then the covariance between bins Et and Et0 of the simulation output 13
263
spectrum is given by C(Yµ (E), Et ) σ 2 (Et ) σ(Et , Et0 ) = . Cov C(Yµ (E), Et0 ) σ(Et0 , Et ) σ 2 (Et0 )
264
σ 2 (Et ) σ(Et , Et0 )
2 2 2 X X N N Y (E ) Ei ,Et Ej ,Et µ i σ 2 (NE ,E )+ = σ i t Y (E ) Y (Ei )Yo (Ej ) µij o i i=1 j=1 o
= σ(Et0 , Et ) =
2 X 2 N X Ei ,Et NEj ,Et0 i=1 j=1
Yo (Ei )Yo (Ej )
σµij
(12)
(13)
(14)
265
Note that σ 2 (Et0 ) can be obtained by replacing Et with Et0 in Eq. (13). The first term in the diagonal com-
266
ponents of this matrix corresponds to the uncertainty associated with NEi ,Et and NEi ,Et0 , but the remaining
267
diagonal terms and the off-diagonal terms arise from the covariance between the simulation output bins.
268
These additional components impact the assumed uncertainty in the resulting simulation output spectrum
269
and must be propagated through to the final covariance matrix.
270
It should be noted that although uncertainty propagation with simulation output spectra calculated via
271
Eq. (6) is significantly more complicated than a simple Gaussian or Poisson statistical approach, a much longer
272
simulation run time can be, but is not required to be, chosen for the R(E, Et ) calculation because there is
273
only one simulation to be run. This guarantees that the statistical uncertainty on the counts in each bin of
274
an output spectrum calculated via Eq. (6) is lower than would be obtained for each of potentially hundreds or
275
thousands of independent Monte Carlo simulations. Additionally, although the resulting covariance matrices
276
are not shown here, the uncertainties obtained with this method are even further reduced when considering the
277
multiplication or division of multiple simulation output spectra. In these cases, the covariance of the result
278
can be significantly lower than that obtained by multiplying or dividing simulation output spectra from
279
multiple independent calculations. Clearly, the covariances that must be considered in order to properly
280
handle the uncertainties on the simulation output produced via the method described in this work can be
281
quite complicated. Nonetheless, analytic expressions for each element of the necessary covariance matrices
282
can be derived, and thus the problem is tractable.
283
6. Conclusions and Future Work
284
The results shown in Section 4 demonstrate that the output of a Monte Carlo simulation for any input
285
distribution can be obtained by calculating the detector response function just once with a single simulation
286
and incorporating the desired input energy distribution in a subsequent integration of the detector response.
287
Furthermore, these output spectra can be calculated a factor of ∼1000× faster than running a new simulation.
288
The amount of time saved with this method will have a profound impact on any experiment that relies on
289
Monte Carlo simulations in which the input distribution is varied, but the experimental geometry remains
290
the same. While the uncertainties on the simulation output spectra created with this method are much
14
291
more complicated than a simple Gaussian or Poisson statistics uncertainty on each bin, the uncertainty on
292
each bin can actually be significantly smaller than would be obtained by running hundreds or thousands of
293
independent simulations. Despite the fact that variants of this method appear to be used primarily, though
294
not commonly, in reactor physics at the present time, its application can significantly impact research in
295
nuclear physics, nuclear astrophysics, and many other fields.
296
7. Acnknowledgments
297
298
The authors would like to thank M.E. Rising and F.B. Brown for fruitful discussions, M. Marcath for assistance with mcnp-related issues, and also P. Talou for helpful comments on this manuscript.
299
[1] T. N. Taddeucci, R. C. Haight, H. Y. Lee, D. Neudecker, J. M. O’Donnell, et al., Multiple-scattering
300
Corrections to Measurements of the Prompt Fission Neutron Spectrum, Nucl. Data Sheets 1232 (2015)
301
135.
302
303
304
305
306
307
308
309
310
311
[2] A. Lajtai, P. P. Dyachenko, V. N. Kononov, E. A. Seregina, Low-Energy Neutron Spectrometer and its Application for
252
Cf Neutron Spectrum Measurements, Nucl. Instrum. and Methods A 293 (1990) 555.
[3] S. A. Pozzi, E. Padovani, M. Marsequerra, A Monte-Carlo Code for Correlation Measurements, Nucl. Instrum. and Methods A 513 (2003) 550. [4] S. A. Pozzi, S. D. Clarke, W. J. Walsh, E. C. Miller, J. L. Dolan, et al., MCNPX-PoliMi for Nuclear Nonproliferation Applications, Nucl. Instrum. and Methods A 694 (2012) 119. [5] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo, et al., Geant4-A Simulation Toolkit, Nucl. Instrum. and Methods A 506 (2003) 250. [6] C. Coceva, R. Simonini, Calculation of the ORELA Neutron Moderator Spectrum and Resolution Function, Nucl. Instrum. and Methods A 211 (1983) 459.
312
[7] M. Mocko, G. Muhrer, T. Ino, M. Ooi, L. L. Daemen, et al., Monte Carlo study of the neutron time
313
emission spectra at the Manuel Lujan Jr. Neutron Scattering Center, Nucl. Instrum. and Methods A
314
594 (2008) 373.
315
316
317
318
319
320
[8] C. Guerrero, A. Tsinganis, E. Berthoumieux, M. Barbagallo, F. Belloni, et al., Performance of the Neutron Time-of-Flight Facility n TOF at CERN, Eur. Phys. J. A 49 (2013) 27. [9] R. G. Gamino, F. B. Brown, M. R. Mendelson, A Monte Carlo Green’s Function Method for ThreeDimensional Neutron Transport, Trans. Am. Nucl. Soc. 65 (1992) 237. [10] F. B. Brown, S. E. Carney, B. C. Kiedrowski, W. R. Martin, Fission Matrix Capability for MCNP Monte Carlo, Los Alamos National Laboratory Report LA-UR-13-26962.
15
321
[11] D. M. Schmidt, R. J. Morrison, M. S. Witherell, A General Method of Estimating Physics Parameters
322
from a Distribution with Acceptance and Smearing Effects, Nucl. Instrum. and Methods A 328 (1993)
323
547.
324
325
[12] G. Christian, N. Frank, S. Ash, T. Baumann, P. A. DeYoung, et al., Spectroscopy of Neutron-Unbound 27,28
F, Phys. Rev. C 85 (2012) 034327.
326
[13] R. C. Haight, C. Y. Wu, H. Y. Lee, T. N. Taddeucci, B. A. Perdue, et al., The LANL/LLNL Prompt
327
Fission Neutron Spectrum Program at LANSCE and Approach to Uncertainties, Nucl. Data Sheets 123
328
(2015) 130.
329
[14] C. Y. Wu, R. A. Henderson, R. C. Haight, H. Y. Lee, T. N. Taddeucci, et al., A Multiple Parallel-plate
330
Avalanche Counter for Fission-Fragment Detection, Nucl. Instrum. and Methods A 794 (2015) 76.
331
[15] H. Y. Lee, T. N. Taddeucci, R. C. Haight, T. A. Bredeweg, A. Chyzh, et al., Li-glass Detector Response 252
332
Study with a
333
703 (2013) 213.
Cf Source for Low-Energy Prompt Fission Neutrons, Nucl. Instrum. and Methods A
334
[16] https://www-nds.iaea.org/standards/.
335
[17] M. B. Chadwick, M. Herman, P. Obl˘ ozinsk´ y, M. E. Dunn, Y. Danon, et al., ENDF/B-VII.1 Nuclear
336
Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields, and Decay Data,
337
Nucl. Data Sheets 112 (2011) 2887.
338
[18] D. Neudecker, T. N. Taddeucci, R. C. Haight, H. Y. Lee, M. C. White, et al., The Need for Precise and
339
Well-documented Experimental Data on Prompt Fission Neutron Spectra from Neutron-induced Fission
340
of
239
Pu, Nucl. Data Sheets 131 (2016) 289.
341
[19] M. Marcath, mcnp Patch to Provide Initial Neutron Energies, private communication (October 2016).
342
[20] https://www-nds.iaea.org/exfor/endf.htm.
343
[21] R. Barlow, C. Beeston, Fitting using Finite Monte Carlo Samples, Comp. Phys. Commun. 77 (1993)
344
345
346
347
348
349
350
219. [22] R. Brun, F. Rademakers, ROOT – An Object Oriented Data Analysis Framework, Nucl. Instrum. Methods A 389 (1997) 81. [23] M. Q. Buckner, C. Iliadis, K. J. Kelly, L. N. Downen, A. E. Champagne, et al., Thermonuclear reaction rate of
18
O(p,γ)19 F, Phys. Rev. C 91 (2015) 015812.
[24] K. J. Kelly, A. E. Champagne, R. Longland, M. Q. Buckner, New Recommended ωγ for the Ecm = 458 r keV Resonance in
22
Ne(p,γ)23 Na, Phys. Rev. C 92 (2015) 035805.
16
351
352
353
354
[25] J. Dermigny, C. Iliadis, M. Q. Buckner, K. J. Kelly, γ-ray Spectroscopy using a Binned Likelihood Approach, Nucl. Instrum. and Methods A 830 (2016) 427. [26] K. J. Kelly, A. E. Champagne, L. N. Downen, J. R. Dermigny, S. Hunt, et al., New Measurements of Low-Energy Resonances in the
22
Ne(p,γ)23 Na Reaction, Phys. Rev. C 95 (2017) 015806.
17