A Functional Form for Liquid Scintillator Pulse Shapes

F.Q.L. Friesen∗, C.R. Howell

Department of Physics, Duke University and the Triangle Universities Nuclear Laboratory
(TUNL). Durham, North Carolina 27708, USA
arXiv:1912.07682v1 [[Link]-det] 16 Dec 2019

Abstract

Digitization of detector signals enables analysis of the original waveform to

extract timing, particle identification, and energy deposition information. Here
we present the use of analytical functions based on sigmoids to model and fit
such pulse shapes from liquid organic scintillators, though the method should
also be applicable to other detector systems. Neutron and gamma interactions
in NE213 detectors were digitized from the phototube anode and fit using a
sigmoid-based function. The acuity of the fit in extracting timing information
and performing neutron-gamma pulse-shape discrimination are presented and
discussed.
Keywords: neutron, scintillator, pulse-shape discrimination, sigmoid, NE213,
BC501A

1. Introduction

Organic liquid scintillators are the dominant technology used in experiments

requiring fast neutron detection. A wide variety of scintillators have been de-
veloped offering different light output characteristics. A common feature is a
rapid rise in amplitude, followed by a slow exponential return to baseline. The
decay tail may consist of multiple exponential components depending on the
formulation of the liquid and particle being detected. For example, NE213 and
the equivalent formulations of EJ301 and BC501A, are designed for fast neutron

∗ Corresponding author: fqf@[Link]

Preprint submitted to Nuclear Instruments and Methods in Physics Research Section ADecember 18, 2019
counting with good pulse-shape discrimination [1] (PSD) capability for distin-
guishing gamma from neutron interactions in the detector [2]. The sensitivity
of the decay time of the emitted light to the type of particle interaction in the
liquid is due to the dependence on the mechanism of excitation. Gamma rays
generally scatter from electrons, which then preferentially cause the liquid to
fluoresce. Neutrons generally scatter from protons in the organic molecules,
and the resulting excitations produce phosphorescence in addition to fluores-
cence [3]. The lifetime associated with phosphorescence is typically orders of
magnitude larger, and pairs of the corresponding triplet states recombine with
each other leading to a final delayed de-excitation by fluorescence [4].
It is desirable to have a functional form to generically represent neutron and
gamma pulse shapes. In addition to describing the decay, the function should
accurately represent the rising edge of the pulse. It should have few free param-
eters, each with a clear physical interpretation. Uncertainties in the quantities
of interest are then directly evaluated by the fit. With these features, it is pos-
sible to perform a fit on a given waveform, alter a parameter in a physically
interpretable way, and regenerate a different but realistic signal. In this way,
neutron and gamma waveforms with arbitrary timing can be quickly generated
for simulating events in the same data format produced by the experimental
apparatus. This ultimately enables processing of simulation output and mea-
surements using identical code.
A frequently used expression for the pulse shape has the form L = A(e−θ(t−t0 ) −
e−τf (t−t0 ) ) + B(e−θ(t−t0 ) − e−τs (t−t0 ) ) [5]. While the number of free parame-
ters is reasonable, obtaining reliable convergence during the fit process can be
challenging. This difficulty increases if the location of the peak in time isn’t
well constrained. The initial (concave up) rise from baseline is also not mod-
eled by this expression. In sufficiently fast scintillators, this initial rise can be
non-negligible. In this work we examine the use of sigmoids for describing de-
tector pulse shapes. They capture the exponential behavior with the benefit of
remaining bounded. The proposed function is the product of several sigmoids
representing the key features of the detector signals, e.g., the rise and decay

2
Figure 1: Example raw waveforms from neutron and photon interactions with similar am-
plitudes. Labels for the time constants are positioned to indicate the regions dominated by
each variable. The key difference between the interactions is the prevalence of the slow decay
component τs as parameterized by R.

