2. Radiography System Optimization

2.1. Proton Radiography System

The basic configuration of the proton radiography system consists of three key parts: the diffuser, matching lens, and Zumbro lens. The proton beam whose direction is from left to right in system is shown in Figure 1. A pencil mono-energy proton beam is firstly scattered in a diffuser. With a broader angular distribution, it transfers in to matching lens. Both the angle and size of the beam are modified in matching lens according to the coordinate-angle correlation required by Zumbro lens. Then, the protons with this correlation can be imaged by Zumbro lens when it reaches the image plane as shown in Figure 1.

Figure 1: Layout of proton beam in radiography system.

2.2. Zumbro Magneto-Optical Lens

The MCS and proton-nucleus elastic scattering result in a random direction of the protons as they leave the object. This effect blurs the recorded image. In addition to this, due to the stopping of object material, the proton energy spreads when it exits the object. This energy spreading leads to chromatic aberration when the proton transfers in the magnetic lens. The chromatic aberration finally blurs the recorded image. The Zumbro lens can remove the blur caused by MCS, as well as the first order of chromatic aberration. According to Ref. [Reference Mottershead and Zumbro1], the rest second-order chromatic aberrations term Err_c has the following form:

(1) $\begin{array}{}Er{r}_c=m_{12}^{\prime }\cdot \varphi \cdot \sigma ,\end{array}$

where ϕ is the angular deviation from the ideal position-angle correlation line and σ is the fractional momentum deviation. Both ϕ and σ scale inversely with the beam momentum. Therefore, the higher the proton energy, the lower the chromatic aberrations will be. Once the proton beam and object are set, Err_c is proportional to the value of $m_{12}^{\prime }$ , regardless of collimation and angle cut. In this case, better resolution is therefore expected with lower $m_{12}^{\prime }$ . The optimization of Zumbro lens is to minimize the value of $m_{12}^{\prime }$ .

2.2.1. Goal of Optimization

As shown in Figure 1, Zumbro lens has two same subsections. Each one consists of eight quadrupoles. In the simplest case, the transfer matrix MZ of one section can be described [Reference Jian-qin21] by four parameters, quadrupole length (L/m) and strength (G/tesla ^∗m⁻¹), the inside drifts (D/m), and the outside drifts (S/m):

(2) $\begin{array}{}MZ=\left(\begin{array}{cc}{m}_{11}\left(L,G,S,D\right)& m_{12}\left(L,G,S,D\right)\\ m_{21}\left(L,G,S,D\right)& m_{22}\left(L,G,S,D\right)\end{array}\right).\end{array}$

Because the expression of $m_{12}^{\prime }$ is an implicit function, an approximate calculation is used that $m_{12}^{\prime }=\left(m_{12}\left(L,G,S,D,E\right)-m_{12}\left(L,G,S,D,E-\Delta E\right)\right)/\Delta E$ . The goal of optimization is to minimize $m_{12}^{\prime }$ .

2.2.2. Optimization Constraints

During the searching of the best parameter set, several constraints should be taken into consideration. The first one is unit-magnification point-to-point inverting identity lens constraint for both x or y coordinate plane (beam along z axis). So, this constraint imposed on the optimization can be expressed by the following:

(3) $\begin{array}{}MZ=\left(\begin{array}{cc}{m}_{11}\left(L,G,S,D\right)& m_{12}\left(L,G,S,D\right)\\ m_{21}\left(L,G,S,D\right)& m_{22}\left(L,G,S,D\right)\end{array}\right)=\left(\begin{array}{cc}-1& 0\\ 0& -1\end{array}\right),\end{array}$

Indeed, the transfer matrix MZ is formed by a pair of the same symmetric submatrix [Reference Maksimov, Tyurin and Fedotov14] that MZ = MZ_sub ^∗MZ_sub. Each submatrix describes the transmission of proton in one half of Zumbro lens. The submatrix is written in the following form:

(4) $\begin{array}{}M{Z}_{\text{sub}}=\left(\begin{array}{cc}{q}_{11}\left(L,G,S,D\right)& q_{12}\left(L,G,S,D\right)\\ q_{21}\left(L,G,S,D\right)& q_{22}\left(L,G,S,D\right)\end{array}\right).\end{array}$

The constraints of inverting identity can be given as follows:

(5) $\begin{array}{}q11\left(L,G,S,D\right)+q22\left(L,G,S,D\right)=0.\end{array}$

