MC计算柱通量的指向概率方法 王瑞宏裴鹿成 (中国原子能科学研究院,北京102413) 摘要本文重点研究MC计算每次碰撞对柱探测器的贡献中对圆柱张角的数值化处理 方法,从源粒子的自然极角与方位角出发直接进行相应的处理,避免了变量替换引起的计算 量增加,变二重积分为单重积分。 1引言 用蒙特卡罗方法对粒子输运进行模拟,几乎都归结为计算通量、其中分计算体通量、面 通量、点通量三种,当体或面较小时,对于体通量或面通量计算都可以很好地转化为点通量 计算问题因此用蒙特卡罗方法计算点通量受到了普遍重视。,2 用蒙特卡罗方法计算点通量时遇到的最大困难是,估计量无界与方差发散,由此而引起 的问题是,计算结果差(统计涨落极大)和无法给了出计算结果的误差。为解决估计量无界 和方差发散的问题,近几年出现了直接计算球或柱的通量的MC方法。34,5,1在MC计算 粒子对有限探测器贡献中统计估计方法是非常有效的降低方差技巧,在粒子的跟踪历史 中,每发生一次散射都计算其对探测器的贡献,于是只需每次根据当时的散射概率分布函数 而无需改变它的权重,如果粒子散射发生在探测器内,粒子历史被终止,在每个历史结束后, 各次散射对探测器的贡献之和即为当前历史对测器的统计概率。123但在计算每次散射对 探测器贡献时都很花时间,这是由于散射点对探测器所张立体角很难确定,按通常方法运算 是很复杂的。 本文从原粒子的自然极角与方位角出发直接进行相应的散值化处理,目的是避免因变 量变换所引起的计算量的增加,更重要的是这种处理方法可以变二重积分为单重积分,使数 值化处理变得非常简单。 在本文第二部分讨论对柱探测器所张立体角的方位角的确定;第三部分讨论散射角的 确定;第四部分讨论统计概率的贡献计箅;第五部分为几点讨论。 2方位角的抽样确定 入射方向 原始问题如右图1所示,其中已知条件是:柱探测器 方向:a,v,to 中心:x0,yo,z 半径:R 长度:2h 入射方向u,v,to 为使计算过程得到简化,以人射方向为z轴,如图2 图1几何示意图
R(xyz P( yz) Pe,y x2y2) 图2几何示意图 0=(u'xo+y0)(1-u2)-l2-(1-v2)12:动o (1) x 同样P1(x1,y,a1),P2(x2,y2,n)以及P1P2所在直线(柱轴)方向余弦u,v,t可以确定 P'(x1,y1),p(x2,y2)坐标确定后由余弦定理 OP12+1OP22-1P'P212 CoS 2|OP|·OP21 方位角∞osθ从(co∞φ,1)中抽样确定,精确地讲,应从于两椭圆相切的切点连线范围内 抽样,同样地在定散射角时,也应在两端做细致考虑,这里仅为介绍方法的思想而做此简化。 ) 3散射角的确定 1.ca从(cosφ,1)抽样确定后,方位角为6,散射为y的出射方向所在直线的方程的确定 图3方位角示意图 图4散射角示意图 设OP“所在直线方向余弦为u,U*,w”,(w=0)则出射方向所在直线的方向余弦 为:( coswig, costart’,cosy)
OB直线方程: tgy 2.下面讨论OP·所在直线的方向余弦U,v·,W“(W*=0)的确定 v '=cos(a+8) 其中 sIna 3.方位角θ确定后,散射角y的变化范围 为求γ的变化范围,只需求出射方向直线和探测器所在圆柱相切的条件。 探测器所在圆柱的方程 v(z-x)-u(y-y)]2+[(x-x0)-(z-x0)]2+[u(y-y0)-v(x-x)}2=R2 将(6)式变为 y=tg> z,令 (8) 将(8)式代入(7)化简得关于z的一元二次方程 (v-wb)2x2+2( +(ua -u)222+2(wa-u)(uzo-wxo)z+(uzo-wxo)2 +(ub-va)222+2(ub-va)(uxo-uyo)z+(0xo-uyo)2-R2=o (9) 令Δ=0得关于a,b的方程 A 2+B62+ Eab Ca+D6 +F=0 (10) 其中A,B,C,D,E,F可计算求得: 将a=1gy·sin(a+θ)b=ty·cos(a+)代入(10)得 A,(cos0, sin])tgy +B,(cos], sin)tgy+C1=0 解之得gy=f1(0)tgy=f2(6) 进而可以确定cosy=[f(0),2(8)的变化范围。 4统计概率的积分计算 Ymax 考虑积分P= P2(6,y)ddo1,2 不妨设光子在碰撞点发生 Compton 1散射,散射前后能量分别为E,E E E compton散射满足Kein- Nishina散射规律,概率密度函数为 f(x/E)=K(E) rE+1-xy21_1+1]1≤x≤1+2E (11) 其它
其中K(E)=[1-2(E1]1n43)+1/2+4/a-nn.1a u为散射角余弦cy叫=1-E+Ex COS (12) 将(12)代人(11) 2(E+1) K(E)L E+1-E 这时 h1K(E)1acwy+b(aoy+b)(acwy+b“(E))÷l E dcosB a=-E,b=E+1,(E)=1-2(E+1 其中 d(E)=( E 显然内重积分是可积的,这样就可以变二重积分为单重积分,使数值化处理变得非常简 单 讨论 尽管可以对圆柱张角进行数值化处理但这种种数值化处理所带来的计算费用的增加, 常常不如用非数值化处理方法好,比如当粒子碰撞点远离圆柱到一定程度时,用什么数值化 处理方法都不如用非数值化处理方法好,而且随着距离不同差别会很大,下一步将研究随碰 撞点距离圆柱不同而不同的处理方法。 参考文献 〔1〕裴鹿成、张孝泽.蒙特卡罗方法及其在粒子输运问题中的应用,北京:科学出版社,1980年 [2〕裴鹿成等.计算机随机模拟,长沙湖南科学技术出版社,1989年。 [3]Gardner P R, Choi H K, Mickael M., Yacout A.M., Jin Y, and Verghese K, Nucl Sci. and Eng, 95, 245 (1987) [4]Mickal M. Gardner R P and Verghesg K Nucl Sci. and Eng. ,99: 251-266(1998) [5]Garder R. P, Mickal M, omby M, and verghess K Nucl. sci, and Engg., 108: 240-246(1991) [6] Gareler R P, Mickal M and Verghese, Nucl. Sci. and Eng,,, 98: 51-63(1988) 〔7〕邓力蔡少辉、黄正丰、黄捷计算物理.1996.13(4) 27
非归一分布随机抽样方法研究 程锦荣)段香梅2)杜爱军1) (1)安徽大学物理系,合肥,2)中国科学院固体物理研究所,合肥23031) 摘要通过对非归一分布随机抽样的 Metropolis方法、期望方法、排队方法和舍选方法 的研究,给出了若个新的结论。 关键词非归一分布, Metroplis抽样方法,期望抽样方法,排队抽样方法,舍选抽样方 法 1引言 物理研究中,不少问题最终常归结为需要计算如下物理系综的平均观察量 A(元)丌(x)d (1) (彩)d 式中,来表示相空间a中的一个点,或一个状态,一个位形;丌(彩表示这个物理系综的分布; A(表示某个微观观察量。 由于(1)式是一个高维积分1,通常的数值方法难以完成对物理系综平均观察量(A)计 算,又由于物理系综分布x(x)仅满足非负条件而不满足归条件{2,即 (x)≥0 丌(x ≠ 故采用一般的 Monte Carlo方法对A进行估计,同样十分困难(2) 由于在物理及其它领域中,存在着大量的类似问题。因此,采用有效方法解决由已知非 归一分布的随机抽样问题,就成了计算物理中急需解决的重要课题之 解决非归一分布随机抽样问题的第一个方法是由N. Metropolis等人于1953年给出的 ,后来人们称之为 Metropolis方法。由于 Metropolis方法存在着若于明显的缺点(2,国内外 学者先后提出了 Barker方法41、 hastings方法、 Heat Bath方法0、Pas法2、期望方法12 和排队方法2等一系列解决非归一分布问题的随机抽样方法 本文从实际计算的角度,选择抽样分布、收敛速度、抽样效率、花费CPU时间、抽样中进 入新状态的样本数目和子样的随机性质等指标,对非归一分布的 Metropolis方法、期望方法、 排队方法和舍选方法进行了评价。数值计算结果表明期望方法和舍选方法明显地优于 Metroplis方法和排队方法,且期望方法更适合于抽取大样本,舍选方法更适合于抽取小样 本 2抽样方法 国家高性能计算基金资助