Geant4模擬1.6 GeV質(zhì)子的輸運(yùn)過程

姚志明 宋顧周 黑東煒 馬繼明 周 鳴 段寶軍 宋 巖 韓長材

（西北核技術(shù)研究所 強(qiáng)脈沖輻射環(huán)境模擬與效應(yīng)國家重點(diǎn)實驗室 西安 710024）

輻射成像技術(shù)是通過測量入射到被測物體上的射線束的參數(shù)變化來確定物體內(nèi)部幾何結(jié)構(gòu)和材料組分等信息的[1]。用于射線檢測的射線源包括X射線、γ射線、中子、電子、α粒子、質(zhì)子等。其中質(zhì)子是帶電強(qiáng)子，高能量的質(zhì)子具有穿透能力強(qiáng)、散射角分布與材料組分相關(guān)、源的單色性好和探測效率高等特點(diǎn)。美國科學(xué)家利用質(zhì)子作為探針做了大量的理論和實驗工作[1?4]。

質(zhì)子的電離能損率與電子密度近似成正比，可以用于 CT技術(shù)[2]；能量損失在射程末端達(dá)到極大值，可以用于治療癌癥[3]；在不同材料中散射角不同，可以獲得高對比度的邊界圖像[1]；與原子核碰撞后發(fā)生衰減，多次庫侖散射角與材料有關(guān)，可以用于高能質(zhì)子照相技術(shù)[4]?？傊|(zhì)子在醫(yī)學(xué)、材料科學(xué)、輻射成像等領(lǐng)域有著廣闊的應(yīng)用前景。

然而，受制于質(zhì)子加速器的昂貴造價，除醫(yī)學(xué)上的應(yīng)用外，國內(nèi)的質(zhì)子照相技術(shù)尚處于理論分析和蒙卡模擬的起步階段。將于2018年建成的中國散裂中子源中的質(zhì)子加速器可以將質(zhì)子加速到1.6GeV[5]。本文將通過Geant4蒙卡模擬軟件對該能量下質(zhì)子在材料中的輸運(yùn)過程進(jìn)行模擬計算，分析個數(shù)衰減、能量損失、散射角分布等參數(shù)隨材料的面密度和組分的變化規(guī)律。

1 質(zhì)子與物質(zhì)相互作用的理論公式

質(zhì)子與物質(zhì)相互作用的物理過程包括電離能損、多次庫侖散射、與原子核的非彈性碰撞以及與原子核的彈性碰撞[4]。其中與原子核的彈性碰撞又可按碰撞后原子核處于基態(tài)或激發(fā)態(tài)劃分為彈性散射和準(zhǔn)彈性散射[6]。對于同一種物理過程不同的文獻(xiàn)中給出的經(jīng)驗公式有所差別，以下各選取了一種進(jìn)行論述。

1.1 電離能損

質(zhì)子與電子發(fā)生庫侖相互作用，能量發(fā)生衰減，能損率由Bethe-Bloch公式[6]計算：

式 中， 能損率 的單位 是 MeV·g?1·cm2，K=0.307075 MeV·g?1·cm2；z為事件粒子電荷量（質(zhì)子z=1）；Z為靶原子電荷量；A為靶原子原子質(zhì)量；b為質(zhì)子速度除以光速，即b=v/c；me為電子靜止質(zhì)量；γ2=1/(1?b2)；Tmax是單次碰撞中質(zhì)子能夠傳遞給電子的最大動能；I是靶原子的平均激發(fā)能；C/Z是殼修正項；d/2是密度修正量。

1.2 多次庫侖散射

質(zhì)子受到原子核庫侖力的作用，發(fā)生多次庫侖散射(Multiple Coulomb Scattering, MCS)。每次散射后，質(zhì)子能量不變，運(yùn)動方向發(fā)生小的改變?？臻g角分布可以用高斯分布來描述[4]：

