Improvement of Crystal Identification Accuracy for Depth-of-Interaction Detector System with Peak-to-Charge Discrimination Method

In positron emission tomography (PET), parallax errors degrade spatial resolution. The depth of interaction (DOI) information provides the position in the depth of the scintillator interacting with the γ-rays, thus reducing parallax errors. A previous study developed a Peak-to-Charge discrimination (PQD), which can separate spontaneous alpha decay in LaBr3:Ce. Since decay constant of GSO:Ce depends on Ce concentration, the PQD is expected to discriminate GSO:Ce scintillators with different Ce concentration. In this study, the PQD-based DOI detector system was developed, which can be processed online and implemented in PET. A detector was composed of four layers of GSO:Ce crystals and a PS-PMT. The four crystals were obtained from both the top and bottom of ingots with a nominal Ce concentration of 0.5 mol% and 1.5 mol%. The PQD was implemented on the Xilinx Zynq-7000 SoC board with 8ch Flash ADC to gain real-time processing, flexibility, and expandability. The results showed that the mean Figure of Merits in 1D between four scintillators are 1.5, 0.99, 0.91 for layers between 1st–2nd, 2nd–3rd, and 3rd–4th respectively, and the mean Error Rate in 1D between four scintillators are 3.50%, 2.96%, 13.3%, and 1.88% for layers 1, 2, 3, and 4, respectively. In addition, the introduction of the 2D PQDs resulted in the mean Figure of Merits in 2D greater than 0.9 and the mean Error Rate in 2D less than 3% in all layers.


Introduction
Small animal positron emission tomography (PET) is used to visualize and study various biological processes in small laboratory animals such as mice and rats [1]. Common examples of the biological process include neuroreceptor imaging [2], metabolic measurement [3], cancer studies [4], gene expression [5], and drug development [6]. To depict the detailed anatomy of small animals, fine spatial resolution is required.
To improve the spatial resolution of PET, administering nuclides with a shorter range of positrons, using a fine pitch detector, and decreasing non-collinearity of annihilation gamma rays can be taken [1,7]. The two annihilation gamma rays have an angular deviation of 180 ± 0.25 • in water [7]. It causes the Gaussian blur in the reconstructed images and degrades spatial resolution. The larger the PET detector ring diameter, the larger the FWHM of the Gaussian blur [1,7]. That is the reason why a smaller ring diameter is preferred, however, it causes parallax errors and degrades the reconstructed image quality, especially at the edge of the ring. The depth of interaction (DOI) PET reduces parallax errors by using the depth information where the annihilation gamma rays interact with the crystal. Time-of-flight (TOF)

Electronics, Data Acquisition and Processing
Current-to-voltage conversion for the PS-PMT output signal is needed for the A/D conversion of the waveforms. Since the frequency response of the operational amplifier might have affected the waveform shape and thus the PQD performance, load resistance conversion was adopted instead of the trans-impedance amplifier (TIA). The larger the load resistance value, the larger the amplitude of the waveform, whereas the frequency response degrades, thus an appropriate load resistance value must be selected for the appropriate acquisition of waveforms. In this study, 10 kΩ load resistors were selected to ensure that the PS-PMT outputs do not exceed the voltage range of a later DAQ board and that the waveforms retain the frequency components required for the PQD. Similarly, if a weighted summing amplifier circuit was used to estimate the position of the luminescent scintillator on the PS-PMT photocathode, the PQD discrimination performance might be degraded due to the frequency response of the amplifier. Therefore, an anger-type logic circuit was employed for position estimation, and no operational amplifiers other than the ADC driver were used in this study. To keep the pulse height, a 10 Ω resistor chain was used. shows the overview image. The scintillators were more distant from the PS-PMT with higher Cedoped concentrations.

