ARTICLE IN PRESS
Nuclear Instruments and Methods in Physics Research A 578 (2007) 261–266 www.elsevier.com/locate/nima
The empirical characterization of organic liquid scintillation detectors by the normalized average of digitized pulse shapes M.D. Aspinalla,, B. D’Mellowa, R.O. Mackina, M.J. Joycea, Z. Jarrahb, A.J. Peytonc a
Department of Engineering, Lancaster University, Bailrigg, Lancaster, Lancashire, LA1 4YR, UK b Hybrid instruments Ltd., Priory Close, St. Mary’s Gate, Lancaster, Lancashire, LA1 4WA, UK c School of Electrical and Electronic Engineering, University of Manchester, P.O. Box 88, Manchester, M60 1QD, UK Received 1 May 2007; accepted 7 May 2007 Available online 16 May 2007
Abstract The application of the digital acquisition of pulses from a liquid scintillator to detector characterization is described. Experimental data for a mixed neutron/g-ray field have been recorded digitally. An empirical method for the characterization of liquid scintillation detectors, in terms of their pulse shape, has been developed which is quick and easy. It provides generic pulses shapes for use in pulse-shape discrimination and that can be used to derive analytical descriptions of each pulse via an accepted six-parameter formulism. The distributions of the neutron and g-ray components arising as a result of discrimination via pulse gradient analysis (PGA) follow a bi-Gaussian trend and exhibit degrees of both kurtosis and skew. r 2007 Elsevier B.V. All rights reserved. PACS: 28.20.v; 28.41.Rc; 29.25.Dz; 29.30.Hs; 29.40.Mc; 29.85.+c; 29.85+v; 07.05.Hd Keywords: Liquid scintillators; Digital; Pulse-shape discrimination
1. Introduction Recently, a fast and efficient method for the digital discrimination of neutrons and g rays in mixed radiation fields was reported [1] known as pulse gradient analysis (PGA). In addition, a six-parameter formulism for the generic mathematical description of pulses from liquid scintillators has been reported [2] which can be used to simulate responses from such detectors. However, the combination of these approaches, such that synthetic pulses can be obtained easily for any of the wide range of liquid scintillator pulse volumes that are possible, has not been straightforward. In this paper, we present an empirical method to obtain a generic pulse from a given detector volume which can be used to extract the six parameters necessary for the production of synthetic pulses. Corresponding author. Tel.: +44 1524 593326, +44 1524 60537.
E-mail address:
[email protected] (M.D. Aspinall). 0168-9002/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2007.05.114
Liquid scintillators have been used for many years as neutron detectors because of their excellent timing performance and good pulse-shape discrimination (PSD) properties. Numerous PSD techniques have been developed to produce clean neutron signatures from mixed neutron and g-ray radiation fields with applications in neutron monitors and neutron recoil spectrometers [3–8]. PSD in liquid scintillators relies on the long lifetime components of the scintillation light produced by high LET events such as proton recoils caused by neutron scatter events. This phenomenon has been extensively studied over the last 40 years, and an assortment of analogue pulse discrimination techniques have been used to extract neutron and g-ray information. Recent advances in highspeed digital electronics have opened a new window of investigation, and studies into digital PSD techniques are now being reported. PSD algorithms benefit in the digital domain from improved neutron/g discrimination compared to analogue techniques. The simple PSD algorithm employed in this study is inherently fast, making it suitable
ARTICLE IN PRESS M.D. Aspinall et al. / Nuclear Instruments and Methods in Physics Research A 578 (2007) 261–266
for the implementation of fast real-time pulse processing in field instruments. In this paper, we present the results of the PSD algorithm with supporting statistical analysis. From our discriminated data sets, average pulse shapes for g-ray and neutron events have been obtained. These can be used as a basis on which to discriminate neutron and g-ray events. Where required, the six-parameter formulism can be applied to each generic pulse shape to extract an analytical description of each pulse shape arising from a given detector. 2. Experimental details The experimental setup used in this research consisted of a 4.5 ml cylindrical cell of EJ-301 liquid scintillator (John Caunt Sci., UK), optically coupled to a Hamamatsu R5611 photomultiplier tube (PMT). The PMT was operated with a negative high-tension (HT) supply voltage of 840 V dc. An 8 m length of 50 O coaxial cable was used to transport the signal from the PMT to an Infiniiums digital oscilloscope (Agilent Tech.). Data were recorded digitally at 4 GHz sampling frequency and 12 bit amplitude resolution. Data acquisition was automated by driving the oscilloscope remotely via a Transmission Control Protocol/Internet Protocol (TCP/IP) connection to a personal computer running a bespoke MATLABs m-file program. The personal computer used for this system was a DELLs Latitude D400 laptop. Events were acquired from an americium–beryllium (Am–Be) neutron source suspended in a moderating water tank and shielded externally on three sides with 30 mm of plate cadmium encased in sheet plywood, as shown schematically in Fig. 1. This source can be moved forwards and backwards to vary the degree of neutron moderation. The source was positioned 180 mm from the outside edge of the front of the tank, and the scintillator was positioned in line with the source (360 mm above ground level) at a distance of 800 mm from the front edge of the tank. A total of 5000 events were collected and the amplitude data were saved to the computer’s hard drive.
Fig. 1. An americium–beryllium (Am–Be) neutron source suspended in a moderating water tank and externally shielded on three sides with 30 mm of plate cadmium encased in sheet plywood. The source is positioned 180 mm from the outside edge of the front of the tank, the scintillator was positioned in line with the source 800 mm from the front edge of the tank.
0
x 104
–0.5 –1 Amplitude
262
–1.5 –2 –2.5 –3 100
120
140
160
180
200
Time, ns Fig. 2. A plot of amplitude versus time for 20 superimposed (50:50) mixed neutron/g pulses acquired from an Am–Be source using an EJ-301 liquid scintillator shielded with a 50 mm thick lead brick wall, at 4 GHz sampling frequency (400 samples per pulse window).
The low neutron/g ratio for this particular source and high contribution of g rays from scattering results in a flux constituted mainly of g-ray events. By including a g-ray absorber data that comprise a more balanced distribution of neutron and g events can be acquired and thus provide data sets with greater versatility. A 50 mm thick lead wall placed immediately in front of the detector attenuates the g component of the mixed field, and this resulted in a neutron/g ratio of (50:50). Fig. 2 is a plot of amplitude versus time for 20 superimposed pulses acquired using this experimental configuration. 2.1. Digital PSD analysis The PGA method of fast PSD for organic liquid scintillators using digital electronics was first presented by D’Mellow et al. [1]. This simple and inherently fast method for the classification between g-ray and neutron events in organic liquid scintillators exploits the principle that the light function of a scintillator due to a neutron interaction decays more slowly than that for a g-ray interaction. Therefore, the comparison of the gradient between the peak amplitude and the amplitude a specified time after the peak amplitude (referred to as the discrimination amplitude) on the decay face of the pulse can be used as a means of pulse discrimination. The sampling frequency of the pulse-amplitude data, collected using the setup described previously, was first reduced to a frequency more consistent with that achievable by an embedded processor. A modular system based on a field-programmable gate array (FPGA) is one possible choice of hardware platform that might be used in the future. For this work the sampling frequency was reduced to 250 MHz (i.e. every sixteenth sample was kept, equivalent to sampling every 4 ns). A digital filter and base-line correction algorithm were applied to each pulse.
ARTICLE IN PRESS M.D. Aspinall et al. / Nuclear Instruments and Methods in Physics Research A 578 (2007) 261–266
Clearly, for symmetric Gaussian data s ¼ 0. For so0 the distribution is skewed to the left and for s40 the distribution is skewed to the right. Kurtosis k indicates whether the data are peaked or flat in comparison with a Gaussian distribution. For Gaussian data k ¼ 3, for peak data k43 and for flat data ko3, where PN ðxi mÞ4 k ¼ i¼1 . (2) ðN 1Þs4 Often, k is normalized by subtracting 3 to constitute the excess kurtosis. Only the excess kurtosis is specified in this paper.
Fit
ms
ss
w2v
Kurtosis
Skew
g Neutron Double
32:9 0:3 85:4 0:5 –
10.87 17.27 –
1.265 1.095 1.157
32.51 4.84 –
7.76 0.57 –
8000 7000 6000 5000 4000 3000
In Fig. 3, the probability distribution histogram of the gradients between the peak amplitudes and the discrimination amplitudes at 10 ns after the peak is presented. This reduces further to a comparison of the discrimination amplitudes alone. The equation for the gradient takes the form m ¼ dy=dt ¼ ððy2 y1 Þ=ðt2 t1 ÞÞ, where y2 is the peak amplitude, y1 is the discrimination amplitude and dt 0.02 Distribution 0.015
Twin Gaussian γ Gaussian n Gaussian
0.01
0.005
0 0
50 100 150 Normalised, quantised sample amplitude
Fig. 3. A probability distribution histogram of the normalized and 8 bit quantized 10 ns discrimination amplitudes complete with Gaussian distributions calculated from the means (mg ¼ 32:9 and mn ¼ 85:4) and standard deviation (sg ¼ 10:9 and sn ¼ 17:3) estimated from these measurements. These data were produced by the digital PSD algorithm operating on 2436 pulses from a 50 : 50 n=g amplitude data set, with a 250 MHz sampling frequency, 8 bit amplitude resolution and a 20 point, running-average FIR filter.
Gamma ray events Neutron events
2000 1000 0
3. Results
Probability
Table 1 Measures of mean ms and standard deviation ss , reduced chi-square w2v , and measures of kurtosis and skew for the Gaussian fits to the experimental distributions shown in Fig. 3
Peak amplitude
The filter was a 20 point running-average finite impulse response (FIR) filter. Each of the pulses were then normalized by dividing every amplitude sample in the pulse by the peak amplitude sample. Finally pulses were quantized to an 8 bit amplitude resolution. In our analysis, skew s is defined as the measure of asymmetry of N data xi with mean m and standard deviation s PN ðxi mÞ3 s ¼ i¼1 . (1) ðN 1Þs3
263
0
1000
2000 3000 Sample amplitude
4000
Fig. 4. A scatter plot of the peak amplitudes against the 10 ns discrimination amplitudes before the amplitude data were normalized and quantized. The experimental arrangement was as described for Fig. 3.
is the time after the peak amplitude that the discrimination amplitude is measured. The time dt is specified and so remains constant for each pulse. Similarly the peak amplitude y2 is constant for each pulse because the pulseamplitude data were normalized and 8 bit quantized. The single varying parameter of the gradient m is the discrimination amplitude y1 . The left peak of the probability distribution in Fig. 3 corresponds to the g-ray events (faster decay face, hence the smaller discrimination amplitude) and the peak to the right corresponds to the neutron events (slower decay face, hence the larger discrimination amplitude). Gaussian fits have been applied to the experimental distribution and the mean ms and standard deviation ss are presented in Table 1. The uncertainties in the means ms have also been estimated. Values of ms and ss were used as estimates of m and s to calculate the Gaussian fits of Fig. 3. The curves were scaled to have the same area as the histogram, and thus they do not represent the true parent distribution but rather a sample distribution based on the experimental distribution. Table 1 includes measures of w2v , kurtosis and skew for the Gaussian fits shown in Fig. 3. In Fig. 4, a scatter plot of the peak amplitude versus discrimination amplitude prior to normalizing and quantizing the pulse-amplitude data is presented. Two characteristic clusters are apparent in these data. These two distinctive
ARTICLE IN PRESS M.D. Aspinall et al. / Nuclear Instruments and Methods in Physics Research A 578 (2007) 261–266
264
Table 3 The six-parameter values for Eq. (3) for the accurate reproduction of the average g ray and neutron pulse shapes measured by the EJ-301 liquid scintillator shown in Fig. 5
Normalised amplitude
1 Marrone's fit to average γ Marrone's fit to average n Average of 1254 gammas Average of 1182 neutrons
0.8
Radiation
A
B
y
ls
ll
t0
Gamma Neutron
10.53 11.32
0.0166 0.0166
4.368 4.317
3.482 3.537
11.52 38.14
0.3311 0.2864
0.6 0.4 0.2 0 0
10
20
30 40 Time,ns
50
60
70
Fig. 5. Average neutron and average g-ray pulse shapes measured by the EJ-301 liquid scintillator with a 1 GHz sampling frequency. The Marrone model has been fit to each pulse shape.
Table 2 Measures of sum of squares due to error (SSE), R-square and root mean squared error (RMSE) for function fit to an average neutron and an average g-ray pulse shapes measured by the EJ-301 liquid scintillator shown in Fig. 5 Fit
SSE
R-square
RMSE
Gamma Neutron
0.081347 0.028072
0.99687 0.99891
0.012629 0.007419
clusters were separated by eye and several points specified mid-way between them. A fit based on a double exponential curve was then applied to them. Having separated and categorized the mixed pulseamplitude data into neutron and g-ray events, we can now determine an average neutron and an average g-ray pulse shape measured by the EJ-301 liquid scintillator. These are shown in Fig. 5. L ¼ Aðeðtt0 Þ=y eðtt0 Þ=ls þ Beðtt0 Þ=ll Þ.
(3)
The six-parameter function that describes pulse shapes in scintillator detectors [2], given by Eq. (3), has been applied to the average neutron and average g-ray pulse shapes, and the corresponding fits are included in Fig. 5. Parameters for the corresponding g-ray and neutron functions are given in Table 3. These parameters were derived using the curve fitting tool provided by MATLABs . Table 2 displays the sum of squares due to error (SSE), R-square and root mean squared error (RMSE) values calculated for the fits made in Fig. 5. Fitting the function to the average pulse shapes was performed using the curve fitting tool available in MATLABs . The fitting algorithm used was Trust-Region and the starting values were set to the parameters derived and presented in Ref. [2]. For both the neutron and g-ray fit SSE relates closely to 0 indicating that there is a small total deviation of the response values from the fit to the response values. Similarly, the RMSE values are close to 0
which further supports the consistency between the fit and the response data. The R-square measure for both fits of 99.7% proves how successful the fit is in explaining the variation of the data. The average neutron and average g-ray pulse shapes are evidently smooth and quantitatively consistent with pulseshape description via the six-parameter model; Eq. (3) and the six-parameter values shown in Table 3. This demonstrates empirically that the noise on the individual pulses used in the averaging process is random as it cancels out. 4. Discussion 4.1. Digital PSD Although a fast digitizing oscilloscope was used to acquire the data in this research, the timing and amplitude resolution are consistent with a fast embedded processor at the high-performance extent of what is currently available. This indicates that in the future the use and significant cost of Flash ADC processors might be avoided as it becomes feasible to prototype single-board digitizers and/or exploit the now significant capability of fast acquisition scopes. With regard to the latter, the potential for built-in personal computing capacity will remove the need for a separate computer in the acquisition setup. The quality of discrimination is governed by several parameters, including the time after the peak at which the discrimination amplitude is measured, the sampling frequency, the amplitude resolution, the filter and the trigger level. It is important to minimize the time after the peak at which the discrimination amplitude is measured in order to limit data loss due to pile-up events. A detailed inspection of Fig. 5 suggests that the optimum sample to use for discrimination lies between 15 and 25 ns after the peak of the pulse, because this is the point of greatest separation between the two pulse shapes. However, it is likely that using a discrimination amplitude so late in the pulse would limit the ability to discriminate events before a possible pile-up event occurs. The sampling frequency, amplitude resolution and complexity of the filter should be kept to a minimum with the intention of short computation time on a stand-alone FPGA device. The trigger level should be set at a value to remove low amplitude, high signal-to-noise pulses and yet ensure a high detected proportion of the mixed field events. However, the reduction of these parameters reduces the
ARTICLE IN PRESS M.D. Aspinall et al. / Nuclear Instruments and Methods in Physics Research A 578 (2007) 261–266
quality of discrimination and it is, therefore, important that they are optimized with this in mind. Whilst the double Gaussian distribution in Fig. 3 clearly demonstrates that there are two separate Gaussian distributions, the left corresponding to the g-ray events and the right corresponding to the neutron events, there is a region of slight excess between the distributions on the g side. This is also corroborated analytically by the measured skew or the g distribution, which indicates an imbalance to the right as a result of the bulge to the right between the distributions. The first impression is that these events are long g events; however, whilst their origin is unknown their position in these data is consistent with their being piled-up g events. The high proportion of g in the mixed field used in this work would further support this hypothesis. Both the double and single Gaussian distributions correspond consistently with the data as evidenced by the w2v for each Gaussian distribution (Table 1). However, the excess kurtosis measured for the g-ray distribution indicates that it is significantly peaked, which is consistent with the relatively limited energy range of the g events associated with this experimental setup compared to that of the neutron component. This is perhaps more obvious in Fig. 4, where the g-ray plume is considerably shorter and exhibits a greater concentration around the lower amplitudes in comparison with the neutron plume. The measure of the kurtosis for the neutron distribution is peaked likewise, but to a lesser degree than that for the g-ray distribution. 4.2. Detector characterization Using the technique described in this work it is possible to obtain estimates for the functional characteristics of a liquid scintillator and PMT, with respect to the six parameters of Eq. (3), without resorting to the parametrization of thousands of pulses and the subsequent averaging of thousands of six-parameter data sets. The extraction of goodness-of-fit parameters provides an indication of the correspondence of the resulting characteristic pulse shapes and, if required, the six parameters can then be determined from the single, mean pulse. This characterization of a scintillator detector is achieved without any prior knowledge of the detector or PMT and is based entirely on events from a source that has sufficient contributions of both neutrons and g rays. Characteristic pulses for both neutron and g components of a mixed field can be obtained in under an hour. 5. Conclusions The characterization of detection responses from liquid scintillators often requires a knowledge of the size of the scintillant cell, the response of the scintillant compound and the influence of the PMT. This can be especially relevant where the scintillant is being applied to neutron detection, and the discrimination of neutrons and g rays is required on the basis of pulse shape.
265
The six-parameter formula [2] can be used to extract the pulse shape if such a fit can be applied to a number of pulses that is large enough to be representative of the field under study and the average of the fit parameters extracted. This is feasible where one or perhaps a few detectors are being used, for example in large-scale nuclear physics experimentation. However, in industrial applications where a diverse range of liquid scintillation detectors often need to be characterized and with the recent resurgence of interest in scintillation arrays for applications in, for example, homeland security, such an approach is likely to be inconvenient and difficult to automate. Furthermore, the extent of the six-parameter description will often exceed the requirements of industry-based system designers who require the generic pulse-shape response for a given detector design in order to establish neutron/g discrimination thresholds. This requirement is becoming more important with the advent of fast digitizers, where the entire pulse shape can be sampled and recorded, and in circumstances where the pulse shape is often corrupted by noise. In the analogue domain it is rarely easy to extract pulseshape information with the level of detail required for detector characterization. Often, only single-valued analogue variables arising as a result of, for example, hard-wired calculus functions can be extracted such as the area and the amplitude of a given pulse. Thus, the extraction of exemplar pulse shapes on which to base discrimination limits and pulse pile-up thresholds is rarely possible and, instead, these levels are often derived in conjunction with the real-time output of the detector. This method requires skill, judgement and is difficult to replicate from measurement to measurement. Fundamentally, this approach is not compatible with the embedded control of large numbers of scintillation counters and the availability and know-how associated with analogue discrimination modules has been in decline for some time. The method described in this work enables the pulseshape response of liquid scintillation counters to be characterized quickly and easily from the average of many digitized pulses. The resulting normalized pulses for g-ray and neutron event constitute generic pulses shapes for a given detector. These can then be used as comparison events with which to determine the nature of the incoming data, either in real time or via post-processing. The shape could be ported to a given detector of the same type and embedded on system-based electronics. It can also be used to extract the Marrone coefficients if these are required for a deeper analytical description of the detector response. This method is a novel example of the flexibility and power that digital processing can bring to the detection and measurement of mixed radiation fields. Acknowledgements The authors would like to acknowledge the help of Mr Paul Rossiter for providing access to and assistance
ARTICLE IN PRESS 266
M.D. Aspinall et al. / Nuclear Instruments and Methods in Physics Research A 578 (2007) 261–266
with the neutron source in the Department of Physics at the University of Liverpool. This research was funded by the Engineering and Physical Sciences Research Council (EPSRC) via the DISTINGUISH research consortium (http://www.distinguish.org.uk/) under the EPSRC Think Crime! initiative. References [1] B. D’Mellow, M.D. Aspinall, R.O. Mackin, M.J. Joyce, A.J. Peyton, Nucl. Instr. and Meth. A, in press, 2007, doi:10.1016/j.nima. 2007.04.174.
[2] S. Marrone, et al., Nucl. Instr. and Meth. A 490 (2002) 299. [3] D. Wolski, M. Moszynski, T. Ludziejewski, A. Johnson, Nucl. Instr. and Meth. A 360 (1995) 584. [4] S.C. Wang, C.C. Hsu, R.W.S. Leung, S.L. Wang, C.Y. Chang, C.P. Chen, K.C. Cheng, T.I. Ho, W.P. Lai, H.M. Liu, Nucl. Instr. and Meth. A 432 (1999) 111. [5] P.L. Reeder, A.J. Peurrung, R.R. Hansen, D.C. Stromswold, W.K. Hensley, C.W. Hubbard, Nucl. Instr. and Meth. A 422 (1999) 84. [6] H. Klein, S. Neumann, Nucl. Instr. and Meth. A 476 (2002) 132. [7] J.S. Neal, W.L. Bryan, C.L. Britton, J.D. Edwards, S.A. Pozzi, J.T. Mihalczo, Institute of Nuclear Materials Management 44th Annual Meeting, Phoenix, Arizona, 2003. [8] S.D. Jastaniah, P.J. Sellin, Nucl. Instr. and Meth. A 517 (2004) 202.