本发明属于粒子输运模拟领域，具体地说涉及一种基于舍选抽样的电子多次散射角抽样方法。
背景技术：
蒙特卡罗方法在粒子输运模拟中应用非常广泛，因为其直接追踪粒子的输运过程，
方法直观；所需物理简化少，所以精度高；同时模拟方法不依赖于具体计算量，故通用性强，输出灵活全面。但是不同于中子、光子等中性粒子，电子等带电粒子因为库仑长程作用，在物质中输运时碰撞非常频繁，比如电子能量从mev量级降低到不能再对介质原子激发之前，会发生105—106次碰撞，其中绝大部分都是库仑散射。所以对于电子，不能像光子、中子的蒙卡模拟那样对每次碰撞进行逐次模拟，通常的方法是“压缩历史”方法，将一段距离上的很多次反应“压缩”成一个准反应来模拟其累积效果，即用“多次散射角”近似多次库仑散射对电子运动方向的角度偏转效应。
多次散射角不是一个具体值，而是分布在[0,π]上的一个变量，计算其分布函数的方法即多次散射理论。目前的多次散射角理论，主要有molière方法、goudsmit-saunderson(以下简称gs)理论和lewis理论，molière方法在早期的蒙卡输运程序中使用较多，但因其步长短时误差大，其所用截面在电子能量mev量级时也不准确，所以逐渐被后两种方法替代。
gs理论和lewis理论形式相近，以gs多次散射理论为例，经距离为t的步长后多次散射角余弦ω的分布函数可以写成一个级数展开式
其中λ为散射平均自由程，pl为l阶勒让德多项式，gl称为散射矩，它依赖于单次散射的微分截面。这种方法虽然精度高，但级数展开式收敛慢，特别是在能量较高或步长很短时。以egs5程序为例，该级数展开式最高要展开到15000阶，所以每计算一次该函数的值都非常费时。但目前的算法为了从该分布函数抽样多次散射角，一般需要计算该分布函数在多个角度的值，然后依据累积概率密度函数抽样多次散射角，这样计算量很大，限制了模拟效率，一个尽可能少地计算该函数的值的抽样方法是非常必要的。
技术实现要素：
本发明所要解决的技术问题是提供一种减少重复抽样的参数，提高抽样效率的基于舍选抽样的电子多次散射角抽样方法。
为了解决上述技术问题，本发明采用如下技术方案：一种基于舍选抽样的电子多次散射角抽样方法，包括以下步骤：
(1)构建多次散射角θ的小角近似的分布函数fsa(θ)；
(2)计算多次散射角θ的理论分布函数fms(θ)与fsa(θ)之商f(θ)的最大值fmax；
(3)根据fsa(θ)直接抽样多次散射角θsample；
(4)计算f(θsample)，并抽样[0,1]内的随机数ζ，如果ζ<f(θsample)/fmax，则接受θsample，否则返回步骤(3)。
进一步地，将原始的精确多次散射角分布fms(θ)分解为两部分之积fms(θ)＝fsa(θ)×f(θ)，前一部分fsa(θ)进行直接抽样，另一部分f(θ)作为舍选函数。
进一步地，使用小角近似得到的分布fsa(θ)作为直接抽样所用的分布函数，要求小角近似分布公式有解析形式且比较简单易于直接抽样。
进一步地，小角近似分布公式的选择要使f＝fms/fsa的最大值在确定的有限几个点如θ＝0和θ＝π，这样避免反复寻找舍选函数f的最大值。
本发明还提供另一种基于舍选抽样的电子多次散射角抽样方法，包括以下步骤：
(1)构建多次散射角余弦ω的小角近似的分布函数fsa(ω)；
(2)计算多次散射角余弦ω的理论分布函数fms(ω)与fsa(ω)之商f(ω)的最大值maxfsa/gs；
(3)根据fsa(ω)直接抽样多次散射角ωgauss；
(4)计算f(ωgauss)，并抽样[0,1]内的随机数ζ，如果ζ<f(ωgauss)/maxfsa/gs，则接受ωgauss，否则返回步骤(3)。
进一步地，另一种方法中，fsa(ω)的表达式为
进一步地，另一种方法中，fms(ω)选用gs理论或lewis理论的多次散射角余弦分布函数。
进一步地，另一种方法中，步骤(3)中，当抽样得到的ωgauss<-1时，则重新抽样ωgauss。
本发明的有益效果为：
多次散射角θ的值域是[0,π]，在此范围内θ和其余弦ω＝cosθ有一一对应关系，所以抽样θ和cosθ是等效的，因此上述两种方法是等效的。
本发明相当于将理论分布函数fms分解为一个小角近似的简化分布函数和一个修正函数，即fms＝fsa*f，其中fsa是小角近似得到的多次散射角分布，修正函数f＝fms/fsa，并使用fms与fsa之商，也即修正函数f，作为拒绝函数来进行舍选抽样，从而减少重复抽样的参数，提高抽样效率。
本发明函数fsa(ω)从小角近似得到，其解析形式简单，便于直接抽样，而且只需要计算在抽样角余弦的fms函数值，大大降低了计算量。
本发明在输运模拟中需要计算多次散射角时，首先从小角近似的简化分布函数抽样得到一个初始多次散射角，将修正函数作为拒绝抽样函数判断是否接受该抽样值。小角近似的分布函数的选择以获得尽可能高的抽样效率为目标，这样能减少重复抽样的参数。可以看出，如果抽样效率是η，则本专利所述的方法只需计算1/η次fms(ω)函数的值。
本发明能快速精确地抽样电子在介质中输运的多次散射角，提高电子蒙特卡罗输运模拟的效率，可广泛应用于辐射输运模拟领域、蒙特卡罗粒子输运模拟等领域。
附图说明
图1是本发明实施例2的流程图。
具体实施方式
下面结合实施例对本发明作进一步描述：
实施例1
一种基于舍选抽样的电子多次散射角抽样方法，包括以下步骤：
(1)根据输运步长、粒子能量、输运所在介质和截面数据库数据，构建多次散射角θ的小角近似的分布函数fsa(θ)；
(2)计算多次散射角θ的理论分布函数fms(θ)与fsa(θ)之商f(θ)的最大值fmax；
(3)根据fsa(θ)直接抽样多次散射角θsample；
(4)计算f(θsample)，并抽样[0,1]内的随机数ζ，如果ζ<f(θsample)/fmax，则接受θsample，否则返回步骤(3)。
本实施例中，将原始的精确多次散射角分布fms(θ)分解为两部分之积fms(θ)＝fsa(θ)×f(θ)，前一部分fsa(θ)进行直接抽样，另一部分f(θ)作为舍选函数。
本实施例中，使用小角近似得到的分布fsa(θ)作为直接抽样所用的分布函数，要求小角近似分布公式有解析形式且比较简单易于直接抽样。
本实施例中，小角近似分布公式的选择要使f＝fms/fsa的最大值在确定的有限几个点如θ＝0和θ＝π，这样避免反复寻找舍选函数f的最大值。
本实施例提供的基于舍选抽样的电子多次散射角抽样方法，多次散射角分布可以分成两部分：一个小角近似的简化分布函数和一个修正函数，即fms(θ)＝fsa(θ)*f(θ)，其中fsa(θ)是小角近似得到的多次散射角分布，修正函数f(θ)＝fms(θ)/fsa(θ)。
对于给定的电子能量、输运步长、介质材料组成：
首先利用小角近似得到的分布函数fsa直接抽样得到多次散射角θ；
使用fms与小角分布函数之商作为拒绝函数来进行舍选抽样：计算精确的多次散射分布fms与fsa之商f的最大值fmax，计算f(θ)，并抽样[0,1]内的随机数ζ，如果ζ<f(θ)/fmax则接受该θ，否则重新抽样一个θ值。
fsa函数从小角近似得到，要求其解析形式比较简单，便于直接抽样，而且只需要计算在抽样角θ的fms函数值，大大降低了计算量。该方法使用小角近似得到的分布要满足下列要求：
小角近似分布公式形式要比较简单，以便于从该分布直接抽样；
小角近似分布公式的选择要使f＝fms/fsa的最大值在两个端点，即θ＝0和θ＝π，以避免反复计算f函数的值；
同时小角近似分布公式也要尽可能接近原始分布，特别是在θ＝0附近，以提高抽样效率。
实施例2
多次散射角θ的值域是[0,π]，在此范围内θ和其余弦ω＝cosθ有一一对应关系，所以抽样θ和cosθ是等效的。而多次散射角余弦ω的分布函数形式更简单，所以本实施例将使用ω的分布函数抽样ω。
与实施例1相同，本实施例也是相当于将ω的理论分布函数fms(ω)分解为一个小角近似的简化分布函数和一个修正函数，即fms(ω)＝fsa(ω)*f(ω)，其中fsa(ω)是小角近似得到的多次散射角分布，f(ω)为修正函数，f(ω)＝fms(ω)/fsa(ω)。
fms(ω)可选用gs理论或lewis理论的多次散射角余弦分布函数，现在以gs理论的多次散射理论为例，参见图1对基于舍选抽样的电子多次散射角抽样方法的具体实施方式进行说明如下：
(1)依据电子能量和材料组成从数据库得到散射矩gl值；
(2)为了加快收敛，通常使用至少发生一次碰撞的g-s分布，即步长t上至少发生一次碰撞的多次散射角余弦ω的分布函数，该分布函数的级数展开式为：
(3)一个最简单的小角近似为高斯近似，即多次散射角的余弦值的分布近似为
(4)小角近似下ω的平均值为
(5)计算
(6)抽样[0,1]上的两个随机数ξ1和ξ2，则按高斯分布抽样得到的多次散射角余弦为
(7)计算f(ωgauss)，也即计算fgs(ωgauss)/fsa(ωgauss)，并抽样[0,1]内的随机数ξ3，如果ξ3<f(ωgauss)/maxfsa/gs则接受该多次散射角余弦抽样值，否则回到步骤6。
本发明未详细阐述的部分属于本领域公知技术。
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。