A study on “position-energy” response correction method based on monolithic crystal coupled SiPM array

In order to investigate this issue, the paper constructs a model of a monolithic scintillator detector using Geant4 software, as depicted in Fig. 1. The crystal employed in the model is a monolithic LaBr₃(Ce) crystal, known for its exceptional properties including density, energy resolution, light yield, and decay time. Compared to other commonly used scintillators, LaBr₃(Ce) crystals possess significant advantages. The specific properties of the LaBr₃(Ce) crystal are provided in Table 1. It is important to note that the simulation study conducted in this paper is also applicable to other types of monolithic scintillator detectors.

Fig. 1.

The monolithic LaBr₃(Ce) crystals of various sizes were employed in the detector model to investigate the relationship between the “position-energy” response of the detector and different detector sizes. The sizes used were 51 mm × 51 mm × 5 mm, 51 mm × 51 mm × 8 mm, 51 mm × 51 mm × 11 mm, and 102 mm × 102 mm × 5 mm. The LaBr₃(Ce) crystal was coated with a Teflon foil (0.3 mm thickness) which acts as the front diffuse reflective layer. The scintillator was surrounded and encapsulated by an aluminum foil (0.5 mm thickness) for the specular reflective layer [24]. A glass layer (3 mm thickness, reflectivity 1.5) and optical grease (0.1 mm thickness, reflectivity 1.41) were placed between the bottom of the crystal and the SiPM array. The SiPM utilized in the model was the Micro-30035-TSV model manufactured by ONSEMI, featuring a 16 × 16 array with a sensitive area of 3.07 × 3.07 mm² and a package area of 3.16 × 3.16 mm².

In this simulation study, the process of transporting visible photons can be divided into two main parts. The first part involves the transportation of visible photons within the medium, and the user can select the material parameters. These material parameters primarily include the emission spectra, scintillation yield, intrinsic resolution, and fast component of the scintillator. They play a crucial role in determining the production of visible photons and their transport through the medium. The second part involves the transportation of visible photons between different media, requiring the user to specify the surface parameters. The primary surface parameters include type, finish, efficiency, reflectivity, and reflectivity type. These parameters control the physical interaction of visible photons at the surface between different media.

The simulation program primarily utilizes the modular physics process FTFP_BERT, which involves replacing the electromagnetic interaction process with the standard electromagnetic process G4EmStandardPhysics_option1. Additionally, the program incorporates the optical process G4Opti-calPhysics. G4EmStandardPhysics_option1 encompasses various fundamental electromagnetic interactions, such as the photoelectric effect, Compton scattering, pair production, Coulomb scattering, and bremsstrahlung. On the other hand, the optical processes implemented in the program include scintillation, Cherenkov radiation, Rayleigh scattering, Mie scattering, refraction, and reflection.

In this study, the effective detection area of the monolithic crystal detector was divided into 25 equal blocks, as illustrated in Fig. 2. Initially, a radioactive source was positioned 180 mm away from the detector plane, and an equal number of parallel beams were emitted along the receiving surface of the detector, starting from the upper right corner of area 0. The signals from the SiPM array were summed to obtain the energy response signal of the detector when radiation was incident from a specific area into the detector. Next, the energy response of each region was determined by irradiating each region with different energy ray beams. This process yielded individual γ-energy spectra for the 25 incident regions. Table 2 provides the details of individual γ-energy spectra and corresponding emitted ray energies. Subsequently, the peak positions of the energy spectra for each of the 25 incident regions were extracted. The peak position of the energy spectrum from incident region 12 was selected as the standard peak position. A linear relationship matrix was then established, relating each region to the peak position of the corresponding γ-energy spectrum from incident region 12. Finally, using the peak position of the energy spectrum from incident region 12 as the standard, the remaining regions were corrected to align their peak positions with that of incident region 12. This correction process ensured that each region’s energy spectrum was consistent with the peak position of incident region 12. As a result, the corrected detector energy spectrum response was obtained.