The FOV is another important constraint. As shown in Figure 1, the proton in FOV should not hit the aperture of magnetic iron or beam pipe when it travels all the way through the lens system. To make it clear for optimization, the Zumbro lens is divided into two fictitious segments arbitrarily. The transfer of an initial proton in the upstream segment can be expressed in the following form:

(6) $\begin{array}{}\left(\begin{array}{c}X\\ X\prime \end{array}\right)={\left(\begin{array}{cc}U{m}_{11}& U{m}_{12}\\ U{m}_{21}& U{m}_{22}\end{array}\right)}_{\text{Upstream}}{\left(\begin{array}{c}{X}_0\\ w\ast X_0+\varphi \end{array}\right)}_{FOV},\end{array}$

where ${\left(\begin{array}{cc}U{m}_{11}& U{m}_{12}\\ U{m}_{21}& U{m}_{22}\end{array}\right)}_{\text{Upstream}}$ is the transfer matrix of upstream segment. At the end of the upstream segment, X and X′ are the coordinate and trajectory angle of the proton, respectively. The coordinate X at any arbitrary divisional plane can be calculated by lens parameters and initial coordinate X ₀. The calculated value of X should be smaller than the inner radius of the quadrupole or beam pipe. As shown in Figure 1, protons located at the edge of FOV have a trace whose coordinate X may reach the maximum. Therefore, the constraint that protons are transmitted without bombarding magnetic iron or beam pipe can be expressed as formula (7). Here, the minimum value of R_quadrupole and P_pipe is set to be 12 cm.

(7) $\begin{array}{}\left(X\left(L,G,S,D,X_0\right)\right)<\min \left(R_{\text{quadrupole}},R_{\text{pipe}}\right)=12cm\left(X_0=R_{FOV}\right).\end{array}$

2.2.3. Optimization Method: Genetic Algorithm

In the aspect of optimization, two constraints for four parameters mean a large space of allowed solutions. The traditional design method of magnetic is based on thin-lens approximation. It takes the size of FOV as a posterior limitation on the parameters. However, this approximation may not provide the best solution.

All the constraints should be taken as object functions to guarantee a reasonable result. The design is expected to be a multiobjective optimization. The genetic algorithm (GA) is recognized as the best globally well-adapted optimization algorithm. It has been successfully employed to deal with the inverse problem in our previous work [Reference Zhang, Hu and Jia19]. Figure 2 shows the flowchart of the genetic algorithm.

Figure 2: The flowchart of optimization with GA.

As shown in Figure 2, four parameters L, G, S, and D of Zumbro lens are represented by a randomly initialized binary vector chromosome. The domains of these variables are set as G ≤ 10 Tesla/m, L∈(0.5/m, 8/m), S∈(3/m, 16/m), and D∈(1.5/m, 8/m). In each generation (evolution iteration procedure), a population with 70 individual chromosomes is used. Then, population evolution processed by three major operators (selection operator, crossover operator, and mutation operator). In the first iteration, the evolution step is skipped because fitness values are not calculated yet.

Fitness function OF1 + OF2 is employed to evaluate the performance (fitness) of each chromosome. The first one, OF1, gives the motivation to minimize the ${\text{m}}_{12}^{\prime }$ . The second one, OF2, functions as constraints. The constant weights C1 and C2 are carefully set to guarantee the effectiveness of the constraints in the optimization.

The fitness value of each individual is then feedback to the evolution step. In this step, the Elitist [Reference Zbigniew18] model used in a selection operator is to make sure the survival of the best individual. The crossover operator performs the inheritance and multiplication of chromosomes. The mutation operator offers the new chromosomes, which may improve the evolution. All in all, with a higher fitness value, the individual has more chance to keep its characteristic gene to next generation. Then, the population evolves from generation to generation until the best fitness is no more improved. A parameter set approximating the best is obtained consequently.

2.2.4. Zumbro Lens for 20–50 GeV/c Protons

The Zumbro systems suitable for 20 to 50 GeV/c protons are optimized. The results are listed in Table 1. The best parameter set practically contains the maximum intensity of magnetic field (G) and the minimum value of D and S. It may be caused by the limitation of the quadrupole inner radius. As shown in Figure 1, in each quadrupole section, the trace of protons should be inside the quadrupole inner radius. The improvement on magnetic intensity G will decrease the offset by which proton trace is turned around. Therefore, the quadrupole with a limited radius desires a high magnetic intensity to avoid the bombarding of protons.

