Geant4模拟1.6 GeV质子的输运过程

姚志明 宋顾周 黑东炜 马继明 周 鸣 段宝军 宋 岩 韩长材

（西北核技术研究所 强脉冲辐射环境模拟与效应国家重点实验室 西安 710024）

辐射成像技术是通过测量入射到被测物体上的射线束的参数变化来确定物体内部几何结构和材料组分等信息的[1]。用于射线检测的射线源包括X射线、γ射线、中子、电子、α粒子、质子等。其中质子是带电强子，高能量的质子具有穿透能力强、散射角分布与材料组分相关、源的单色性好和探测效率高等特点。美国科学家利用质子作为探针做了大量的理论和实验工作[1−4]。

质子的电离能损率与电子密度近似成正比，可以用于 CT技术[2]；能量损失在射程末端达到极大值，可以用于治疗癌症[3]；在不同材料中散射角不同，可以获得高对比度的边界图像[1]；与原子核碰撞后发生衰减，多次库仑散射角与材料有关，可以用于高能质子照相技术[4]。总之，质子在医学、材料科学、辐射成像等领域有着广阔的应用前景。

然而，受制于质子加速器的昂贵造价，除医学上的应用外，国内的质子照相技术尚处于理论分析和蒙卡模拟的起步阶段。将于2018年建成的中国散裂中子源中的质子加速器可以将质子加速到1.6GeV[5]。本文将通过Geant4蒙卡模拟软件对该能量下质子在材料中的输运过程进行模拟计算，分析个数衰减、能量损失、散射角分布等参数随材料的面密度和组分的变化规律。

1 质子与物质相互作用的理论公式

质子与物质相互作用的物理过程包括电离能损、多次库仑散射、与原子核的非弹性碰撞以及与原子核的弹性碰撞[4]。其中与原子核的弹性碰撞又可按碰撞后原子核处于基态或激发态划分为弹性散射和准弹性散射[6]。对于同一种物理过程不同的文献中给出的经验公式有所差别，以下各选取了一种进行论述。

1.1 电离能损

质子与电子发生库仑相互作用，能量发生衰减，能损率由Bethe-Bloch公式[6]计算：

式 中， 能损率 的单位 是 MeV·g−1·cm2，K=0.307075 MeV·g−1·cm2；z为事件粒子电荷量（质子z=1）；Z为靶原子电荷量；A为靶原子原子质量；b为质子速度除以光速，即b=v/c；me为电子静止质量；γ2=1/(1−b2)；Tmax是单次碰撞中质子能够传递给电子的最大动能；I是靶原子的平均激发能；C/Z是壳修正项；d/2是密度修正量。

1.2 多次库仑散射

质子受到原子核库仑力的作用，发生多次库仑散射(Multiple Coulomb Scattering, MCS)。每次散射后，质子能量不变，运动方向发生小的改变。空间角分布可以用高斯分布来描述[4]：

式中，θ是质子偏离初始方向的角度；θ0是多次库仑散射角的均方根值，可用下式估计：

式中，p是质子动量；β为质子速度与光速的比值；li为面密度；Ri是材料的辐射长度，经验公式为：

式中，Ri单位是g·cm−2；Z为靶原子电荷量；A为靶原子原子质量。

1.3 核反应

高能质子与原子核中的质子和中子发生非弹性碰撞。质子个数以指数形式衰减[4]：

式中，N为透射质子个数；N0为入射质子个数；li为面密度；λi是平均自由程，单位g·cm−2。λi可用下式估计：

式中，A为相对原子质量；NA为阿伏伽德罗常数；A/NA为每个原子的质量，g；σi是核反应截面，可用下式估计：

式中，ri=r0×A1/3，r0≈1.2 fm。

1.4 弹性散射

质子与原子核发生碰撞，能量不发生损失或者损失非常小的部分，运动方向发生改变。与电离和MCS不同，核反应和弹性散射并不是每个质子都会发生，而是有一定的概率会发生。

如果碰撞时质子能量不发生损失，原子核处于基态，不产生粒子。微分截面可以由下式估计[6]：

式中，p是质子动量；h是普朗克常量；J1为一阶贝塞尔函数；R为原子核黑体有效半径k=p/h。

如果质子能量损失非常小的部分，原子核将处于激发态，或者发射出粒子。微分截面可以由下式估计[6]：

2 Geant4的参数设置

Geant4是CERN开发的一款用于模拟粒子输运过程的软件工具包。与MCNP、EGS等蒙卡模拟软件相比，Geant4具有源代码完全开放的优势，用户可以根据实际需要改进和扩展程序。Geant4用户可以设置的参数主要包括物体的几何结构、粒子与物质反应的物理过程模型、粒子源的信息、粒子记录的信息等[7]。