Fig. 2.

In this simulation study, γ-rays with 10 common energies were emitted from a distance of 180 mm away from the detector. These radiation energies are listed in Table 2. Each energy was randomly sampled in equal proportions, and the vertical beams of these 10 energies were used to irradiate each region of the detector. As a result, the corresponding γ-energy spectra were obtained for each incident region of the entire detector, as depicted in Fig. 3.

Fig. 3.

From Fig. 3, it is evident that the energy response of the detector is satisfactory for different detector sizes, with 10 distinct peaks observed in each energy spectrum. This indicates the feasibility of the detector model developed in this study. However, in the same detector, the peak positions of the energy spectrum in the incident region 0 exhibit a noticeable shift compared to region 12, particularly at energies above 365 keV. This shift can be attributed to several factors. Firstly, when γ-rays interact with the crystal, they ionize and excite the atoms within the crystal. During the subsequent de-excitation process, photons are generated, and subsequently visible light is absorbed, refracted, and reflected within the crystal. The edge position of the energy spectrum is more susceptible to the effects of backscattering. Secondly, the photons within the crystal are influenced by the crystal’s structure and physical properties, which can affect the direction of their propagation and result in light losses. These factors contribute to spatial inconsistencies in the detector’s intrinsic properties, leading to a position-dependent energy response to γ-rays.

As the thickness of the crystal increases, the detector acquires the capability of depositing higher energies, and the impact of backscattering and refraction on the energy spectrum decreases. Consequently, the energy spectrum becomes less influenced by the position-dependent effects as mentioned earlier. This is supported by the observation that as the crystal thickness increases, the shift in peak positions becomes less pronounced.

The simulation study takes the energy spectrum peaks from the 25 incident regions of detectors with different sizes in order to address the significant peak shift observed in the edge region for medium-high energy (Eγ ≥365 keV) peaks. The energy spectrum peak in the incident region 12 is considered as the standard peak. Using Eq. (1), the difference between the energy spectrum peaks in the other 24 incident regions and the standard peak is calculated and summed up to obtain the sum of offsets for each region, and the obtained results are presented in Fig. 4. (1) Zi=∑j=1j=5| Y12−Yi | \[{{Z}_{i}}=\sum\nolimits_{j=1}^{j=5}{\left| {{Y}_{12}}-{{Y}_{i}} \right|}\] where Z_i is the sum of the peak offsets in the incident region i, Y₁₂ is the energy spectrum peaks above 365 keV in the incident region 12, Yi is the energy spectrum peaks above 365 keV in incident region i, and j = 1~5 is the energy spectrum peaks above 365 keV in each of the five incident regions.

It can be seen from the results that the peak inconsistency caused by different incident positions is common in the measurement of monolithic crystal detectors and the amount of peak offset varies for each incident region. According to Fig. 4(a), the 51 mm × 51 mm × 5 mm detector, for example, shows a small offset at the center and a large offset at the edges, which reflects a proportional relationship between the offset and the distance from the center. The more marginal the peak shift, the more severe it is. This is because γ-rays are more likely to backscatter and refract with backscattering materials when interacting at the edge of the crystal, and secondly, the physical structure and nature of the crystal also affects the propagation of photons within the crystal, resulting in a loss of ray energy, incomplete collection, and a positive relationship between offset and distance from the center. The other three sizes of detector offsets have the same distribution law, the phenomenon of a positive relationship between the offset and the distance to the center point. The distribution law also reflects the fact that when the crystal thickness is constant, the sum of the energy spectrum peak offset in the same detection region decreases as the crystal area increases; when the crystal area is constant, the sum of the energy spectrum peak offset in the same detection region decreases as the crystal thickness increases.

Fig. 4.

