Energy calibration of the silicon strip detectors

Strip-by-strip energy calibration is required for SSDs because of their fabrication limitations and the nonuniformity of the associated electronics. Due to the linear energy response of SSDs, the combination of the α source and pulser calibration method was used. Each telescope was equipped with an aluminum Mylar foil entrance window (2.0 μm Mylar + 0.06 μm Al) to protect the detector and provide electromagnetic shielding. For source calibration, a three-component α source comprising ²³⁹Pu, ²⁴¹Am and ²⁴⁴Cm with energies of 5.147, 5.480, and 5.795 MeV, respectively, was used. The energy deposited in the silicon strip detector was corrected according to LISE++ calculations, resulting in calibrated values of 4.904 MeV for ²³⁹Pu, 5.248 MeV for ²⁴¹Am, and 5.571 MeV for ²⁴⁴Cm. An interactive program was developed to fit the peak of the ADC channel (α) and to determine its relationship with E(α), as shown in Fig. 3.

Fig. 3Full-screen View

(Color online) Example energy spectrum of α particles fitting process (left panel) and results checking (right panel) from an individual strip on the single-sided silicon detector at SSDT1

Pulser calibration was used to ascertain the response characteristics of the electronics. An automatic peak-finding program was developed to obtain peak values, as shown in Fig. 4(a). Linear fitting [Eq. (1)] of the standardized pulse amplitudes (V) and the ADC channel (Ch) was carried out, as shown in 5(b). $V=k^{\prime }\times Ch+b^{\prime }$ (1) Assuming a linear relationship between V and the energy deposition (ε) through $V=\beta \epsilon$, the energy calibration was finalized, as shown in Eq. (2): $\epsilon =\frac{k^{\prime }}{\beta }\times Ch+\frac{b^{\prime }}{\beta }$ (2)

Fig. 4Full-screen View

(Color online) Pulser calibration results for peak-finding (a) and linear response fitting (b)

Particle identification

Particle identification is not only essential to calibrate the CsI(Tl) energy response to different types of incident particles, but also plays a crucial role for hits matching in the track reconstruction process. In the scatter plot ΔE - E, particles with different charges and masses tend to gather in different bands. Typically, PID is performed through graphical cuts around each band. Despite its versatility and ability to accommodate diverse scenarios, this approach suffers from two major limitations. First, in regions characterized by low particle statistics, it fails to extrapolate reliable identification, as there is insufficient data to establish clear cut boundaries. Second, manual contouring of each particle deposited in every Si-CsI(Tl) detection cell is a time-consuming and tedious process. As the number of detection cells increases, the time-consuming problem becomes more pronounced [22]. To address the aforementioned issues, Tassan-Got and Gawlikowicz et al. [23, 24] proposed a novel identification function based on the Bethe–Bloch formula [23, 24]. This function enables accurate modeling of the data and safe extrapolations in data-sparse regions by fitting the ridge lines of each isotope using 14 parameters, as shown in Eq. (3). $\begin{array}{ll}\text{Δ}E=[{\left(p_{1}E\right)}^{p_2+p_3+1}+{\left(p_4{Z}^{p_5}{A}^{p_6}\right)}^{p_2+p_3+1}\hfill & \hfill \\ {+p_7{Z}^2{A}^{p_2}{\left(p_{1}E\right)}^{p_3}]}^{\frac{1}{p_2+p_3+1}}-p_{1}E\hfill & (a)\hfill \\ E=p_{8}W(Z)L+p_{9}Z\sqrt{A}\ln \left(1+p_{10}W(Z)L\right)\hfill & (b)\hfill \\ W(Z)=p_{11}\frac{\ln \left(1+p_{12}Z\right)}{Z}\left(1-\exp \left(-p_{13}{Z}^{p_{14}}\right)\right)\hfill & (c)\hfill \end{array}$ (3) Here, L is the light output of particles in CsI(Tl) (equivalent to the ADC channel) and $W(Z)$ is a gain factor.

