Smooth filter

Because the detector signals are always distorted by random noise, the goal of this procedure is to reduce the high-frequency noise while not significantly changing the detector signals to preserve the signal characteristics and improve the data SNR. This procedure was performed initially because it may slightly alter the raw waveform and affect the performance of the other procedures.

The moving average filter (MAF) is the most commonly used smooth filter in the time domain owing to its simplicity, ease of implementation, and rapid execution speed. In addition to the aforementioned advantages, it offers the lowest noise for a given edge sharpness for any linear filter [18]. Other optional filters include the Savitzky–Golay [29], binomial [28, 30], Whittaker [31], and the Kalman filters [32, 33]. Although these smooth filters perform well in multiple fields such as Raman and Mössbauer spectroscopy, complex algorithms make them difficult to implement in the FPGA. For this reason, the MAF was chosen as a smooth filter for the CsI(Tl) array.

The only parameter for MAF is the number of samples averaged over. Its value should be chosen with care because as the value increases, the noise decreases, while the edges become less sharp. The aforementioned test bench, which was operated using cosmic rays, was used to acquire raw waveform data, and the MAF algorithms were executed with different parameter values in the off-line analysis. The results are presented in Fig. 2. The reason for using a power of two for the parameter values is that the division is simpler to implement in the FPGA using only shift operations. Finally, owing to its good performance in high-frequency noise reduction and small waveform changes, the value of the parameter was set as equal to 8.

Fig. 2Full-screen View

(Color online) Comparison among MAF results with different values of N. Inset illustrates leading edges in the picture.

Baseline restoration

In many cases, the baseline is assumed constant at all times. This enabled the baseline value to be measured at any time while the DAQ system was idle during the experiment. However, this assumption is only true when the signal length is short and the count rate is relatively low. This is certainly not the case in this instance, because the CSA signal extracted has a long tail. A good solution is to use optimal filters [34]. However, complicated arithmetic discourages us from proceeding further. To simplify the procedure, it was assumed that the baseline remained constant within each selected dataset, which contained only one signal waveform with several baseline data points in front. This baseline data can subsequently be used to calculate the baseline value of each individual event. Reference [35] showed that two values, the average and median of the baseline data, can be used to evaluate the baseline level. The method for obtaining average baseline data is clear. The median is the exact quantity in the middle of the baseline dataset when ordered. Reference [35] showed that the median is a better estimate of the baseline level than the average over a wide range of count rate loads. The other two methods, which are called "averaging over the selection set" and "averaging over the flat chunk selection set" respectively, are also introduced in [35]. These two methods are identical in terms of how to proceed, albeit differ in terms of data selection. Further detail regarding both methods can be found in [35], and the main procedures are summarized in the flowchart shown in Fig. 3. Over a wide range of count rate loads, these two methods provide better estimates of the baseline than either the average or the median alone. In this case, because the baseline dataset has already been selected, the procedure list in Fig. 3 can be processed directly. Therefore, both methods are treated as a single method and are hereafter called iterative methods.

Fig. 3Full-screen View

The main procedures of the methods, averaging over the selection set and averaging over the flat chunk selection set, from [35].

Here, there are three methods: the average of the baseline dataset, median of the baseline dataset, and the iterative method. To evaluate performance, the baseline samples in the data used in the previous section were processed using these methods. The results of the baseline recoveries are shown in Fig. 4. The discrete and continuous histograms are displayed using the median and average methods, respectively, because the data types for these two methods are "int" and "float.” The histogram obtained using the iterative method is similar to that obtained using the median method. The only difference is the width of each discrete part, which is equal to twice the minimum difference parameter set in the algorithm(δz in Fig. 3, here the value is 0.1, further details can be found in [35]). Gaussian functions were used to fit the envelopes of the three histograms; the results are listed in Table 1. Because the data are all integers for the median method, the most probable value(MPV) of the histogram is treated as zero and not the mean fitting parameter. From a comparison of the MPVs of the three methods, it can be concluded that all three methods are good estimators of the baseline levels, and therefore, fully meet the performance requirements.