times.
The widespread use of digitizers, supported by cheap commodity storage,
have made it common practice to keep detailed records of every critical signal in
an experiment with nanosecond precision. With this capability, PSD analysis no
longer requires dedicated hardware, e.g., using fixed width integration windows
[2]. More detailed signal analyses using fits becomes possible.
The traditional, time-tested PSD analysis method uses the charge compari-
son technique by computing the ratio of charge collected from the detector across
two different fixed windows which correspond to the total charge in the pulse
and the charge in the fast component (i.e., the total charge minus the charge in
the tail)[1] [6]. A drawback of this method is that it can be difficult to find a
pair of windows which work robustly over a wide range of pulse heights. Even

3
though decay times are independent of amplitude, the pulse height of a signal
determines how quickly the tail will become comparable to noise. This means
that large signals are best characterized with longer windows; the opposite is
true for small signals. This can lead to undesirable trade-offs in experiments
when the detector threshold is low and neutrons and gammas need to be reliably
distinguished across a large dynamic range.
This work was performed in support of 3 He(γ, pn) measurements at TUNL.
The functions are primarily used to generate digitizer-like data from the output
of GEANT4 simulations. The functions can also be used in signal analysis by
using them to fit the raw digitizer data. However, the use of these functions in
data analysis would come at the expense of computation time relative to other
methods.

2. Experimental Setup

Neutron detector efficiency calibration data were obtained using the 3 H(p, n)
reaction at a polar angle of 136.4◦ relative to the incident proton beam axis. The
digitized PMT signals were fitted using the proposed function. A proton beam
was incident on a tritiated titanium foil target, producing neutrons at 0.8 MeV,
1.0 MeV, 1.2 MeV, 1.4 MeV, and 1.6 MeV. The 1.6 MeV data are presented
here, as they provide the largest range of pulse heights. The backward angle
was chosen by requiring a higher proton beam energy than used for forward
angle neutron production to reduce the energy spread of the neutrons.
The detectors used in this study were filled with BC501A and coupled to
Hamamatsu R1250 photomultiplier tubes. The interior of the liquid scintillator
cells were cylindrical with diameter 12.7 cm and depth of 5.1 cm. The PMT
anodes were connected directly to sis3316 digitizers. The digitizers have 12
bit resolution with a 4 ns sampling period. Waveforms were recorded for 240
ns such that the first 76 ns consisted of baseline as shown in Fig.1. The full
duration of sampled data was used for fitting.

4
Figure 2: The red solid line shows a representative neutron pulse shape as constructed from
constituent sigmoids with R = 0.914. The slow decay contribution is non-negligible.

3. Pulse-Shape Analysis

For BC501A, the relative weighting of the importance of the two decay rates
is a free parameter in the fit. This ratio R is the fraction of decay light output
produced by the fast mechanism. The expectation is that long decay times
should be much less prominent in gamma events than neutron events. There
are 6 independent parameters to fit for a PMT signal from a scintillation event
in BC501A. The first 3 are fixed for a given PMT and scintillator combination,
the remaining 3 are fit to each individual waveform (see Fig.1):

• τr = Rise time (∼1 ns)

• τf = Fast decay time (∼5 ns) representing fluorescence

• τs = Slow decay time (∼50 ns) representing return from phosphorescence

• R = Relative weighting of decay times 0 ≤ R ≤ 1, e.g. R=1 means that

the amplitude of the slow decay component is zero

5
 • A = Pulse amplitude

• t0 = Pulse arrival time

y(t) = Af (t)[Rg(t) + (1 − R)h(t)] (1)

where

1
f (t) = t−t0 (2)
−
e τr +1

1 1
g(t) = 1 − t−t = t−t0 (3)
− τ 0 τf
e f +1 e +1