Figure 5(a) displays the ΔE-E scattering distribution for one hit events, in which the energy loss ΔE is from the front strips of the DSSSD transmission detector (layer 2) and the energy E is deposited in one CsI(Tl) unit of SSDT1. The sampling points on the isotopes' ridge curves were manually selected and indicated by the red dotted lines. The best fit was achieved with an individual set of parameters constrained by the simultaneous fitting of all visible isotopes. The mass and charge of each particle with given values of ΔE and E can be extracted by inverting the multiparametric identification function. When the extracted mass was not zero, the events were shown in Fig. 5(b). As defined in Ref. [23], the particle identification variable PID_n is calculated as ${\text{PID}}_n=Z_n+0.2\left(A_n-2Z_n\right)$, where n represents different particles. Figure 5(c) presents the ${\text{PID}}_n$ spectrum. By fitting the ⁴He peak to a Gaussian function, the relative mass resolution for the α particle can be derived.

Fig. 5Full-screen View

(Color online) Selected events with only one particle hitting in the DSSSD transmission detector (layer 2) and one CsI(Tl) unit of the SSDT1. Panel (a) shows the $\text{Δ}{E}_2-E_{\text{CsI}}$ scattering plot and the manually sampling points on the ridge curves (the dotted red line). Panel (b) displays the results when the extracted mass from inverting the 14-parameter identification function is not zero. Panel (c) is the PID_n spectrum

The correction of non-uniformity in the thickness of thin SSDs is of crucial importance for PID. Before correction, the isotopic resolution of the SSSSD-DSSSD component in SSDT2 was comparatively poor, as shown in Fig. 6(a). This is due to the random energy deposition when the particles pass through the detector at different positions with varying path lengths. For thickness correction, the sensitive area was divided into 8 × 8 bins. In each bin, we plotted the ΔE₁-E₂ histogram individually and selected a specific bin showing clear PID bands as a reference. Subsequently, for the remaining bins, the thickness correction factor η for the SSSSD was adjusted from -0.05 to 0.10 with a step size of 0.005 and the energy loss ΔE₁ in the SSSSD was corrected to ΔE₁(1+η). For each value η, a two-dimensional histogram was generated and compared with the reference histogram. Subsequently, the thickness correction for each bin was determined according to the value η, which provided the best consistency with the reference histogram. Figure 6(b) shows the two-dimensional graph after thickness correction for the 104 μm SSSSD, demonstrating a significant improvement in identification resolution compared to the graph before correction in Fig. 6(a). Figure 6(c) illustrates the “relative thickness variation” for the 104 μm SSSSD, ranging from -1% to 7.5%. This approach offers a relative correction for the thickness of a thin ΔE detector with respect to the reference bin. The absolute thickness correction relies on in-beam tests, as described in Refs. [25, 26].

Fig. 6Full-screen View

(Color online) Panel (a) represents a ΔE₁-E₂ scattering plot of SSDT2 before thickness correction for the 104 μm ΔE₁ SSSSD (the thickness of E₂ detector is 1010 μm). Panel (b) illustrates the two dimensional plot after the thickness correction for the ΔE₁ detector. Panel (c) displays the “relative” thickness correction for the ΔE₁ detector