Fig. 4Full-screen View

Typical results of baseline restorations obtained with the following algorithms: (a) Median, (b) Average, and (c) Iteration.

Algorithms subsequently become the focus of the selection criteria. The iterative method was the most complex of the three methods. However, this method is powerful because it can continue to provide reasonable baseline values when the signal samples are also included the dataset [35]. The median method is simpler, albeit the data ranking algorithm is not "FPGA friendly”. Therefore, for the CsI(Tl) array DAQ system, the average method is preferred.

Arrival time

According to the ETF design scheme, the CsI(Tl) array was intended to generate a trigger for the entire system. With the complete digitization of the input waveform and the powerful calculation capability of the FPGA, the CsI(Tl) array trigger signal can be generated and controlled by software. This simplifies the electronic system by eliminating the necessity for additional electronics, such as splitters and time discriminators.

An ETF trigger system consists of 2-levels trigger generators [36]. At the front end, primary trigger signals were generated according to the signal shapes and logical relationships between the readout elements of a particular detector included in the trigger system. These primary trigger signals were subsequently fed into the global trigger logical unit to generate an event trigger signal for the entire system according to the physical interests. To improve logic operation effectiveness for signals with different delay times and time jitter, each primary trigger signal can be delayed and widened separately in a global trigger logical unit. From this perspective, time resolution is not the key element for the CsI(Tl) array.

Several methods are available to determine the arrival time of a signal. Common methods, which are identical to those in the analog scheme, are leading-edge discrimination (LED) and constant fraction discrimination(CFD). LED is the simplest method; however, it has a large time jitter owing to the time-walk effect. To achieve better performance, correction should be performed, which is a significant task for the FPGA programmer. Therefore, CFD was introduced to eliminate the time walk effect. There are two methods for implementing this algorithm in an FPGA: constant-fraction zero-crossing (CFDzc) and digital contact-fraction discrimination (dCFD) [37]. CFDzc is a digital version of the classic analog CFD [38], whereas dCFD is similar to an LED with a different threshold value, equal to the constant fraction of the signal amplitude [39]. Better performance in terms of time resolution was obtained with the CFDzc method than with the dCFD method [40]. Other methods exist for determining the arrival time of a signal, such as the RC-CR2 filter [41] and pulse-shape fitting. However, utilization of these algorithms is resource-intensive in the FPGA and considered outside the scope of this research.

Because time resolution is not the key criterion for the CsI(Tl) array, the LED method without interpolation was chosen to determine the arrival time, considering the expected high computational resource consumption of the energy extraction and PID procedures. Figure 5 shows the relationship between the amplitude of the waveform and arrival time using the previous data. The arrival time is the difference between the reference time and the time when the first sample point is above the threshold. The time reference was selected as the point at which the ADC waveform passed through 90% of its full amplitude at the leading edge. Linear interpolation was used to reduce the reference time jitter. A clear dependency between the waveform amplitude and the arrival time is shown in Fig. 5. Another phenomenon is that the lower the waveform amplitude, the greater the time jitter. This is because a low amplitude causes the leading edge of the waveform to become flat, and it is difficult to determine the time at which the waveform crosses the threshold and the time reference. Although the time jitter is as large as approximately 3 μs when the waveform amplitudes are extremely small, this may not be useful because a threshold can be set in the FPGA to determine the recorded signals. The time jitter of the waveforms with relatively large amplitudes was approximately 342 ns (the σ value with a Gaussian fit).

Fig. 5Full-screen View

(Color online) The relationship between waveform amplitude and time obtained with LED method. Red line shows the fit result with a cubic polynomial, used to correct the time jitter.

Further analysis can be performed with time-walk correction, which was not implemented in the FPGA in the proposed scheme. A cubic polynomial, shown in Fig. 5 with a red solid line, was used to fit the 2-dimension histogram. The arrival time histograms with and without time-walk correction are shown in Fig. 6. Good correction can be obtained with small waveforms, and the time jitter is approximately 344 ns(the σ value with a Gaussian fit), which is almost identical the value obtained with large waveforms.