Geant4的几何结构设置如图1。整个“world”几何体由真空填充，在“world”中心放置一定厚度的足够大平板材料，距离平板5 mm位置放置足够大的平板探测平面，用于记录穿过物体后的粒子信息。粒子与物质反应的物理过程模型选取为G4hMultipleScattering、G4hIonisation、G4LElastic和G4ProtonInelasticProcess[8]，分别对应质子与物质的四种相互作用。能量为1.6 GeV的质子在距离物体一定距离的位置沿z轴方向发射。探测平面记录下透射粒子的能量、动量方向和位置信息。

图1 几何结构模型Fig.1 Geometry model.

3 蒙卡模拟结果与讨论

3.1 电离能损

3.1.1 单一材料中的能量损失

Geant4中物理模型设置为G4hIonisation，计算单个质子在足够厚的钨材料中的输运过程。图2统计出每经过1 cm厚钨材料的能量沉积大小和剩余能量。

可以看出，质子在单一材料中的能量衰减规律是：在射程的前86%，单位长度的能量损失变化不大，而在射程的末端，能量损失率突然增大，称为Bragg峰。利用这一特性，质子被用于肿瘤治疗[3]。

由式(1)，穿过单位面密度的能量损失率与材料的Z/A成正比。穿过一定面密度材料后的能量损失大小可以由下式给出：

式中，ρ是体密度；S(E)是能损率中除去Z/A的部分；ηe是电子密度。

进而可以得到以下表达式：

综上所述，生长抑素辅助治疗急性胰腺炎患者的疗效显著，可有效改善Try-2、CER及AMY水平，同时降低TNF-α、IL-6水平。具有较高的临床推广应用价值。

即测得入射和出射质子能量，就可以重建电子密度，该方法被应用于质子CT技术[2]。

图2 质子在钨中的能量损失Fig.2 Proton energy loss in tungsten.

3.1.2 不同材料中的能量损失

对于不同面密度的铍、铜和钨材料，模拟计算了1.6 GeV质子透射物体后的剩余能量，得到质子的能量损失。三种材料的能量损失随面密度的变化曲线如图3所示。

图3 能量损失随面密度变化曲线Fig.3 Energy loss vs. thickness.

由图 3，能量损失与面密度近似成正比关系，能量损失的大小可以反映已知材料的面密度信息；相同面密度下，不同材料的能量损失大小也不同，因而测得能量损失也可以用于区别不同的材料。

3.2 多次库仑散射

3.2.1 单一材料中的散射分布

Geant4物理模型设置为G4hMultipleScattering，选取面密度为5 g·cm−2、10 g·cm−2和15 g·cm−2的钨材料，分别入射105个质子，记录下每个质子偏离入射方向的角度，统计出角度分布情况，结果如图4、5所示。

图4 多次库仑散射角分布图(dN/dW)Fig.4 MCS angular distribution (dN/dW).

图5 多次库仑散射角分布图(dN/dθ)Fig.5 MCS angular distribution (dN/dθ).

这是高斯分布在空间立体角内的积分，是高斯分布和正弦函数的乘积，纵坐标中的个数除以角度的正弦值，就得到了高斯分布。

式(3)表明散射角的均方根值与材料的辐射长度有关。对于不同面密度的物体，统计出了铍、铜和钨的多次库仑散射角的均方根值，如图6所示。

图6 θ0随面密度变化曲线Fig.6 θ0 vs. thickness.

可以看出，对于单一材料，θ0随面密度而增大，测得θ0的大小就反映了面密度的大小。当已知材料面密度时，不同材料的θ0的值不同，测得θ0的大小就反映了材料的组分信息。然而引起散射的因素还包括与原子核的弹性碰撞以及与原子核非弹性碰撞产生的次级质子，θ0能否精确测量需要更多的探索。

3.3 非弹性碰撞

Geant4物理模型设置为G4ProtonInelasticProc，对于不同面密度的铍、铜和钨，统计出了能量未发生改变，即未发生非弹性碰撞的透射质子个数，图7画出了透射率随面密度的变化曲线。

图7 透射率随面密度变化曲线Fig.7 Transmission vs. thickness.

由图 7，对于单一材料，透射率随面密度的增加呈指数衰减，与式(5)符合得较好。测得透射率大小就反映了面密度的大小。当已知材料面密度时，不同材料的透射率的值不同，测得透射率的大小就反映了材料的组分信息。