Electronics, Data Acquisition and Processing
Current-to-voltage conversion for the PS-PMT output signal is needed for the A/D conversion of the waveforms. Since the frequency response of the operational amplifier might have affected the waveform shape and thus the PQD performance, load resistance conversion was adopted instead of the trans-impedance amplifier (TIA). The larger the load resistance value, the larger the amplitude of the waveform, whereas the frequency response degrades, thus an appropriate load resistance value must be selected for the appropriate acquisition of waveforms. In this study, 10 kΩ load resistors were selected to ensure that the PS-PMT outputs do not exceed the voltage range of a later DAQ board and that the waveforms retain the frequency components required for the PQD. Similarly, if a weighted summing amplifier circuit was used to estimate the position of the luminescent scintillator on the PS-PMT photocathode, the PQD discrimination performance might be degraded due to the frequency response of the amplifier. Therefore, an anger-type logic circuit was employed for position estimation, and no operational amplifiers other than the ADC driver were used in this study. To keep the pulse height, a 10 Ω resistor chain was used.
The PQD performance mainly depends on the accurate peak voltage acquisition of the waveform. Since developed DAQ system does not write waveforms to files, waveforms from aforementioned anger-type logic circuit acquired in a preliminary experiment using an oscilloscope (WaveRunner 64xi, Teledyne Japan Corporation, Fuchu, Japan) is shown in Figure 2 for reference. The output signals of the PS-PMT contain a lot of highfrequency noise, therefore a low-pass filter (LPF) is essential. The sampling rate and the bandwidth are also the dominant factors in determining the performance of the PQD since the accuracy of peak voltage acquisition affects the deviation of Vp/Q values. The PQD The PQD performance mainly depends on the accurate peak voltage acquisition of the waveform. Since developed DAQ system does not write waveforms to files, waveforms from aforementioned anger-type logic circuit acquired in a preliminary experiment using an oscilloscope (WaveRunner 64xi, Teledyne Japan Corporation, Fuchu, Japan) is shown in Figure 2 for reference. The output signals of the PS-PMT contain a lot of high-frequency noise, therefore a low-pass filter (LPF) is essential. The sampling rate and the bandwidth are also the dominant factors in determining the performance of the PQD since the accuracy of peak voltage acquisition affects the deviation of Vp/Q values. The PQD needs approximately 100 MSPS to discriminate Ce concentration difference of the GSO:Ce scintillators.
A Cosmo-Z board (TokushuDensiKairo Inc., Sumida-ku, Tokyo, Japan) has an 8ch 12-bit 100 MSPS Flash ADC connected with a Xilinx Zynq-7000 SoC. The bandwidth of the analog front-end (AFE) circuit is 50 MHz, the input voltage range is ±0.5 Vpp, and the input impedance is 50 Ω. Since arbitrary signal processing can be performed online by rewriting the Zynq-7000 firmware, this board was used as the DAQ board in this study. The Programmable Logic (PL) and the software working on the Processing System (PS) were developed in-house to specialize in the detector for the PQD.
The overview of the readout system and signal processing flow are shown in Figure 3. Output signals from the PS-PMT were divided into four signals, X+, X−, Y+, and Y− for position estimation. The four waveforms were amplified by a preamplifier and A/D was converted by the Flash ADC. Then the digitalized waveforms were sent to the 14thorder finite impulse response (FIR) LPF implemented on the Programmable Logic to acquire accurate peak voltage. The cut-off frequency of the LPF was set at approximately 35 MHz. The frequency response of the FIR filter is shown in Figure 4. The coefficients of the FIR filter to be implemented on FPGA must be integers, therefore, they were scaled by 2 15 in this system. Since this made the outputs of the FIR filter 2 15 times larger than the inputs, the outputs of the FIR filter were shifted 15 bits to the right, and 12 bits were extracted from the least significant bit (LSB) side. The waveforms applied to the LPF were sent to a trigger detector and a PQD engine. The trigger detector outputted a selectable length of the trigger signal for the waveforms exceeding the threshold level. The PQD engine calculated a Vp, Q, and pedestal during the trigger and a parameter encoder encoded them, combined with a 44-bit timestamp. The size of an event packet was 20 bytes. A round-robin arbiter read the event packets of each channel sequentially and stored them in Block RAM. An AXI DMA transferred the event packets from Block RAM to the SDRAM, and Processing System sent them to the back-end PC via Gigabit Ethernet (GbE). Graphical User Interface (GUI) software running on the back-end PC communicated with the Cosmo-Z board over a TCP socket, decoded the measurement data, calculated Vp/Q values, and performed histogram analysis and visualization.  A Cosmo-Z board (TokushuDensiKairo Inc., Sumida-ku, Tokyo, Japan) has an 8ch 12-bit 100 MSPS Flash ADC connected with a Xilinx Zynq-7000 SoC. The bandwidth of the analog front-end (AFE) circuit is 50 MHz, the input voltage range is ±0.5 Vpp, and the input impedance is 50 Ω. Since arbitrary signal processing can be performed online by rewriting the Zynq-7000 firmware, this board was used as the DAQ board in this study. The Programmable Logic (PL) and the software working on the Processing System (PS) were developed in-house to specialize in the detector for the PQD.
The overview of the readout system and signal processing flow are shown in Figure  3. Output signals from the PS-PMT were divided into four signals, X+, X−, Y+, and Y− for position estimation. The four waveforms were amplified by a preamplifier and A/D was converted by the Flash ADC. Then the digitalized waveforms were sent to the 14thorder finite impulse response (FIR) LPF implemented on the Programmable Logic to acquire