Fig. 6Full-screen View

Time spectra obtained using LED technique. Blue dash line is the spectrum without time correction, red solid line is the spectrum with time correction.

Energy

In general, two quantities are used to extract the energy loss from the detector: the amplitude and total charge of the signal. These two quantities were measured using ADCs and charge-to-digital converters (QDCs) in a conventional DAQ system. Using a digital approach, these measurements are replaced by appropriate algorithms in the FPGA. However, the signal extracted from the CSA was too wide, indicating that the digital QDC method was not suitable for the CsI(Tl) array.

The most direct approach to extract the signal amplitude is to determine the maximum (positive signal) or minimum (negative signal) waveform. However, the measured amplitudes were significantly influenced by noise. To improve the energy resolution, digital filters, which perform the same function as shaping amplifiers in an analog measurement system, were used to shape the signal before amplitude extraction. Theoretically, the best signal shape that maximizes the SNR is the infinite-width cusp [42]. However, a practical filter of this type is the finite-width cusp filter that limits the amplitude measurement of a single signal to a specific time. The algorithm uses a different function instead of an infinite exponential function and performs truncation [43, 44], which means that the performance is reduced. Other commonly used shaping filters include a series of trapezoidal [45-51] and CR-RC^m filters [52-54]. To identify a suitable shaping filter, the γ-rays produced by a ⁶⁰Co source were measured on a test bench operating in waveform mode, and the recorded data were processed with various digital filters. The waveforms before and after applying the shaping filters are shown in Fig. 7. All parameters in each filter were optimized using repeated trials. The energy resolutions of the full-energy peaks for 1.33 MeV γrays were calculated as the criterion, and the results are shown in Table 2. The energy resolutions achieved by each filter were comparable with the performance of the finite cusp filter being slightly better. This is because the noise in the output test was extremely low, and the SNR values did not improve significantly. It should be noted that the noise level in the test is of the same order of magnitude as that of the ETF, indicating that all filters meet the requirements in terms of performance.

Fig. 7Full-screen View

(Color online) Shapes of raw waveform (a) and those modified by the following filters: (b) cusp, (c) trapezoidal with a single exponential function, (d) trapezoidal with a bi-exponential function, (e) MWD, (f) CRpz-RC2, and (g) CRpz-RC4.

In [55], it is concluded that the family of trapezoidal filters offers good performance and simpler implementation among many other shaping filters, which influences this choice. Many algorithms exist for the implementation of trapezoidal filters, of which moving-window deconvolution (MWD) [49-51] is the simplest for FPGA implementation. Therefore, the MWD filter was selected as the shaping filter for the CsI(Tl) array, and the corresponding ⁶⁰Co energy spectrum is indicated by the red dashed line in Fig. 12.

Fig. 12Full-screen View

Energy spectra of ⁶⁰Co radioactive source. The histogram shown by the blue solid line is obtained by the DAQ system algorithm mode and the on-line MWD algorithm, while the histogram shown by the red dashed line is obtained by the raw waveform mode and the off-line MWD algorithm.

PID algorithms

As previously mentioned, performing a pulse shape analysis for the CsI(Tl) array can reduce the background of charged particles in the energy spectrum of γ-rays and improve the PID performance through a combination with Δ E-E method. The basis of pulse shape analysis is that the ratio of the fast and slow components of light generated by the CsI(Tl) crystal depends on the type of incident particles, resulting in different shapes of the output waveform. Multiple methods exist for extracting these differences, among which the digital charge comparison [56-59] and rise-time comparison methods [59, 60] are widely used and straightforward to implement in the FPGA. Another method worth focusing on is called reconstructive particle identification (RPID) [51]. One reason is that this method was developed for CALIFA, whose construction is similar to the CsI(Tl) array. The RPID method can be successfully migrated into the DAQ system with significant potential. Another reason is that RPID can directly extract the fast and slow components of the CsI(Tl) crystal, thereby achieving a more improved performance.