After thickness correction of the SSSSD, each isotope can be effectively separated. Due to the impact of punching through on the scattering plot ΔE₁-E₂, the traditional method was used instead of the multiparametric formula. In Fig. 7(a), the distinct particles are identified by the green inclusive contour. Additionally, the center of each isotope ridge was fitted by a polynomial, $\begin{array}{c}{f}_n\left(E\right)=a_n^0\cdot E^{-1}+a_n^1\cdot E^0+a_n^2\cdot E^1+a_n^3\cdot E^2\\ +a_n^4\cdot E^3+a_n^5\cdot E^4+a_n^6\cdot E^5+a_n^7\cdot E^6\end{array}$ (4) where E is the energy deposited in the front strips of the DSSSD detector, $a_n^0$ to $a_n^7$ are eight parameters, and n represents the different particles as previously defined. The $\text{PID}\left(E_2\mathrm{,}\text{Δ}{E}_1\right)$ of each particle can be calculated using Eq. (5). Figure 7(b) and (c) display the results after identifying the particles. Table 3 summarizes the relative mass resolution of α particles for each telescope used in the experiment. $\begin{array}{c}\text{PID}\left(E_2\mathrm{,}\text{Δ}{E}_1\right)={\text{PID}}_n\\ +\frac{\text{Δ}{E}_1-f_n\left(E_2\right)}{f_{n+1}\left(E_2\right)-f_n\left(E_2\right)}\cdot \left({\text{PID}}_{n+1}-{\text{PID}}_n\right)\end{array}$ (5)

Fig. 7Full-screen View

(Color online) Particle identification of SSDT2 using the traditional method. In panel (a), distinct particles are identified by the green inclusive contour, and each isotope ridge is marked by a red dotted line. Panel (b) displays the results of events in the graphical cuts. Panel (c) shows the PID_n spectrum

Energy calibration of the CsI(Tl) crystals

Leveraging the high-quality calibrations of the DSSSD and the multiparametric formula obtained in the PID, the energy loss method can be utilized to perform the energy calibration on the CsI(Tl) crystal. Based on the two-dimensional plot for $\text{Δ}{E}_2-E_{\text{CsI}}^{\text{ch}}$, the range of $E_{\text{CsI}}^{\text{ch}}$ for each isotope was manually selected. For example, in Fig. 5(a), the range for protons is $90\le E_{\text{CsI}}^{\text{ch}}\le 3200$. Next, we selected 30 independent $E_{\text{CsI}}^{\text{ch}}$ points uniformly within this range. Using the multiparametric formula, we determine the corresponding energy loss (ΔE₂) in the DSSSD. Subsequently, the total kinetic energy of the incident particle Ek was deduced by numerically inverting Ziegler's energy loss table [27] based on ΔE₂ and the thickness of the DSSSD. Finally, the energy impinging on the CsI(Tl) crystal ($E_{\text{CsI}}^{\text{MeV}}$) was calculated by subtracting the energy losses in the SSD and Mylar foil that wraps the entrance face of the CsI(Tl) crystal, $E_{\text{CsI}}^{\text{MeV}}=E_k-\text{Δ}{E}_2-\text{Δ}{E}_{\text{Mylar}}$, which corresponds to the energy associated with the initial raw ADC channel $E_{\text{CsI}}^{\text{ch}}$.

Figure 8 shows the results of the energy calibration of the CsI(Tl) crystals for each isotope. The CsI(Tl) energy response for the incident particle exhibited insignificant nonlinearity and a slight difference between the isotopes in the energy range of 0 to 100 MeV, in accordance with the results in Ref. [28]. Two separate formulas were used to describe the calibration process. For Z = 1 isotopes, the calibration function was written as [28] $\begin{array}{l}L(E\mathrm{,}Z=\mathrm{1,}A)=a_0+a_1{E}^{\left(a_2+A\right)/\left(a_3+A\right)}\hfill \end{array}$ (6) where a₀ is an offset, a₁ is a gain factor, A is the mass number of Z = 1 isotope, and a₂, a₃ are empirical non-linearity parameters. The fitting results are presented in Fig. 8(a), where three solid lines are depicted: the red line represents the proton, green line represents the deuteron, and purple line represents the triton. For isotopes with Z ≥ 2, the standard Horn's formula was used for calibration [29], $\begin{array}{l}L(E\mathrm{,}Z\ge \mathrm{2,}A)=a_0+a_1\left(E-a_{2}A{Z}^2\text{ln}\left|\frac{E+a_{2}A{Z}^2}{a_{2}A{Z}^2}\right|\right)\hfill \end{array}$ (7) where a₀, a₁ and a₂ are the parameters obtained from the simultaneous fitting of all heavy isotopes. The perfect fitting results are shown in Fig. 8(b).