1 1
h(t) = 1 − t−t0 = t−t0 (4)
−
e τs +1 e τs +1
The constituent functions of Eq.1 are standard logistic functions [7]. Those
representing the decay of the pulse have been subtracted from unity. Eq.2
describes the early exponential rise from baseline at the start of the pulse and
part of the round-off as the amplitude approaches its maximum. Note that there
is no time offset between functions f , g, and h. This means that when the rise
time is much smaller than the decay times, then function f has negligible effect
on the tail.
Example waveforms are shown in Fig.1 with labels indicating the features
determined by each parameter. Waveforms were prepared for the fit by first
ensuring that the baseline level averages to 0. The subtraction was performed
on an event-by-event basis. The curve fit function from the [Link] [8]
python library was used to perform the fit with bounded variables. Examples
of the functions described by Eqs.1-4 are shown in Fig.2 for experimentally
determined parameter values. The following bounds (in ns) were used to help
the fit converge quickly:

0≤R≤1 A≥0 0 ≤ τr ≤ 6 0 ≤ τf ≤ 24 40 ≤ τs ≤ 100 (5)

6
 Pulse Height (MeV )
ee
3
10
1.8
1.6
1.4
1.2 102

1
0.8
0.6 10

0.4
0.2
0 1
0.86 0.88 0.9 0.92 0.94 0.96 0.98 1
R

Figure 3: The R parameter plotted vs. waveform pulse height with τf and τs held constant
during the fits to template generated samples. The lower R values are more neutron-like. The
lower pulse height threshold corresponds to 60 keVee (1/8 137 Cs edge). This histogram was
accumulated with an incident neutron energy of 1.6 MeV.

The waveforms taken from data were divided into three categories: tem-
plates, general samples, and template generated samples. The templates were
137
those with a pulse height equivalent to the Cs edge or larger. Traditional
PSD works very well in this region, so gamma pulses and neutron pulses can
be unambiguously tagged by type before being fit. The general samples were
drawn from the same pool of waveforms as the templates, but with a much
137
lower pulse height cut (1/8 of the Cs edge). General sample waveforms were
used to evaluate the timing capability of the fit. Template generated samples
were constructed from templates by scaling down a template waveform ampli-
tude to the obtain the desired pulse height, and then adding Gaussian noise to
reproduce the original detector noise level (the pulse height spectrum was set to
duplicate that of the general samples). The template generated samples were

7
Figure 4: Example of normalized mean gamma and neutron pulse shapes obtained from fits
to template waveforms.

used to evaluate PSD capability.

The templates were fit individually using the bounds shown in Eq.5. The
mean τf and τs values from neutron template fits were used to fix these pa-
rameters for all subsequent fits (all γ waveforms, general samples, and template
generated samples). The parameters R, A, τr and t0 were left free. The R can
be interpreted as a PSD parameter (Fig.3). The PSD capability of R is im-
proved if τr is kept as a free parameter in the fit. Representative pulse shapes
for neutron and photon events were generated as shown in Fig.4.
Finally, the template scaling method provides a means to test PSD perfor-
mance. Templates with known identities are scaled down and added to back-
ground noise to simulate lower pulse height events. Pools of realistic waveforms
with low pulse-heights can be produced from templates with known identities.
Using these generated events, the pulse identification capability can be quan-
titatively evaluated. PSD results for two sets of waveforms at different pulse
height thresholds are shown in Fig.5. The R value from the fit is compared
to a traditional PSD calculation in Fig.6. The traditional PSD parameter was
calculated as the quotient of two baseline-subtracted waveform integrals. They

8
Figure 5: Histogram of fit (PSD) parameter R from a pool of template generated waveforms.
The upper lines are for a neutron pulse height cut of 1/8 Cs, the lower lines are at 1/4 Cs. The
left distributions show known neutron inputs. The right distributions show known gamma
inputs. The gaussian fits used to calculate the FOM at 1/8 137 Cs are also shown.