非弹性碰撞还会产生次级质子，为了分析次级质子的特性，选取面密度为5 g·cm−2、10 g·cm−2和15 g·cm−2的钨材料进行计算，物理模型设置为G4ProtonInelasticProcess，每次入射质子个数105。模拟计算中发现有能量很低的质子产生，而这部分质子由于能量损失将不能透射物体。因而程序中同时加入G4ProtonInelasticProcess和G4hIonisation物理过程，统计出了透射质子总个数，减去单独用G4ProtonInelasticProcess过程计算时能量未发生改变的直穿质子的个数，就得到了透射的次级质子个数。对于次级质子，图8画出了散射角分布图，次级质子个数较少，bin的划分较大，为50 mrad。表1给出了产生的次级质子个数和多次库伦散射角均方根值大小。

图8 次级质子的散射角分布Fig.8 Scattering angular distribution of secondary protons.

三种面密度下散射角的大小要比MCS的散射角大很多，因而调整了横坐标的取值区间，分布函数都与高斯函数类似，具有中间多，周围少的特点。此外，三种面密度下的分布没有明显的差别。

3.4 弹性散射和准弹性散射

Geant4中物理模型设置为 G4LElastic，G4LElastic模型包括了弹性散射和非弹性散射过程。选取面密度为5 g·cm−2、10 g·cm−2和15 g·cm−2的钨材料，分别入射 105个质子，统计出发生弹性散射的比例(表1)。散射质子的角度分布见图9。

表1 5 g·cm−2、10 g·cm−2和15 g·cm−2钨中的散射角均方根值Table 1 RMS of scattering angle in 5 g·cm−2, 10 g·cm−2, 15 g·cm−2 tungsten.

图9 弹性散射的质子角分布Fig.9 Elastic scattering angular distribution of protons.

弹性碰撞的散射角也与高斯函数类似，具有中间多，周围少的特点，且三种面密度下的散射角分布没有明显的差别。

3.5 讨论

质子穿过物体过程中，能量损失、个数衰减、运动方向发生偏转。式(11)表明，通过测量能量损失可以给出电子的密度；在射程的前86%，能量损失与面密度近似成正比，且相同面密度的不同材料的能量损失大小不同。图7表明，不发生非弹性碰撞的质子个数随面密度的增加呈指数衰减，且相同面密度的不同材料的个数损失不同。由表 1，引起质子运动方向发生偏转的因素包括多次库仑散射、非弹性碰撞产生的次级质子以及与原子核的弹性散射。随着面密度的增加，产生的次级质子和发生弹性散射的质子个数在逐渐增多；多次库仑散射角逐渐增大，次级质子的散射角逐渐减小，弹性散射角没有明显的变化规律。图6表明，对于不同的材料，散射角的均方根值也不同。

4 结语

对于某种已知材料，能量损失的大小、个数损失的多少以及散射角的大小都可以反映材料的面密度。对于相同面密度的物体，不同材料的能量损失、个数损失和散射角的大小也不同，可以通过测量这些物理量来反映材料的组分信息。

1 李家伟. 无损检测手册[M]. 北京: 机械工业出版社, 2002

LI Jiawei. Non-destructive inspection handbook[M]. Beijing: China Machine Press, 2002

2 Schulte R W, Bashkirov V, Loss Klock M C, et al. Density resolution of proton computed tomography[J]. Medical Physics, 2005, 24(4): 1035−1046

3 樊明武. 用于医学诊断和治疗的质子回旋加速器[J].中国工程科学, 2000, 12(2): 9−15

FAN Mingwu. Medical cyclotron used for diagnostic or therapy[J]. Engineering Science, 2000, 12(2): 9−15

4 Morris C L, Ables E, Alrick K R, et al. Flash radiography with 24 GeV/c protons[J]. Journal of Applied Physics, 2011, 109(10): 104905_1−104905_10

5 陈延伟. 中国散裂中子源(CSNS)[J]. 中国科学院院刊, 2011, 26(6): 726−729

CHEN Yanwei. China spallation neutron source[J]. Bulletin of Chinese Academy of Sciences, 2011, 26(6): 726−729

6 刘进, 章林文, 刘军, 等. 快速高能质子照相程序QMCPrad的研制[J]. 强激光与粒子束, 2012, 24(12): 2959−2964

LIU Jin, Zhang Linwen, LIU Jun, et al. Development of code QMCPrad for fast high-energy proton radiography[J]. High Power Laser and Particle Beams, 2012, 24(12): 2959−2964

7 Geant4 User’s Guide For Application Developers [EB/OL]. http://geant4.web.cern.ch/geant4/

8 Physics Reference Manual[EB/OL]. http://geant4.web. cern.ch/geant4/