Fig. 8Full-screen View

(Color online) Energy calibration of one CsI(Tl) unit in SSDT1. Panel (a) shows the results for hydrogen isotopes: proton (red points), deuteron (green squares), and triton (purple triangles). The corresponding curves represent the fitting results with Eq. (6). Panel (b) displays the calibration result for heavy ions ^3,4,6He, and ^6,7Li. The corresponding curves represent the fitting results with Eq. (7)

Constraints analysis

The geometric relationships between different layers are determined solely by the structure of the SSDTs, which is related to the design and installation of the detectors. A valid track with hits in each detector layer should approximately be a straight line. Table 4 provides a detailed geometric mapping of the CSHINE SSDT. For particles that are stopped in a specific unit of the CsI(Tl) crystal, denoted by N_3A, only certain strips of DSSSD, marked by N_2F and N_2B for the front and back strips, respectively, satisfy the geometric relationships G_3A,2B and G_3A,2F. The incident event of one particle can be used to confirm geometric relationships, as illustrated in Fig. 9. On the surface of each SSD, there were 32 strips with a width of 2 mm, and the gap between adjacent strips was 100 μm. As in the previous experiment, to save electronic channels, every two strips are merged into one channel, and thus, there is a total number of 16 channels on each side. However, an improved merging technique was implemented, where two strips near the edge (first and thirty-second strips) were merged into one channel to minimize the shadow effect in the edge channels. Furthermore, geometric matching between SSSSD strips and DSSSD back strips requires G_2B,1S. According to the configuration of the CSHINE SSDT, the SSSSD and the back strips of the DSSSD are in parallel and numbered in the same direction. Because the distance between the SSDs is relatively shorter compared to the distance from the target, hits of a single track on the SSSSD and the backside of the DSSSD will occur on the same or neighboring strip numbers, as restricted by G_2B,1S of $|N_{\text{2B}}-N_{\text{1S}}|\le 1$.

Fig. 9Full-screen View

(Color online) Panels (a) and (c) display scattering plots before applying geometric constraints for one particle incident events in SSDT2, showing the geometric relationship between the front and back sides of the L2 DSSSD and L3 layer CsI(Tl) crystals. Panels (b) and (d) display scattering plots after applying geometric constraints

Table 5 lists the energy constraints that should be met for the particles stopped in DSSSD or CsI(Tl). During the particle identification process, the relationship between ΔE and E was determined. For DSSSD-CsI, the multiparametric identification function can be used to determine the mass of a particle using the information of the deposited energies (E_3A,Δ E_2F). The extracted mass should not be zero, denoted as $E_{\text{3A,2F}}^{\text{iso}}$. For SSSSD-DSSSD, the deposited energies of the particle (E_2F,Δ E_1S) should fall into the banana-shaped graphical cut corresponding to a specific isotope, which is represented by $E_{\text{2F,1S}}^{\text{iso}}$. Another important constraint is the relationship between the energies detected on the front (ΔE_2F) and back sides (ΔE_2B) of the DSSSD, indicated by E_2B,2F, and utilized for both the DSSSD-CsI and SSSSD-DSSSD cases. Figure 10 presents a processing example demonstrating how to determine the constraints of E_2B,2F for SSDT1. Events of single particle incidence, where the particle stops at CsI(Tl) and is constrained by three geometric constraints, were selected to investigate the relationship between ΔE_2F and ΔE_2B. This selection was performed to avoid the influence of inter-strip events. Figure 10(a) shows that Δ E_2F and ΔE_2B are approximately equal. To improve the precision of the constraint, the relative difference $\left(\text{Δ}{E}_{\text{2B}}-\text{Δ}{E}_{\text{2F}}\right)/\text{Δ}{E}_{\text{2F}}$ is defined and its absolute value should be less than 0.2, as shown in Fig. 10(c). Finally, events that stop at CsI(Tl) can be constrained by the energy loss ratio ΔE₁/ΔE₂ between the two silicon detectors. Based on the thicknesses of the SSDs, the energy loss of the particles as they penetrate the two SSDs was calculated using LISE++. The comparison between the calculated values and the experimental results is shown in Fig. 11. The left panel is a scattering plot of ΔE_1S-Δ E_2F, which confirms the reliability of the calculated values. The right panel displays the distributions of ΔE₁/ΔE₂, where the black and green curves represent the experimental and calculated results, respectively. In this study, only the calculated and experimental distribution ranges of ΔE₁/ΔE₂ were utilized as an additional criterion for track recognition. However, the difference in shape between the calculated and experimental spectra, which stems from the simplified assumption of a uniform incident energy distribution used in the simulations, was not taken into account. Taking SSDT1 as an example, the energy constraint E_2F,1S was determined as $0.10\le \text{Δ}{E}_{\text{1S}}/\text{Δ}{E}_{\text{2F}}\le 0.40$.