Table 1: Optimized parameters.

A tungsten slab with an areal density of 200 g/cm² is taken as the object. The corresponding chromatic aberrations for difference cases are calculated and given in Table 1. The fourth column offers the calculated full-width at half maximum (FWHM) of the chromatic aberration induced blur. As a comparison, Geant4 is employed to simulate FWHM which can be found in the fifth column.

2.3. Diffuser and Matching Lens

As shown at far left of Figure 1, the upstream diffuser works together with the matching lens. The diffuser provides a desired position-angle correlation beam for matching lens. The basic structure of matching lens consists of three quadrupoles combining with four drifts. The matching lens has two functions. The first one is also to provide a proton beam with desired coordinate-angle correlation for Zumbro lens. The second is to expand the beam size to a fully illuminate object. There are 9 parameters for matching lens: D _2,3,4,, T _1,2,3, and G _1,2,3. The transfer of a proton in matching lens of X-plane has the form:

(8) $\begin{array}{}\left(\begin{array}{c}A\cdot X_0\\ w\cdot A\cdot X_0\end{array}\right)={\left(\begin{array}{cc}M{m}_{11}& M{m}_{12}\\ M{m}_{21}& M{m}_{22}\end{array}\right)}_{\text{Matching}}\left(\begin{array}{c}{X}_0\\ w_0\cdot X_0+{\varphi }_{\text{diffuser}}\end{array}\right),\end{array}$

where A is the amplification factor of matching lens; w is the required correlation at the object plane; w ₀ is the correlation when the proton is at the entrance of the matching section; and ϕ _diffuser is average scattering angle of proton in diffuser. The value of w ₀ not only involves in the matching lens optimization but also closely relates to the thickness of the diffuser.

Compared to the distance from matching lens to the diffuser, the beam size on diffuser is small which can be approximated as a point. Therefore, w ₀ equals to the reciprocal of the draft from a diffuser to matching lens (D ₁, as shown in Figure 1). w ₀ is expected to have a positive correlation with the thickness of the diffuser:

(9) $\begin{array}{}{w}_0\propto T_{\text{diffuser}}, w_0=\frac{1}{D_1},\end{array}$

To decrease the chromatic aberration caused by the diffuser, a thinner diffuser should be used. However, it results in a large distance between the matching lens and diffuser. Therefore, apart from the correlation and amplification, a reasonable value of w ₀ should be carefully selected in the optimization about matching lens.

With a thin diffuser, ϕ _diffuser is close to zero, and formula (8) can be expressed as follows:

(10) $\begin{array}{}\left(\begin{array}{c}M{m}_{11}+w_0\cdot M{m}_{12}=-A,\\ M{m}_{21}+w_0\cdot M{m}_{22}=-w\cdot A.\end{array}\right)\end{array}$

The equation for y plane is similar, but the amplification factor is A. Formula (10) acts as the constraints in the optimization for the matching lens and diffuser.

The value of w is expected as 0.205 for the optimized Zumbro lens. The goal of the matching lens design is to obtain a higher amplification factor. To ensure a thin diffuser and an available drift length, the domain of w ₀ is set to be from 0.1 to 0.13. The overall length is limited within 30 meters. GA is employed to obtain a matching lens pursuing a 10 times amplification. All the optimized parameters are listed in Table 2.

Table 2: Parameters of an optimized diffuser and matching lens.

Two issues should be noted for the data in Table 2. As mentioned for equation (9) that $w_0=1/D_1$ , it is made based on point source assumption. Therefore, 1/w ₀ is just an initial value of D ₁. Geant4 is employed to obtain the exact value of D ₁ by repetition simulations. The secondary issue is relevant to the thickness of the diffuser. The designed amplification factor is 10 and the radius of FOV equals to 65 mm. Therefore, the beam scale at the entrance plane of matching lens is expected to equal to or larger than 6.5 mm in radius. This beam scale is produced by the diffuser.

The object is often located at the center of FOV. To use the proton beam efficiently, the proton intensity at the center of illuminating beam should be higher than that at the edge of FOV. Here, proton intensity at the FOV edge is assumed to be 20% of the maximum value. Then, angular deviation θ ₀ induced by a diffuser can be calculated by the following:

(11) $\begin{array}{}{\theta }_0=\sqrt{\frac{{\theta }^2}{-2\cdot \ln \left(0.2\right)}-{\delta }^2},\end{array}$

where δ is the angular deviation of the initial beam. Its value is assumed as 0.05 mrad. Then, the expected θ ₀ is about 0.44 mrad. In a pure material, the angular deviation θ ₀ induced by MCS can be calculated by using the following formula [Reference Yan, Sheng and Huang22]:

(12) $\begin{array}{}{\theta }_0=\frac{13.6MeV}{\beta cp}\sqrt{\frac{T_{\text{diff}}}{X}}\left(1+0.038\frac{T_{\text{diff}}}{X}\right),\end{array}$

where p is the proton momentum, βc is the velocity, T _diff is diffuser thickness, and X is radiation length in material. With θ ₀ = 0.44 mrad, thickness of the tungsten diffuser T _diff is calculated as 1.5 mm. Geant4 simulation shows that the fraction of momentum spreading caused by this diffuser is on order of 0.01%.

Designed diffuser, matching lens, and Zumbro lens were assembled as a whole system and modeled by Geant4. Then, the proton passing through the diffuser and magnetic lens is simulated. The obtained coordinate-angle correlations for both X and Y plane at the object plane are given in Figure 3(a). It coincides with the design. The simulation also shows that over 95% of protons can pass through this whole system. Figure 3(b) presents the beam distribution at the object plane. At the edge of FOV, the proton flux is about 20% to the maximum value which is also consistent with the design.

Figure 3: Proton beam at object plane. (a) Coordinate-angle correlations. (b) beam intensity distribution.

3. Auxiliary Lens for the Zumbro System

As shown in formula (1), the chromatic aberration is determined by the fractional momentum deviation (σ) and the angular deviation (ϕ) from the ideal position-angle correlation line and the value of $m_{12}^{\prime }$ . The Zumbro lens is optimized according to a reference momentum. To reduce the value of σ, the reference momentum can be selected based on the momentum distribution when the proton leaves the object. For example, when a proton has penetrated a 200 g/cm² tungsten slab, its momentum will distribute in a narrow range. Smaller chromatic aberration can be achieved if the reference momentum is adjusted to be at or nearby this range. However, for object with an irregular shape and areal density, adjusting the reference momentum cannot have the same effect on suppressing chromatic aberration as it is used for a slab shape object.

A new lens system named “auxiliary” is first proposed and shown in Figure 4. It has four different quadrupoles and located downstream of Zumbro lens (reference energy E ₀). The auxiliary lens is designed for the protons with E ₁ energy. Transfer matrix of these protons in Zumbro and auxiliary system has the following form:

(13) $\begin{array}{}{M}_{\left(E_1\right)}=M{Z}_{\left(E_{0,}{E}_1\right)}^n\cdot M{A}_{\left(E_1\right)}=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\end{array}$

where n means the number of Zumbro sections and MA and MZ are the transfer matrices of the incident proton in the auxiliary lens and the designed Zumbro lens, respectively. The auxiliary lens has 13 parameters including four magnetic field grads, four magnetic field lengths, and five drifts. The goal of GA optimization is to searching best parameter set. With these parameters, the transfer matrix of protons for both X and Y plane equals to identity matrix, as given in Formula (13).

Theoretically, by the Zumbro lens system as shown in Figure 4, inverting identity imaging without chromatic aberration can be achieved for E ₀ proton at 1Th image plane. At the 2Th image plane, E ₁ proton will be refocused without chromatic aberration by auxiliary lens.

Figure 4: The structure of auxiliary lens and beam trajectory.

Because proton energy varies from place to place on the image due to the areal densities of object and consequence proton energy loss, the E ₀ proton dominant area on 1Th image has the minimum chromatic aberration. On 2Th image plane produced by auxiliary lens, the area corresponding to E ₁ proton is not overlapped with that of E ₀ proton on 1Th image. It is means that auxiliary lens can offer an extra area with minimum chromatic aberration.

Two objects are used to show the performance of auxiliary. The first sample is a tungsten slab with 28 vacancies. To avoid the thickness effect, the slab thickness is set to 1 cm. As shown in Figure 5, this slab consists of two parts, the top and the bottom half. Densities of these two parts are designed artificially as 200 g/cm³ (red) and 170 g/cm³(blue), respectively. Both parts have the same geometry and vacancy pattern. In each part, two sets of 7 vacancies can be found on left and right. The vacancy widths are 0.5 mm and 1 mm, respectively. The thickness of the 7 vacancies varies from 1 mm to 10 mm by 1.5 mm intervals.

