用LabVIEW实现NaI单晶的γ能谱模拟

本项目在LabVIEW的环境下,利用Monte Carlo方法,实现NaI单晶对γ光子沉积能量的模拟。

说明:

  1. 本网页初衷是学习交流,故所有程序开源,如有兴趣可阅读并修改开。
  2. 以下程序均在Labview2012上编写,更低版本无法打。
  3. 下载VI时请按照所列顺序依次下载,因为之后的程序可能会调用之前的程序。
  4. 主要参考文献:许淑艳. 蒙特卡罗方法在实验核物理中的应用[M]. 北京:北京原子能出版社, 2006. 下面未注明出处的分布公式及该分布的抽样方法均出自此书。
  • γ光子由放射源产生,其能量取决于具体放射源,而在空间中的位置和运动方向满足一定的分布关系,需要由抽样生成,并且只有一部分光子能进入晶体并对沉积能量有所贡献。[详见:光子的位置和角度抽样]
  • 进入晶体的光子可能与不同的核碰撞,发生不同的反应,具体和哪种核碰撞,发生什么反应是一个概率问题,用反应截面来描述。应计算晶体中两种原子各种反应类型的反应截面并抽样决定碰撞核的种类和反应类型。[详见:反应截面计算及碰撞核和类型的确定]
  • 发生的反应类型可能有光电效应、Compton散射和对生成。这里由于时间关系,将目标限于能量较低的光子,故只考虑光电效应和Compton散射。[详见:光电效应过程 Compton散射过程]
  • 不论是这两种过程中的哪一种都会产生次级电子,应计算次级电子射程并判断其是否逃离晶体,从而确定出电子在晶体中运动的有效距离。[详见:电子射程计算]
  • 电子在运动过程中会辐射出轫致光子,这些轫致光子的总能量和上述有效距离有关,应确定这些光子的分布和能量。[详见:电子的模拟]
  • 因此一个光子入射,会有包括次级光子和轫致光子在内的许多具有不同能量,不同位置、角度分布的光子,这些光子全都应该被记录下来,重新参与到光子的模拟过程中去。[详见:一个光子的模拟]
  • 上述一个光子的模拟过程作为一次沉积能量被NaI记录,沉积能谱实际是对许多次入射光子的许多个沉积能量所做的一个统计。[详见:沉积能量的计算]

在位置抽样程序中,我们假设放射源是一个点源且位于NaI晶体的正下方,以放射源为球心做一个刚好包含NaI晶体的球,进一步假设γ光子在该球内均匀分布。因此抽样生成位于这个球内的光子的位置,只有当这个光子位于晶体内,才认为这个光子进入了晶体而对能量沉积有贡献,否则重新抽样。

事实上,单从放射源发射光子的角度来讲,如果光子的位置确定了,那么由源指向光子的矢量便确定了这个光子的运动方向。因此如果仅考虑球内均匀分布假设,我们不需要新的程序来确定光子的方向。但考虑到晶体会对光子发生作用,上述程序中生成的位于晶体内的光子被视作已发生多次作用的光子,不妨也假设它的方向也完全随机。

此处考虑的反应是光电效应和Compton散射,它们的反应截面都和能量大小有关。具体公式相当复杂,但一个为人熟知的结论是:若入射光子能量远小于电子静能量(0.511MeV),Compton散射过渡为Thomson散射,其散射截面有简单的表达式sigma_{Thomson}={8/3}{pi}{Z}{r_{0}}^2;在当入射光能量远大于电子静能量是,Compton散射截面也能有简单的表达式sigma_{photo}={Z/{h nu}}。上述程序在Compton散射截面的计算上仅运用了这些近似公式,因此显然对0.511MeV附近γ光子,模拟会出现很大误差。

而对于光电效应,Wikipedia上也给出了以电子静能量为分界的几个反应截面的渐进公式,程序中同样采用的是这些渐进公式。

应用已经计算反应截面的子程序可以判断碰撞核的种类和反应类型。给定入射能量,对每种核而言,光电效应截面和Compton散射截面之和为反应总截面,首先比较反应总截面大小,可以抽样产生碰撞的核的类型。其次比较该种核的光电效应截面和Compton散射截面,可以确定反应类型。

如果发生光电效应,那么光子会把能量交给碰撞的核的K层电子使其获得能量,而K层电子原先的结合能即作为沉积能量被晶体吸收。光电子的散射极角(光电子的运动方向与光子运动方向的夹角)满足如下分布:

alpha={sqrt{1-beta^2}/{3beta}}delim{[}{{{3}({beta}{cos theta})^2-3beta{cos theta}+1-{beta}^2}/{(1-{beta}{cos theta})^2}-{1-2beta}/{(1-beta)^2}}{]}+{(1-sqrt{1-beta^2})(1-2sqrt{1-beta^2})}/{4beta{(1-beta^2)}}delim{[}{{3-beta}/{1-beta}-{beta^2+3-{4beta}{cos theta}}/{(1-{beta}{cos theta})^2}-2ln({1-{beta}{cos theta}}/{1-beta})}{]}

其中beta=lbrace{1-{delim{[}{{{m_{0}}{c^2}}{/}({E_{gamma}-E_{K}+{m_{0}}{c^2}})}{]}}^2}rbrace^{1/2}

Compton散射会产生一个散射光子和一个散射电子,令alpha_{1}alpha_{2}分别为入射前后光子的能量,x=alpha_{1}/alpha_{2},则光子的能量分布满足:

f(x)=delim{lbrace}{matrix{2}{2}{{{1/K(alpha_{1})}delim{[}{({alpha_{1}+1-x}/{{alpha_{1}}{x}})^2+1/x-1/{x^2}+1/{x^3}}{]}}{{when}{1<=x<=2alpha_{1}+1}}{0}{else}}}{} ,其中K(alpha_{1})=ln(2alpha_{1}+1)delim{[}{1-{2(alpha_{1}+1)}/{{alpha_{1}}^2}}{]}+1/2+4/{alpha_{1}}-1/{2(2alpha_{1}+1)^2}

抽样产生光子能量后,根据能量守恒可以得到电子能量。

考虑散射光子和电子的运动方向。光子的散射角的余弦值分布满足:

f(cos theta)={delta}{(1-1/{alpha_{2}}+1/{alpha_{1}}-cos theta)}

抽样产生光子的散射角后,根据动量守恒可以推导出电子运动方向和散射前Compton运动方向夹角满足如下关系:

cos theta={(1+alpha_{1}){tan theta_{gamma}/2}}/{sqrt{1+{(1+alpha_{1})^2}{tan^2 theta_{gamma}/2}}}

考虑电子射程应考虑晶体对电子运动过程的影响,即计算NaI晶体对电子的散射截面,其中应用的公式参考了Rutherford散射实验中散射截面的关系,即与靶材料质量数Z平方成正比,与入射粒子能量平方成反比。

电子自由程rho分布满足指数关系:f(rho)=e^{-rho},对于均匀介质,电子射程L和上述散射截面存在关系:L=rho/{sigma N},其中sigma是NaI对电子的散射截面,N是NaI的数密度。计算电子射程是为了用它来计算之后辐射的轫致光子。因此有必要判断电子经历了这段射程后是否逃离晶体。如果离开,那么有效的射程只应计算位于晶体内的一段。

电子在运动的过程中会辐射轫致光子,这些光子会参与上述模拟过程,不得不进行考虑。但从文献《电子的轫致辐射光子谱计算》中可以看出,这些光子具有相当复杂的能量分布,要细致地考虑它们,一是对计算要求相对较大,二是存储这些轫致光子能量、位置、角度的信息对空间的要求也更大。这里做了一个简化,我们利用文献中给出的公式计算出辐射光子总能量,同时文献也给出了辐射光子能量分布的范围,我们假设这些光子具有相同的能量,随机生成落在上述能量范围内的一个能量大小作为一个轫致光子的能量大小,用总能量除以这个能量所得结果取整得到轫致光子数。

至此,一个电子的过程已经可以被完全分析清楚了。该电子的能量和运动方向取决于是Compton效应还是光电效应,碰撞的核是Na还是I,能量确定了以后,它的射程,以及辐射出的轫致光子数和轫致光子能量也可以由抽样决定。当然为了让这些轫致光子可以投入到进一步的模拟中,我们必须规定这些轫致光子的位置和运动方向。我们假设这些轫致光子是在电子运动过程中均匀产生的,即沿着电子的射程,每运动相同的距离辐射出一个光子;而运动方向和最初的抽样一样,具有任意性。当然我们产生的轫致光子总能量和电子能量是不相同的,两者的差值是电子在晶体运动中交给晶体的能量,即沉积能量,在这个程序中也应该被给出。