Experimental Methods
The supply voltage for PS-PMT was set at −990 V and the output signals were acquired by the Cosmo-Z board with in-house firmware and software. The gain of the PS-PMT depends on the position on the photocathode, which is considered to affect the performance of PQD. To evaluate the PQD performance depending on the position of scintillators, 64 positions were assigned corresponding to an 8 × 8 array with a 2.7 mm pitch, and the scintillator block was shifted to each position for measurement. combined with a 44-bit timestamp. The size of an event packet was 20 bytes. A roundrobin arbiter read the event packets of each channel sequentially and stored them in Block RAM. An AXI DMA transferred the event packets from Block RAM to the SDRAM, and Processing System sent them to the back-end PC via Gigabit Ethernet (GbE). Graphical User Interface (GUI) software running on the back-end PC communicated with the Cosmo-Z board over a TCP socket, decoded the measurement data, calculated Vp/Q values, and performed histogram analysis and visualization.

Experimental Methods
The supply voltage for PS-PMT was set at −990 V and the output signals were acquired by the Cosmo-Z board with in-house firmware and software. The gain of the PS-PMT depends on the position on the photocathode, which is considered to affect the performance of PQD. To evaluate the PQD performance depending on the position of scintillators, 64 positions were assigned corresponding to an 8 × 8 array with a 2.7 mm pitch, The performance of the DOI detector was compared between collimated and noncollimated irradiation. For the non-collimated irradiation, the 123.9 MBq 137 Cs source was Sensors 2023, 23, 4584 6 of 15 placed 15 cm away from the detector and 40,000 events were collected. For the collimated irradiation, the 137 Cs source was placed at the far end of a 15 cm long tungsten collimator, 10,000 events were collected for each layer, and all layers were scanned. The experimental setup is shown in Figure 1b.
Jigs to position the scintillator block and the collimator were created by using a 3D printer (Anycubic Photon M3 Plus, Anycubic, Shenzhen, China).