式中，θ是質(zhì)子偏離初始方向的角度；θ0是多次庫侖散射角的均方根值，可用下式估計：

式中，p是質(zhì)子動量；β為質(zhì)子速度與光速的比值；li為面密度；Ri是材料的輻射長度，經(jīng)驗公式為：

式中，Ri單位是g·cm?2；Z為靶原子電荷量；A為靶原子原子質(zhì)量。

1.3 核反應(yīng)

高能質(zhì)子與原子核中的質(zhì)子和中子發(fā)生非彈性碰撞。質(zhì)子個數(shù)以指數(shù)形式衰減[4]：

式中，N為透射質(zhì)子個數(shù)；N0為入射質(zhì)子個數(shù)；li為面密度；λi是平均自由程，單位g·cm?2。λi可用下式估計：

式中，A為相對原子質(zhì)量；NA為阿伏伽德羅常數(shù)；A/NA為每個原子的質(zhì)量，g；σi是核反應(yīng)截面，可用下式估計：

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

1.4 彈性散射

質(zhì)子與原子核發(fā)生碰撞，能量不發(fā)生損失或者損失非常小的部分，運(yùn)動方向發(fā)生改變。與電離和MCS不同，核反應(yīng)和彈性散射并不是每個質(zhì)子都會發(fā)生，而是有一定的概率會發(fā)生。

如果碰撞時質(zhì)子能量不發(fā)生損失，原子核處于基態(tài)，不產(chǎn)生粒子。微分截面可以由下式估計[6]：

式中，p是質(zhì)子動量；h是普朗克常量；J1為一階貝塞爾函數(shù)；R為原子核黑體有效半徑k=p/h。

如果質(zhì)子能量損失非常小的部分，原子核將處于激發(fā)態(tài)，或者發(fā)射出粒子。微分截面可以由下式估計[6]：

2 Geant4的參數(shù)設(shè)置

Geant4是CERN開發(fā)的一款用于模擬粒子輸運(yùn)過程的軟件工具包。與MCNP、EGS等蒙卡模擬軟件相比，Geant4具有源代碼完全開放的優(yōu)勢，用戶可以根據(jù)實際需要改進(jìn)和擴(kuò)展程序。Geant4用戶可以設(shè)置的參數(shù)主要包括物體的幾何結(jié)構(gòu)、粒子與物質(zhì)反應(yīng)的物理過程模型、粒子源的信息、粒子記錄的信息等[7]。

Geant4的幾何結(jié)構(gòu)設(shè)置如圖1。整個“world”幾何體由真空填充，在“world”中心放置一定厚度的足夠大平板材料，距離平板5 mm位置放置足夠大的平板探測平面，用于記錄穿過物體后的粒子信息。粒子與物質(zhì)反應(yīng)的物理過程模型選取為G4hMultipleScattering、G4hIonisation、G4LElastic和G4ProtonInelasticProcess[8]，分別對應(yīng)質(zhì)子與物質(zhì)的四種相互作用。能量為1.6 GeV的質(zhì)子在距離物體一定距離的位置沿z軸方向發(fā)射。探測平面記錄下透射粒子的能量、動量方向和位置信息。

圖1 幾何結(jié)構(gòu)模型Fig.1 Geometry model.

3 蒙卡模擬結(jié)果與討論

3.1 電離能損

3.1.1 單一材料中的能量損失

Geant4中物理模型設(shè)置為G4hIonisation，計算單個質(zhì)子在足夠厚的鎢材料中的輸運(yùn)過程。圖2統(tǒng)計出每經(jīng)過1 cm厚鎢材料的能量沉積大小和剩余能量。

可以看出，質(zhì)子在單一材料中的能量衰減規(guī)律是：在射程的前86%，單位長度的能量損失變化不大，而在射程的末端，能量損失率突然增大，稱為Bragg峰。利用這一特性，質(zhì)子被用于腫瘤治療[3]。