To compare these PID algorithms, a triple α source (²⁴⁴Cm, ²³⁹Pu, and ²⁴¹Am) and a ⁶⁰Co source were used separately to irradiate the elemental detector of the CsI(Tl) array under identical conditions. The CsI(Tl) crystal envelope was pierced with a small hole to allow the α-particles to access the crystal. With respect to the digital charge comparison method, because the input signals are extracted from the CSA, which results in information on the incident particles being contained in the leading edges of the digitized waveforms, both time windows for charge integration were set to include the leading edge of the waveform with the same starting point, as shown in Fig. 8. The starting point is determined using the dCFD method without interpolation. The integration windows size was determined through repeated trials. With respect to the rise-time comparison method, the rise time is defined as the time frame between the points where the waveforms cross 12.5 and 87.5% of the full amplitude in the leading edges. Because the method performance is correlated with the time accuracy, the linear interpolation technique was used to determine the cross-points. Regarding the RPID method, some parameters are identical to those of the MWD filter previously described, with three additional parameters: the fast and slow decay times of the CsI(Tl) crystal and the time window for the second MWD procedure. The fast and slow decay times are in accordance with [51], because they are almost constants. The length of the second time window was determined using repeated trials.

Fig. 8Full-screen View

Two integration windows in digital charge comparison method.

The two-dimensional PID spectra for each method are shown in Fig. 9. Superficially, α and γ-rays were well identified for all three methods. The spectrum shown in Fig. 9(a) illustrates that the rise time and the amplitude of the waveform are weakly correlated and can be considered independent. The red dashed lines in Figs. 9(b) and (c) indicate that the spectra approach zero when both the horizontal and vertical axes are reduced. These relationships indicate that the PID parameters listed in Table 3 can be used to convert two-dimensional spectra into one-dimensional histograms, and the figure-of-merit (FoM), as defined in [61, 62], can be used to quantify the separation performance of the PID methods. These one-dimensional histograms are shown in Fig. 10 and the relevant FoMs are illustrated. Although data in the low-energy regions, which can make the FoM values larger than those in the beam experiments, were not included, the FoM values shown in Fig. 10 remain good references for evaluating these three methods. Of the three methods, the rise time comparison achieved the worst score, which is consistent with the results in [28]. This is because the accuracy of the rise time deteriorated as the amplitude of the input signal decreased. Surprisingly, the RPID method does not work as well as the digital charge comparison method. One reason for this is the approximate treatments in the algorithm implementation process, such as calculating the exponential function. Moreover, it was identified that the decay time of the fast component in the CsI(Tl) array depends on the type of incident particles (more precisely, on the average energy density deposited in the crystal) [22, 63] and even on the total energy of the particles [20, 21], while the slow component decay time does not. If this is the case, RPID may not be an accurate method because of the underlying assumption of constant decay time. However, this remains sufficient for the separation of γ-rays and light-charged particles.

Fig. 9Full-screen View

(Color online) 2-dimensions PID spectra for three methods: (a) rise-time comparison, (b) digital charge comparison, and (c) RPID. The Q_{Long Gate} and Q_{Short Gate} are the two integral values in the long and short gates as shown in Fig. 8. The N_s and N_f are the relative amplitude values of the slow and fast decay components of the CsI(Tl) crystal.Red dashed lines in (b) and (c) indicate that the spectra approach zero when both horizontal and vertical axis values are reduced. The rise-time comparison spectrum shows that rise time and amplitude of the waveform are weakly correlated and can be considered as independent. This indicates that the list of PID parameters in Table 3 is well defined owing to the nearly constant values of the PID parameters.

Fig. 10Full-screen View

Histograms of PID parameters including FoMs for three methods: (a) rise-time comparison, (b) digital charge comparison, and (c) RPID.

Overall, the aforementioned descriptions support the conclusion that the digital charge comparison method is preferred owing to its good performance and ease of FPGA implementation. Performance of the RPID method lags slightly; however, the complexity of this algorithm rules this option out. The rise-time comparison method had the worst performance among the three methods; however, the result remains acceptable, making it an alternative option to the digital charge comparison method.