Performance Evaluation
In accordance with the central limit theorem, 662 keV total absorption peaks of the energy histogram and Vp/Q values each follow a normal distribution. Since the Vp/Q value does not depend on deposit energy but on only the time constant of the scintillator [24], the 662 keV total absorption peak can be represented by a normal distribution of two variables, energy and Vp/Q value.
Energy resolutions for each layer were the FWHM obtained by fitting analysis with a one-dimensional normal distribution for the 662 keV total absorption peak of 137 Cs.
To determine the analysis range for the PQD and to evaluate the DOI discrimination performance, total absorption peaks were fitted with one-dimensional or two-dimensional normal distributions. Figure 5 shows the one-dimensional and two-dimensional histograms and their analysis ranges.
where and are the lower and higher boundaries for m-th crystal along x-axis, respectively, and are the lower and higher boundaries for m-th crystal along The boundaries of the analysis range for each crystal in the histograms are defined as where B lxm and B hxm are the lower and higher boundaries for m-th crystal along x-axis, respectively, B lym and B hym are the lower and higher boundaries for m-th crystal along y-axis, respectively, µ x and µ y are the mean values of x and y of the distribution, respectively, σ x and σ y are the standard deviations of x and y of the distribution, respectively, and I(a, b) is the intersection of the a-th and b-th normal distributions, which is the point between the two averages. The Vp/Q histograms for DOI discrimination performance were calculated by extracting events inside the boundaries (1)-(4). ER nD (Error Ratio) and FOM nD (Figure of Merit) were defined as n-dimensional DOI identifiability by fitting with n-dimensional normal distribution. ER nD , which is the discrimination ability of each crystal, is computed as where n EnD is the number of error events occurred at the other crystals and n TnD is the number of total events inside the boundaries. FOM nD , which is the discrimination ability between 2 crystals, is computed as where mean nDa is the mean value of one, mean nDb is the mean value of the other, FWHM nDa is the FWHM of one, and FWHM nDb is the FWHM of the other. Figure 6 shows an illustration of the one-dimensional ER and FOM.
where meannDa is the mean value of one, meannDb is the mean value of the other, FWHMnDa is the FWHM of one, and FWHMnDb is the FWHM of the other. Figure 6 shows an illustration of the one-dimensional ER and FOM. Calculations of the , and , are given in the Appendix A.

Results
A flood histogram acquired by the collimated irradiation is shown in Figure 7. It shows that the targeted layer is irradiated accurately. Some leakage events are observed in the layer on the photocathode side, however, few leakage events existed in the layer on the opposite side. Therefore, these leakages are due to the collimator diameter and experimental setup, not due to crosstalk. Calculations of the n EnD , n TnD and FW HM 2Da,b are given in the Appendix A.

Results
A flood histogram acquired by the collimated irradiation is shown in Figure 7. It shows that the targeted layer is irradiated accurately. Some leakage events are observed in the layer on the photocathode side, however, few leakage events existed in the layer While Figure 8 shows the results for collimated irradiation, Figure 9 shows the results for non-collimated irradiation, and Figure 10 shows the differences. The mean values and SDs of the performance index for each layer are summarized in Table 2.  While Figure 8 shows the results for collimated irradiation, Figure 9 shows the results for non-collimated irradiation, and Figure 10 shows the differences. The mean values and SDs of the performance index for each layer are summarized in Table 2.
The energy resolution of the first layer was worse than the other layers, approximately 18% on average. The other layers had an average energy resolution of 12.4-14.4%. The difference in mean values of energy resolution with and without collimation was approximately 2%.
The normalized gain was higher in the center of the photocathode or closer to the photocathode and lower at the edge of the photocathode or farther from the photocathode. There was little difference in normalized gain with or without collimation.
FOMs were higher and ERs were lower for crystals closer to the photocathode, which means that crystals closer to the photocathode have better DOI discrimination ability. Since each crystal layer has a different 662 keV peak position due to the difference in the amount of light reaching the PS-PMT, two-dimensional discrimination using both Q and Vp/Q values showed improved discrimination performances.
The differences in FOMs with and without collimation were slight, within 1.2%, however, the differences in ERs appeared to be large, up to approximately 54%. This is due to the very small values of ERs. means that crystals closer to the photocathode have better DOI discrimination ability. Since each crystal layer has a different 662 keV peak position due to the difference in the amount of light reaching the PS-PMT, two-dimensional discrimination using both Q and Vp/Q values showed improved discrimination performances.
The differences in FOMs with and without collimation were slight, within 1.2%, however, the differences in ERs appeared to be large, up to approximately 54%. This is due to the very small values of ERs.