extend to 20 and 4 samples past the leading edge, respectively (i.e., 172 ns and
108 ns in Fig.1). There is good agreement over the classification of the vast
majority of waveforms.
The result of fixing the discrimination threshold on R at 0.98 is shown in
table 1. The table was generated from a pool of 1.2 x105 tagged waveforms
(generated by template scaling) with a 1/8 Cs (60 keVee ) pulse height threshold
limit, where half are neutrons and the other half are photons. At this threshold,
2 γ-ray signals leak through the PSD cut for every 1000 photons detected. For
every 1000 incident neutrons, 3 are misidentified as photons. The reduction
in total number of neutron and gamma events as the threshold is increased
corresponds to the reduced detection efficiency. The values and uncertainties
of parameters obtained from fitting the template signals are given in Table 2.
A figure of merit (FOM) value was also calculated for the R discrimination

9
Table 1: Neutron/γ discrimination capability as a function of pulse height threshold. The
discriminator threshold for R was set to 0.98 for this analysis. Increasing it reduces rejection
of true neutron events at the expense of mis-identifying more γ events as neutrons. The
uncertainties are statistical.
Threshold Correctly Correctly # n above # γ above
(MeVee) Identified n Identified γ Threshold Threshold
0.0598 0.997±0.006 0.998±0.006 6.0x104 6.0x104
0.12 0.997±0.008 0.998±0.009 3.4x104 4.3x104
0.239 1.00±0.01 1.00±0.03 1.1x104 2.5x104
0.478 1.00±0.03 1.0±0.1 2.1x103 1.2x104

parameter by fitting the neutron and gamma peaks with gaussian functions and
using the definition FOM= (µγ − µn )/(σγ + σn ) where σ refers to the FWHM
of each gaussian [9]. The traditional charge comparison PSD method had FOM
values of 1.34 and 1.43 for pulse height thesholds of 1/8 and 1/4 Cs, respectively.
The FOM values for the R fit method were 1.39 and 1.52 calculated on the same
data.

Table 2: Mean values and standard deviation of the parameters obtained by fitting templates.
The pool consisted of 12,000 waveforms with an equal number of neutron and photon events.
The uncertainty in t0 is given as an estimate of the time resolution of this method because the
absolute value of t0 depends on the source of the particle. The τf and τs values for photons
were constrained to be the same as the mean values obtained from neutron template fits.
Therefore, these fixed parameters do not have uncertainties for photons.
Parameter Neutrons γ-rays
τr (ns) 1.46 ± 0.29 1.37 ± 0.26
τf (ns) 5.45 ± 0.63 5.45
τs (ns) 52 ± 11 52
R 0.949 ± 0.035 0.992 ± 0.050
t0 (ns) ± 0.17 ± 0.11

10
 1
R
0.98
3
10

0.96

0.94
102

0.92

0.9
10

0.88

0.86
1
0.5 0.6 0.7 0.8 0.9 1
Traditional PSD

Figure 6: The PSD parameter calculated by the window integration method plotted against
the R fit parameter for the same data. The peak in the upper right hand corner corresponds
to γ events.

4. Timing Analysis

The t0 fit parameter was compared with a software constant fraction discrim-
inator [10] (CFD) value computed with a trapezoidal filter. The data sampling
period was 4 ns. The peak and gap times were both 2 samples long (8 ns), and
the time was interpolated on the falling edge of the filter at 1/2 of the maximum
output. The relationship between the two, shown in Fig.7, indicates that the
raw t0 value from the fit closely follows the CFD value independent of event
type or amplitude. Slices taken either perpendicular or parallel to the CFD
axis in Fig.7 corresponding to 85.6 ≤ t0 ≤ 89.6 ns (100 to 104 ns on the CFD
axis) were used to compare timing resolution. These are experimental data,
so different arrival times have different physical interpretations, and different
slices should not be directly combined. The standard deviation of the points
in the vertical slices from the fit is 0.33 ns within the range, compared to a
standard deviation of 0.35 ns as determined with the CFD method (horizontal

11
slices passing through the same points).