由式(1)，穿過單位面密度的能量損失率與材料的Z/A成正比。穿過一定面密度材料后的能量損失大小可以由下式給出：

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

進(jìn)而可以得到以下表達(dá)式：

綜上所述，生長抑素輔助治療急性胰腺炎患者的療效顯著，可有效改善Try-2、CER及AMY水平，同時降低TNF-α、IL-6水平。具有較高的臨床推廣應(yīng)用價值。

即測得入射和出射質(zhì)子能量，就可以重建電子密度，該方法被應(yīng)用于質(zhì)子CT技術(shù)[2]。

圖2 質(zhì)子在鎢中的能量損失Fig.2 Proton energy loss in tungsten.

3.1.2 不同材料中的能量損失

對于不同面密度的鈹、銅和鎢材料，模擬計算了1.6 GeV質(zhì)子透射物體后的剩余能量，得到質(zhì)子的能量損失。三種材料的能量損失隨面密度的變化曲線如圖3所示。

圖3 能量損失隨面密度變化曲線Fig.3 Energy loss vs. thickness.

由圖 3，能量損失與面密度近似成正比關(guān)系，能量損失的大小可以反映已知材料的面密度信息；相同面密度下，不同材料的能量損失大小也不同，因而測得能量損失也可以用于區(qū)別不同的材料。

3.2 多次庫侖散射

3.2.1 單一材料中的散射分布

Geant4物理模型設(shè)置為G4hMultipleScattering，選取面密度為5 g·cm?2、10 g·cm?2和15 g·cm?2的鎢材料，分別入射105個質(zhì)子，記錄下每個質(zhì)子偏離入射方向的角度，統(tǒng)計出角度分布情況，結(jié)果如圖4、5所示。

圖4 多次庫侖散射角分布圖(dN/dW)Fig.4 MCS angular distribution (dN/dW).

圖5 多次庫侖散射角分布圖(dN/dθ)Fig.5 MCS angular distribution (dN/dθ).

這是高斯分布在空間立體角內(nèi)的積分，是高斯分布和正弦函數(shù)的乘積，縱坐標(biāo)中的個數(shù)除以角度的正弦值，就得到了高斯分布。

式(3)表明散射角的均方根值與材料的輻射長度有關(guān)。對于不同面密度的物體，統(tǒng)計出了鈹、銅和鎢的多次庫侖散射角的均方根值，如圖6所示。

圖6 θ0隨面密度變化曲線Fig.6 θ0 vs. thickness.

可以看出，對于單一材料，θ0隨面密度而增大，測得θ0的大小就反映了面密度的大小。當(dāng)已知材料面密度時，不同材料的θ0的值不同，測得θ0的大小就反映了材料的組分信息。然而引起散射的因素還包括與原子核的彈性碰撞以及與原子核非彈性碰撞產(chǎn)生的次級質(zhì)子，θ0能否精確測量需要更多的探索。

3.3 非彈性碰撞

Geant4物理模型設(shè)置為G4ProtonInelasticProc，對于不同面密度的鈹、銅和鎢，統(tǒng)計出了能量未發(fā)生改變，即未發(fā)生非彈性碰撞的透射質(zhì)子個數(shù)，圖7畫出了透射率隨面密度的變化曲線。

圖7 透射率隨面密度變化曲線Fig.7 Transmission vs. thickness.

由圖 7，對于單一材料，透射率隨面密度的增加呈指數(shù)衰減，與式(5)符合得較好。測得透射率大小就反映了面密度的大小。當(dāng)已知材料面密度時，不同材料的透射率的值不同，測得透射率的大小就反映了材料的組分信息。