When γ-rays interact with the crystal, several processes such as scattering and refraction occur. Additionally, the properties of the SiPM in terms of light absorption and transport contribute to the characteristics of the γ-energy spectrum obtained by the detector. These factors result in several effects, including a shift in peak position, poor consistency of peaks, and compression of the full energy peak at the edge of the detector position. As a consequence, the energy resolution of the detector is compromised, as depicted in Fig. 5.

Fig. 5.

Figure 5 demonstrates that the energy resolution is very poor in the edge region compared to the center region for different detector sizes. For instance, when the detector area is 51 mm × 51 mm with a thickness of 5 mm, the energy resolution of the 662 keV characteristic peak in the incident region 0 is 4.58%, while in the incident region 12, it is 3.92%. Similarly, when the detector area is 102 mm × 102 mm with a thickness of 5 mm, the energy resolution of the 662 keV characteristic peak in region 0 is 4.18%, and it is 3.92% in region 12. The larger crystal area leads to higher detection efficiency and better energy resolution, resulting in the former having a poorer energy resolution than the latter. When the detector area is 51 mm × 51 mm with a thickness of 8 mm, the energy resolution of the 662 keV characteristic peak in region 0 is 4.62%, and it is 4.31% in region 12. Moreover, when the detector area is 51 mm × 51 mm with a thickness of 11 mm, the energy resolution of the 662 keV characteristic peak in region 0 is 4.53%, and it is 4.26% in region 12. The latter case exhibits better energy resolution than the former because thicker crystals allow for more complete energy deposition of the rays, thereby improving energy resolution within a certain thickness range where the rays scatter and refract with the crystal. Consequently, due to the mentioned factors, as the thickness and area of the detector crystal increase, the detected ray energies become closer to the energies of the emitted rays.

In this study, a “position-energy” response model was established to analyze the γ-energy response spectra of the 25 incident regions in the entire detector using crystals of different sizes. Here the focus of correction was on the peak positions of energy spectra above 365 keV. The peaks above 365 keV were extracted from each energy spectrum, with the peaks in the incident region 12 serving as the reference standard peaks. It was found that the peak position of the energy spectrum in each region of the same size detector is linearly correlated to the peak position at the center y = kx + b. For instance, in the case of a 51 mm × 51 mm × 5 mm crystal, the peak level of the energy spectrum (Eγ >365 keV) in the incident region 0 exhibited a linear correlation with the peak level of region 12, as depicted in Fig. 6.

Fig. 6.

Figure 6 exhibits the linear relationship between the peak levels of the energy spectrum in incident region 0 and the peak levels of the energy spectrum in the incident region 12. The linear relationship is y = 1.11x + 3.07, where k = 1.11, b = 3.07, and R²(COD) = 0.9999. Therefore, this simulation study establishes the correction function matrices with respect to different detector sizes separately. A binomial parametric fitted straight line is established between the energy spectrum peak level in region i (Eγ >365 keV) and its corresponding detector peak level in the incident region 12 (Eγ >365 keV), and the linear relationship is expressed as: (2) Y=ki×Xi+bi \[Y={{k}_{i}}\times {{X}_{i}}+{{b}_{i}}\] where Y is the energy spectrum peak position in the incident region 12, X_i is the peak position in the incident region i, k_i is the slope in region i, and b_i is the intercept in the incident region i.

The position and energy information of the incident particle is obtained according to the positioning algorithm. When the incident ray is irradiated in region i, the (k, b) of the corresponding region is matched according to the established correction function matrix, and the peak position of the corresponding region is corrected, thus correcting the whole detector energy spectrum. (3) [ (k11,b11)…(k15,b15)(k21,b21)(k25,b25)(k31,b31)⋱(k35,b35)(k41,b41)(k45,b45)(k51,b51)…(k55,b55) ] \[\left[ \begin{array}{*{35}{l}} \left( {{k}_{11}},{{b}_{11}} \right) & \ldots & \left( {{k}_{15}},{{b}_{15}} \right) \\ \left( {{k}_{21}},{{b}_{21}} \right) & {} & \left( {{k}_{25}},{{b}_{25}} \right) \\ \left( {{k}_{31}},{{b}_{31}} \right) & \ddots & \left( {{k}_{35}},{{b}_{35}} \right) \\ \left( {{k}_{41}},{{b}_{41}} \right) & {} & \left( {{k}_{45}},{{b}_{45}} \right) \\ \left( {{k}_{51}},{{b}_{51}} \right) & \ldots & \left( {{k}_{55}},{{b}_{55}} \right) \\ \end{array} \right]\]