92
t 0 from Fit

500

90
400
88

300
86

200
84

82 100

80 0
96 98 100 102 104 106
Trapezoidal Filter CFD

Figure 7: Trapezoidal filter CFD vs t0 from the sigmoid-based fit applied to general samples.
The pulse height threshold corresponds to 60 keVee (1/8 137 Cs edge). Slices passing through
the bin with the most counts (102.6 ns, 88.1 ns) were used to compare timing resolution.

5. Conclusions

A sigmoid-based functional form was developed for representing pulse shapes

from liquid scintillation detectors. The rise time, amplitude, time offset, and
decay components are represented in a physically interpretable way. Contribu-
tions from different decay times can be reliably identified by the fit, thereby
enabling PSD analysis of each waveform. The bounded nature of sigmoids sim-
plifies construction of plausible functions, and the function presented here can
be further generalized to allow for more decay components as needed. The
function parametrized in Eqs.1-4 can be used to generate signals to realistically
represent pulse shapes from liquid scintillator detectors. Waveforms associated
with BC501A can be generated using values given in table 2. Other sigmoids
could also be explored for constructing the fit function.

12
References

[1] J. Adams, G. White, A versatile pulse shape discriminator for charged

particle separation and its application to fast neutron time-of-flight spec-
troscopy, Nuclear Instruments and Methods 156 (3) (1978) 459 – 476.

[2] A. Comrie, A. Buffler, F. Smit, H. Wrtche, Digital neutron/gamma discrim-

ination with an organic scintillator at energies between 1MeV and 100MeV,
Nuclear Instruments and Methods in Physics Research Section A: Acceler-
ators, Spectrometers, Detectors and Associated Equipment 772 (2015) 43
– 49.

[3] J. Birks, The theory and practice of scintillation counting, International

series of monographs on electronics and instrumentation 27.

[4] P. Schuster, E. Brubaker, Investigating the anisotropic scintillation re-

sponse in anthracene through neutron, gamma-ray, and muon measure-
ments, IEEE Transactions on Nuclear Science 63 (3) (2016) 1942–1954.

[5] S. Marrone, D. Cano-Ott, N. Colonna, C. Domingo, F. Gramegna, E. Gon-

zalez, F. Gunsing, M. Heil, F. Kppeler, P. Mastinu, P. Milazzo, T. Pa-
paevangelou, P. Pavlopoulos, R. Plag, R. Reifarth, G. Tagliente, J. Tain,
K. Wisshak, Pulse shape analysis of liquid scintillators for neutron studies,
Nuclear Instruments and Methods in Physics Research Section A: Acceler-
ators, Spectrometers, Detectors and Associated Equipment 490 (1) (2002)
299 – 307.

[6] J. Heltsley, L. Brandon, A. Galonsky, L. Heilbronn, B. Remington,

S. Langer, A. V. Molen, J. Yurkon, J. Kasagi, Particle identification via
pulse-shape discrimination with a charge-integrating adc, Nuclear Instru-
ments and Methods in Physics Research Section A: Accelerators, Spec-
trometers, Detectors and Associated Equipment 263 (2) (1988) 441 – 445.

[7] D. H. von Seggern, CRC Standard Curves and Surfaces with Mathematica,
Chapman and Hall/CRC, 2016.

13
 [8] E. Jones, et al., SciPy: Open source scientific tools for Python (2001).
URL [Link]

[9] S. Bose, M. Chatterjee, B. Sinha, R. Bhattacharyya, Neutron-gamma dis-

crimination based on leading edge shape measurement, Nuclear Instru-
ments and Methods in Physics Research Section A: Accelerators, Spec-
trometers, Detectors and Associated Equipment 270 (2) (1988) 487 – 491.

[10] M. Nakhostin, Signal Processing for Radiation Detectors, John Wiley and
Sons, 2017.

14