通过对电子的模拟,已经将不论什么过程中的电子转化为轫致光子进行分析。接下来我们只需要考虑光子的模拟即可,但并非所有的光子都对沉积能量有所贡献,因为光子运动的自由程分布决定了一部分光子可能未经碰撞直接逃离系统。判断光子是否逃离系统利用的方法是,根据反应截面,抽样生成光子的自由程,再结合光子的初始位置和速度方向,计算经历这段自由程后,光子是否还在晶体内。

和电子模拟类似,光子模拟主要在电子模拟的基础上,对光电效应过程,将一个光子分成若干轫致光子;对Compton散射过程,将一个光子分成一个Compton光子和若干轫致光子。不论是Compton光子还是轫致光子,都需要输出它的能量、位置和角度。另外,初始光子究竟发生Compton散射还是光电效应这一点我们也是需要关注的,因为从反应截面的源程序可以看出,光子发生Compton散射和光电效应的概率以电子静能量为划分。如果一开始光子发生光电效应,那么显然,所有产生的轫致光子不可能发生Compton散射。只有对那些一开始就发生Compton散射的光子,才需要判断散射光子以及轫致光子会不会再发生Compton散射而产生更多不同类型的光子。

我们最后要算沉积能量,但之前的程序中只有光电效应和电子模拟过程给出了沉积能量的具体求法。要求发生Compton散射的光子对沉积能量的贡献,应该考虑一个光子经历了多次碰撞产生了许多光子,这个过程中次级电子在每次模拟时会有一个沉积能量的贡献。而最后程序终止的条件是:生成的所有光子中的每个光子要么逃离晶体、要么能量小到只能发生光电效应。然后再考虑这些只能发生光电效应的电子就可以了。在程序实现的过程中,考虑到如果一次Compton散射后的光子还能发生Compton散射,那么它又将有一个电子产生,进而辐射出若干轫致光子。这些轫致光子的数目、能量、位置都和第一次Compton散射不同。因此这里程序相当于用了一个2维数组来记录这些轫致光子的信息。这个2维数组共有5行,第一行记录轫致光子数,第二行记录轫致光子能量,最后三行记录轫致光子位置。这个2维数组在今后的调用过程中是按列调用的,在程序中引入一个相当于指针的变量,便于以后调用。

我们已经为一个光子的模拟过程做了充足的准备。现在只需调用之前的一系列子VI,便可以产生一次抽样得到的沉积能量。

反复抽样,抽样次数由用户决定。将入射能量分为100份,可以对应实验中的道址。每一次的沉积能量落入这100个道址中的一个,对次数做一统计,便得到沉积能量谱。

  1. 建设网页中,部分显然用到的假设已用粗体标出,但并不排除程序中用到了其他未列出的假设。这些假设不全都合理。
  2. 所用到的公式很多都为近似公式,一方面是时间原因,一方面有些公式确实很难找到。
  3. 未考虑真实测量时的统计涨落,因此从下一节的演示视频可以看出,全能峰是单一峰值,没有展宽。

当然时间仓促,问题远不止以上这些,烦请各位读者将阅读程序时发现的诸多问题在留言区进行指正,造福群众。

  • 结果演示Flash

模拟结果演示

  • 期末报告PPT下载PPT

感谢复旦大学物理教学实验中心俞熹副教授对网站建设和LabVIEW编程的帮助。

感谢复旦大学物理教学实验中心乐永康副教授对NaI单晶γ能谱实验的悉心指导。

感谢复旦大学现代物理研究所王旭飞副教授在截面计算公式代入方面给出的帮助。

留言区

这个项目已经做完了么?感觉好有创意啊! — 沈金辉 2015/01/06 12:12
这两天在建这个网页,作为近代物理期末的大报告。其实也不算做完,因为很多细节考虑很不仔细,所以我准备把子程序都放上来,并且把我写的思路、做了哪些假设写明,以方便如果有人有兴趣的话可以修改。这么快就有人发觉好震惊!其实也不算很有创意啦~当初查gamma能谱实验文献的时候看到有人在别的环境下实现了,然后正好最后选做实验选的是LabVIEW,就想尝试一下用Labview来实现看看~ —吴征东
很不错的尝试。 — 乐永康 2015/01/06 23:58
谢谢老师鼓励~ — 吴征东
再改进一下,可以发文章 — 俞熹 2015/01/09 21:49
  • home/whyx/proj/labview_gamma.txt
  • 最后更改: 2015/01/09 21:50
  • (外部编辑)