Figure 5: The geometry of slab object. (a) Beam eye view. (b) Oblique view.

20 GeV/c proton beam is employed. The reference momentum for Zumbro lens is designed as 19.78 GeV/c according to the simulated energy loss of this proton in the slab. At the 1Th image plane (see Figure 4), the obtained image provided by the Zumbro lens is shown in Figure 6(a). The auxiliary lens is designed for 19.76 GeV/c protons. The corresponding image captured at 2Th image plane can be found in Figure 6(b).

Figure 6: Recorded image of the slab object with the same proton beam. (a) Given by the Zumbro system (19.78 GeV). (b) Provided by auxiliary lens (19.76 GeV).

As shown in Figure 6(a), if a proton passes through the thickest vacancy (10 mm like a hole), its scattered angle and consequent chromatic aberration blur are close to zero. At the top part, when the protons transmit 8.5 mm thick vacancy, as it penetrates the 1.5 mm slab material, the chromatic aberration blur is the largest. Apart from the thickest vacancy, minimum blur is found when vacancy is 2.5 mm (7.5 mm material with 150 g/cm² areal density) in thickness. It is because the momentum of the proton, which has penetrated 7.5 mm thick material, is close to the designed reference value of Zumbro lens. Compared to the top part, at the bottom, 1 mm thick vacancy (leftmost one) is more clearly. It is because the difference of areal density between top and bottom part.

In Figure 6(b), the thinnest vacancy at the leftmost of the top half is more clearly imaged compared to that in Figure 6(a). It is because auxiliary lens is designed for 19.76 GeV/c protons whose momentum is smaller than the reference value 19.78 GeV/c of Zumbro lens. Therefore, the proton which losses more energy in the object is expected to be better imaged by the auxiliary lens.

Both Zumbro and auxiliary lens have good imaging performance for theirs reference protons, respectively. Therefore, a new image can be generated using the information of the two images given by Zumbro lens and auxiliary lens. The obtained chromatic aberration will be evidently smaller than that obtained by only Zumbro lens.

Another object with complicated areal density distribution is employed. The shape of this object is an iron sphere with bubble shape vacancy at the center. The momentum distribution of the emitted protons is more complicated than the slab object. According to the energy loss of 20 GeV/c proton in this object, the reference momentum of the Zumbro and auxiliary lens are designed as 19.75 GeV/c and 19.78 GeV/c, respectively. Geant4 code is employed for simulation and cross-sectional transmission distributions obtained by Zumbro and auxiliary lens are presented in Figure 7.

Figure 7: Cross-sectional distribution of sphere object. −250 MeV means the reference momentum is 250 MeV lower than that of the incident proton (20 GeV/c). (a) Cross-sectional distribution. (b) Residual error.

The actual attenuation length is determined by the removal cross section. Because the removal cross section varies with material components, object thickness, and the acceptance angle of the lens [Reference Morris, Ables and Alrick8], it is hard to obtain by a simple calculation. Geant4 is therefore employed to calculate the ideal transmission at different positions which can be found in Figure 7(a) (labeled as Geant4 without lens).

As a comparison, another Zumbro lens with 20 GeV/c reference is simulated and given in Figure 7(a) as well. Figure 7(b) shows the residual of these distributions compared to the ideal one (Geant4 without lens).

The energy loss of the proton is different from place to place. The reference momentum of the lens should be varying with the position. For example, to the thickest part locating at 10–20 mm, the reference momentum should be set as 19.75 GeV/c (250 MeV lower than20 GeV/c). At the center or other place shielded by less areal density, a higher momentum ranging from 19.77 GeV/c to 19.80 GeV/c acts better as the reference value.

As shown in Figure 7(b), the transmission residual of 19.75 Zumbro lens is the smallest at the position locating at 10–20 mm away from the center. At the rest positions with 0–10 mm and 20–60 mm distance away from the center, the transmission residual provided by the auxiliary lens is the smallest. If two images given by Zumbro and auxiliary lens are processed properly according to the momentum deviation and corresponding area, the residual error of the combined transmission is only few times 0.1%.