Fig. 10Full-screen View

(Color online) (a) The correlation between the energies detected in the front- (ΔE_2F) and back-sides (ΔE_2B) of the DSSSD in SSDT1. (b) The relative difference between ΔE_2F and ΔE_2B as a function of ΔE_2F. (c) The spectrum of relative differences in a logarithmic scale

Fig. 11Full-screen View

(Color online) Panel (a) displays the experimental scattering plot of ΔE_1S-Δ E_2F compared with the calculations using LISE++. In panel (b), the distributions of ΔE₁/ΔE₂ are presented. The black curve denotes the experimental result, while the green curve represents the calculated result

Track decoding

Two types of track decoding start from DSSSD and CsI(Tl), respectively. By considering the multiplicities extracted from each detector layer and matching them in a loop, all the possible cases can be comprehensively covered. Screening out candidate tracks through geometric constraint is crucial, as this process can eliminate a significant number of false tracks.

Table 6 lists some typical track modes that display charge-sharing and multi-hit effects, providing the basic characteristics for decoding complex track modes. In the track decoding algorithm, it is reasonable to assume that charge sharing only occurs between two adjacent channels in the same layer, denoted by $G_{\text{O}}^{{\text{sharing}}_{ij}}$, where i and j are the candidate track numbers and the subscript O represents the layer. In the cases of SSSSD-DSSSD, we exemplify the charge-sharing effect by taking one Tr2-211 event. The sum of the energy losses in the two adjacent strips, strip 11 and strip 12, in the “2B” layer has been shown to be approximately equal to the energy loss in the single strip 11 on the “2F” layer. This result meets the energy constraint $E_{{\text{2B}}_{{\text{sum}}_{12}}\mathrm{,}\text{2F}}$. If $E_{{\text{2F}}_1\mathrm{,}{\text{1S}}_1}^{\text{iso}}$ is also true, then one real track is finally determined. The track decoding flow diagram for this mode is shown in Fig. 12(a). In Table 6, Tr2-212 represents an event in which two valid candidate tracks exist. In this event, the two particles strike the same strip on the “2F” layer. In such cases, E_2B is used instead of E_2F for particle identification when double hits occur on the front side of the DSSSD. Both tracks must satisfy the energy constraint $E_{{\text{2F}}_i\mathrm{,}{\text{1S}}_i}^{\text{iso}}$ simultaneously, denoted as $E_{{\text{2B}}_{{\text{scale}}_{12}}\mathrm{,}{\text{1S}}_{12}}^{\text{iso}-{\text{2Hit}}_{12}}$. Fig. 12(b) shows the complete decoding process of Tr2-212, which also needs to be checked for the inter-strip effect if the multi-hit energy constraint mentioned above is not satisfied.

Fig. 12Full-screen View

(Color online) Schematic flowchart of the track decoding process for modes TrM2-211 [panel (a)] and TrM2-212 [panel (b)] in the L1L2 track reconstruction