Based on the correction function matrix obtained from the simulation study for different sizes of crystal detectors, the energy spectrum (Eγ >365 keV) peak positions were corrected for each incident region at different sizes, and the correction results are shown in Fig. 7.

Fig. 7.

Figure 7 clearly demonstrates the reduction in the offset of the energy spectrum peak positions achieved through the correction function matrix, thus highlighting its effectiveness in correcting the energy spectrum peak positions for different detector sizes. The results indicate significant improvements in offset reduction and are as follows:

For the 51 mm × 51 mm × 11 mm detector (Fig. 7a), the total offset is reduced from 48.43 to 17.48.

Comparing the results of Fig. 7(a), (b), (c), and (d), it is evident that the correction function matrix is particularly effective in correcting areas with larger offsets, with the greatest effectiveness observed in the 51 mm × 51 mm × 5 mm detector.

In this simulation study, the Cs-137 point source is utilized as an illustrative example. Using the correction function matrix derived from the simulation study, the peak positions of the energy spectrum in different incident regions are individually corrected. This correction process yields the energy spectrum of the Cs-137 point source, as well as the energy resolution at 662 keV before and after the overall detector correction. Figure 8 exhibits the visual representation of the obtained results.

Fig. 8.

As can be seen in Fig. 8, the detector’s own matching correction function matrix is universally applicable, and the correction function matrix not only corrects the energy spectrum peak position well, but also improves the energy resolution of the detector with a Cs-137 point source of 662 keV. So, for a Cs-137 point source of 662 keV:

The energy resolution of the 51 mm × 51 mm × 5 mm detector can be improved from 4.54% to 3.92% after the peak correction;

Notably, the correction function matrix exhibits a weaker correction effect for positions with smaller peak offsets. This can be attributed to the fact that larger crystal sizes result in more efficient energy absorption by the crystal during ray interactions, leading to energy spectrum peak positions that are closer to the simulated ray energy.

The linearity of the energy of γ-rays and the pulse amplitude characterizes the energy linearity index of the γ-spectrometer [25], and the linearity directly affects the determination of the energy of γ-rays and the ability and accuracy of the identification of nuclide species by the spectrometer. Generally speaking, the relationship between the energy of γ-rays and the address of the all-energy peak height channel can be expressed as a linear function, which is (4) E=A+Bh \[E=A+Bh\]

In this simulation study, the Monte Carlo procedure emits rays in the energy range 0–1332 keV, therefore this simulated detector energy linear function relationship is only applicable in the range 0–1332 keV. In the equation, E (keV) is the energy emitted by the simulation, h (keV) is the energy of the rays collected by the simulated detector, and the effect of this correction method on the energy linearity is determined by comparing the deviation from the energy linearity of different detector sizes before and after the correction. Table 3 shows a comparison of the deviation from energy linearity before and after correction for different detector sizes.

According to Table 3, it is clear that the peak energy corrected by the correction function matrix is closer to the energy of the emitted rays. By comparing the energy linearity of the four different detector sizes, it can be concluded that the method is effective in improving the energy linearity match of the detectors, resulting in a significant reduction in the linearity deviation. For incident ray energies in the range 0–1332 keV, the linearity deviation can be reduced from at least 1.6% to 1.4%, reducing the linearity error by 0.2%, and the linearity deviation can be reduced from at most 2.1% to 1.2%, shrinking the linearity error by 0.9%.