非彈性碰撞還會產(chǎn)生次級質(zhì)子，為了分析次級質(zhì)子的特性，選取面密度為5 g·cm?2、10 g·cm?2和15 g·cm?2的鎢材料進(jìn)行計算，物理模型設(shè)置為G4ProtonInelasticProcess，每次入射質(zhì)子個數(shù)105。模擬計算中發(fā)現(xiàn)有能量很低的質(zhì)子產(chǎn)生，而這部分質(zhì)子由于能量損失將不能透射物體。因而程序中同時加入G4ProtonInelasticProcess和G4hIonisation物理過程，統(tǒng)計出了透射質(zhì)子總個數(shù)，減去單獨(dú)用G4ProtonInelasticProcess過程計算時能量未發(fā)生改變的直穿質(zhì)子的個數(shù)，就得到了透射的次級質(zhì)子個數(shù)。對于次級質(zhì)子，圖8畫出了散射角分布圖，次級質(zhì)子個數(shù)較少，bin的劃分較大，為50 mrad。表1給出了產(chǎn)生的次級質(zhì)子個數(shù)和多次庫倫散射角均方根值大小。

圖8 次級質(zhì)子的散射角分布Fig.8 Scattering angular distribution of secondary protons.

三種面密度下散射角的大小要比MCS的散射角大很多，因而調(diào)整了橫坐標(biāo)的取值區(qū)間，分布函數(shù)都與高斯函數(shù)類似，具有中間多，周圍少的特點(diǎn)。此外，三種面密度下的分布沒有明顯的差別。

3.4 彈性散射和準(zhǔn)彈性散射

Geant4中物理模型設(shè)置為 G4LElastic，G4LElastic模型包括了彈性散射和非彈性散射過程。選取面密度為5 g·cm?2、10 g·cm?2和15 g·cm?2的鎢材料，分別入射 105個質(zhì)子，統(tǒng)計出發(fā)生彈性散射的比例(表1)。散射質(zhì)子的角度分布見圖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 彈性散射的質(zhì)子角分布Fig.9 Elastic scattering angular distribution of protons.

彈性碰撞的散射角也與高斯函數(shù)類似，具有中間多，周圍少的特點(diǎn)，且三種面密度下的散射角分布沒有明顯的差別。

3.5 討論

質(zhì)子穿過物體過程中，能量損失、個數(shù)衰減、運(yùn)動方向發(fā)生偏轉(zhuǎn)。式(11)表明，通過測量能量損失可以給出電子的密度；在射程的前86%，能量損失與面密度近似成正比，且相同面密度的不同材料的能量損失大小不同。圖7表明，不發(fā)生非彈性碰撞的質(zhì)子個數(shù)隨面密度的增加呈指數(shù)衰減，且相同面密度的不同材料的個數(shù)損失不同。由表 1，引起質(zhì)子運(yùn)動方向發(fā)生偏轉(zhuǎn)的因素包括多次庫侖散射、非彈性碰撞產(chǎn)生的次級質(zhì)子以及與原子核的彈性散射。隨著面密度的增加，產(chǎn)生的次級質(zhì)子和發(fā)生彈性散射的質(zhì)子個數(shù)在逐漸增多；多次庫侖散射角逐漸增大，次級質(zhì)子的散射角逐漸減小，彈性散射角沒有明顯的變化規(guī)律。圖6表明，對于不同的材料，散射角的均方根值也不同。

4 結(jié)語

對于某種已知材料，能量損失的大小、個數(shù)損失的多少以及散射角的大小都可以反映材料的面密度。對于相同面密度的物體，不同材料的能量損失、個數(shù)損失和散射角的大小也不同，可以通過測量這些物理量來反映材料的組分信息。

1 李家偉. 無損檢測手冊[M]. 北京: 機(jī)械工業(yè)出版社, 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 樊明武. 用于醫(yī)學(xué)診斷和治療的質(zhì)子回旋加速器[J].中國工程科學(xué), 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]. 中國科學(xué)院院刊, 2011, 26(6): 726?729

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

6 劉進(jìn), 章林文, 劉軍, 等. 快速高能質(zhì)子照相程序QMCPrad的研制[J]. 強(qiáng)激光與粒子束, 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/