In the DSSSD-CsI case, considering track mode Tr2-1121^a as an example, which includes the inter-strip effect that occurs in the “2F” layer, the sum of the energy losses in the two neighboring strips 12 and 13 in the “2F” layer and the energy deposited in the crystal should satisfy the energy constraint $E_{{\text{3A}}_1\mathrm{,}{\text{2F}}_{{\text{sum}}_{\mathrm{1,2}}}}^{\text{iso}}$. In addition, two constraints, $G_{\text{2F}}^{{\text{sharing}}_{12}}$ and $E_{{\text{2B}}_1\mathrm{,}{\text{2F}}_{{\text{sum}}_{12}}}$, were applied to guarantee the correct determination of the energy of the inter-strip event. Moreover, to complete the reconstruction of the track, we must cross-check the energy loss ΔE_1S in the “1S” layer. This is achieved using the energy constraint $E_{\text{2F,1S}}$. If this constraint is not satisfied, ΔE_1S can be recalculated by the following procedures: (1) the charge (Z) and mass (A) of the particle can be obtained from $E_{{\text{3A}}_1\mathrm{,}{\text{2F}}_{{\text{sum}}_{\mathrm{1,2}}}}^{\text{iso}}$; (2) the energy of the particle impinging the DSSSD, $E_2^{\text{imp}}$, can be deduced by the LISE++; (3) the energy of the particle impinging the SSSSD, $E_1^{\text{imp}}$, can also be deduced similarly; (4) the particle's energy loss in the SSSSD is finally obtained with $E_2^{\text{imp}}-E_1^{\text{imp}}$. Figure 13(a) shows the flow diagram of the Tr2-1121 track mode.

Fig. 13Full-screen View

(Color online) Schematic flowchart of the track decoding process for modes TrM2-1121 [panel (a)] and TrM2-2121 [panel (b)] in L2L3 track reconstruction

When events involve two or more particles hitting the same CsI(Tl) crystal, the real deposited energy in the crystal cannot be determined accurately, and these events are typically discarded. In the case of Tr2-2121^a in Table 6 as an example of multiple hits, two candidate tracks can satisfy the energy constraint $E_{{\text{3A}}_i\mathrm{,}{\text{2F}}_i}^{\text{iso}}$, and the sum of the energy loss in the two strips in the “2F” layer is equal to the energy loss in a single strip in the “2B” layer. This indicates that two particles hit the same strip in the “2B” layer simultaneously. Furthermore, since both particles also hit the same strip in the “1S” layer at the same time, their energy must be calculated using LISE++. The Tr2-2121 mode tracking decoding process follows the flow diagram shown in Fig. 8(b). In most cases, the candidate tracks do not satisfy the charge-sharing or multi-hit constraints, but are still considered valid tracks, as demonstrated by Tr2-1121^b and Tr2-2121^b in Table 6. These tracks can be found with $E_{{\text{2B}}_i\mathrm{,}{\text{2F}}_i}$, $E_{{\text{3A}}_i\mathrm{,}{\text{2F}}_i}^{\text{iso}}$ and $E_{{\text{2F}}_i\mathrm{,}{\text{1S}}_i}$ directly. When dealing with a larger value of M, each pair of tracks must be tested to determine the positions of the hits. Taking Tr4-2222 in Table 6 as an example, it is clear that candidate tracks No. 2 and No. 3 are two valid tracks. However, in the scenario where two particles with nearly identical energies strike the DSSSD, it becomes extremely challenging to identify the actual tracks. All decoding diagrams for the track modes with $M\le 4$ were systematically compiled and are provided as supplementary material to this article [21].

wTrack recognition efficiency is defined as the ratio of the number of successful track recognition events to the number of candidate events that meet the geometric constraints. In this study, we introduced reconstruction information, such as the final effective identification event number and energy spectra, into the track recognition algorithm to optimize the detector signal threshold. With this improvement, track recognition efficiency was further enhanced compared to our previous work [6]. As shown in Table 7, the track recognition efficiencies of all telescopes approached 90% in this experiment.