Discussion
The DOI detector system based on the PS-PMT was developed by using the PQD. Less than 2.2% differences with and without collimation were shown for the mean values of detector performances other than ERs. Therefore, collimation does not affect the PQD discriminability.

Discussion
The DOI detector system based on the PS-PMT was developed by using the PQD. Less than 2.2% differences with and without collimation were shown for the mean val-ues of detector performances other than ERs. Therefore, collimation does not affect the PQD discriminability.
The first layer showed worse energy resolution than the other layers. Optical photons generated in each layer travel in two paths, one heading directly to the PMT, and the other traveling to the fourth layer and then being reflected back to the PMT. Photons generated in the first layer of the scintillator tend to have worse energy resolution due to the large difference in the length of these two paths. However, the energy resolution of the detector system developed in this study is comparable to previous studies [9,10,[12][13][14][15]17,18,20,22,23] and is considered sufficient for DOI-PET for small animals.
In this study, a 2D PQD was introduced that utilizes the unique property of the phoswich detector to have different normalized gains for each layer. Although the ERs 1D of the third layer was relatively poor, the 2D PQD resulted in the mean ERs 2D of less than 3% and the mean FOMs 2D of more than 0.9 for all layers. Further optimization of the decay constant and crystal position might improve the discrimination ability by increasing the 2D Mahalanobis distance in Q vs. Vp/Q 2D histogram. The optimization might include the material of the optical grease.
The higher the normalization gain, the higher the DOI discrimination performance was observed. This might be caused by the improved accuracy in obtaining Vp and smaller variation in Vp/Q as the gain increases. Using the larger value load resistor for current-tovoltage conversion was expected to increase the gain and to improve PQD performance, however, the input range of the Cosmo-Z and the scintillator position dependency of the gain of PS-PMT make changing the resistor value difficult. The input range of the Cosmo-Z is ±500 mV, and the pedestal is approximately 0 mV. The pulse height at the high gain position of the PS-PMT is about 400 mV, which is near the limit of the input range of the Cosmo-Z. If a larger load resistor is used, the waveform might be saturated. On the other hand, at the low gain position of the PS-PMT, such as the edge of the photocathode, the 662 keV total absorption peak is likely to be buried in noise due to low pulse height. Thus, using lower load resistor seems to be impractical. Therefore, designing a low-noise, wide-bandwidth trans-impedance amplifier that amplifies the PMT output to barely the input range of the Cosmo-Z might improve PQD performance. Optimization of the voltage divider circuit for position estimation might also contribute to improved performance. However, in this study, the scintillators were not arranged in an array. The scintillator array might have caused performance degradation due to crosstalk.
Since the PQD is a kind of PSD, the crystal shape is not constrained. Some previous studies [28,29] propose unique arrangements of the crystals such as taper shape to improve the detection efficiency. The detector system developed in this study could be implemented for such arrangements. Additionally, a previous study [23] combined PSD and LS to achieve four layer discrimination. Since PQD is a kind of PSD and is considered to be unaffected by reflective materials, PQD could also be combined with LS or laser engraving techniques. It might allow for more multi-layered discrimination.
The count rate when the source was placed close to the detector without the collimator was saturated at approximately 12 kcps. Since a single event was measured on four channels, the maximum event rate that can be processed by the developed system was approximately 3 k events/sec. The mass attenuation coefficient of GSO is 0.0814 cm 2 /g for 662 keV and the density of GSO is 6.7 g/cm 3 , thus the probability of interaction with 2.5 mm thickness GSO:Ce is 12.7%. The geometrical efficiency with the tungsten collimator was 6.95 × 10 −3 %. Therefore, the detection efficiency of collimated measurement was 8.82 × 10 −4 % and the event rate with 123.9 MBq 137 Cs source is 1 k events/sec. It means that counting losses might occur in PET applications with a higher count rate although all events could be measured for the collimated experiment in this study. The signal processing circuit implemented on the Programmable Logic was pipelined and running at 100 MHz, which is the same as the sampling frequency of the Flash ADC, and the event processing rate of the Programmable Logic was 5 M events/sec. Since the DMA bandwidth is also larger than the Gigabit Ethernet (GbE) data transfer bandwidth [30], GbE communication and postprocessing for visualization are expected to be the rate-limiting factor. The backend PC was equipped with 8GB of RAM and Intel Core i5-5278U CPU (2.9 GHz × 4 cores). The mean software processing time and SD were measured to be 99.07 ± 16.17 ms for socket communication, 37.94 ± 2.477 ms for decoding, 3.352 ± 10.09 ms for histogram calculation, and 97.90 ± 7.029 ms for visualization. A DAQ developed in a previous study has processed 250 k events/sec [31], and compared to this, the system developed in this study is slow. However, the processible data rate improvement is out of the scope of this study. The histogram calculation and visualization can be omitted for some applications, although socket communication and decoding are required. Thus, to achieve higher processible event rates, omitting post-processing for visualization, hard-coded communication, full FPGA implementation of PQD, event packet reduction, and improvement of back-end PC performance would be effective. These improvements should be taken for future work.
In this study, 662 keV gamma rays from 137 Cs source were measured. In detecting 511 keV gamma rays for PET applications, a large load resistor for current-to-voltage conversion might be needed since the lower pulse height than that of 662 keV gamma rays might reduce the PQD performance. LSO:Ce, which is one of the most commonly used scintillator for PET, has approximately four times higher light yield than GSO:Ce. Therefore, if LSO:Ce is used with the developed system in this study to achieve more layers of scintillators, gain adjustment might be required. The developed system could discriminate it since the decay time of LSO:Ce is approximately 40 ns.
In this study, the PQD-based DOI-PET detector system was developed. The admirable AFE circuit of the Cosmo-Z and the FIR filter implemented in the Programmable Logic showed good discrimination ability as DOI detector even using the 1D PQD. Also, the 2D PQD introduced in this study showed more accurate DOI discriminations. Therefore, more multi-layered detectors could be implementable. Since the array detectors were not used in this study, the effect of crosstalk could not be evaluated. Future work will evaluate the PQD performance with a multi-layered GSO:Ce array. This system could be applied to various PSD applications since the PQD is simple and has few restrictions on detectors.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
When the detector consists of the M layered Gd 2 SiO 5 (GSO) crystals, M one-dimensional normal distributions, f m (x), appear in the one-dimensional Vp/Q histogram. Additionally, M two-dimensional normal distributions, f m (x, y), appear in the two-dimensional histogram for Q and Vp/Q.
The boundary of each crystal in the histograms is defined as where B lxm and B hxm are the lower and higher boundaries for m-th crystal along x-axis, respectively, B lym and B hym are the lower and higher boundaries for m-th crystal along y-axis, respectively, µ x and µ y are the mean values of x and y of the distribution, respectively, σ x and σ y are the standard deviations of x and y of the distribution, respectively, and I(a, b) is the intersection of the a-th and b-th normal distributions, which is the point between the two averages. For one-dimensional normal distribution, the number of error events occurred at the other crystals n E1D and the number of total events inside the boundaries n T1D are computed as For two-dimensional normal distribution, n E2D and n T2D are computed as The probability density function of the two-dimensional normal distribution, P r (x, y), is P r (x, y) = 1 where ρ is the correlation between x and y.
For n-dimensional normal distribution, the probability density function is determined by the Mahalanobis distance d. For two-dimensional normal distribution, Mahalanobis distance is Since the square of Mahalanobis distance follows the chi-squared distribution with n degrees of freedom, the cumulative distribution function of two-dimensional normal distribution, F(d), is Thus, when HWHM 2D is expressed in Mahalanobis distance, d HW HM 2D , (A11) becomes d HW HM 2D = 2ln(2). (A13) From (A11) and (A13), the HWHM 2D for arbitrary angle ϕ is calculated as where Thus, FW HM 2D for arbitrary angle ϕ is FW HM 2D (ϕ) = 2 × HW HM 2D (ϕ).