-
产生随机数的方法有很多种,本文中使用了一种新的产生方法。该方法不仅适用于产生指数衰减的随机数,还适用于产生大部分已知概率密度函数的随机数。本方法的基本思路是:由于概率密度函数的累积分布(应变量
$ y $ )服 从(0,1)的均匀分布,通过公式转换成自变量$ x $ (或$ t $ )的随机数,即为对应函数的随机数。以下我们产生本文中所用的指数衰减的随机数(衰变时间),并展示该方法的操作流程。模拟中,我们产生了200 000个均匀分布的随机数
$ y $ ,分布如图1(a)所示。假设衰变时刻$ t $ 是指数分布的自变量,衰变时刻$ t $ 指数分布的累积分布函数:$$ F(t) = 1 - {\rm{e}}^{-\lambda t}, \;\;\;\;0\leqslant t < \infty\,, $$ (1) 式(1)中,
$ F(t) $ 服从(0,1)的均匀分布。假设$ y = F(t) $ 为(0,1) 均匀分布的随机数。则式(1)可变为$ y = 1 - {\rm{e}}^{-\lambda t} $ ,转换可得:$$ t = -\frac{{\ln}(1-y)}{\lambda},\;\;\;\; 0< {y} < 1{\text{。}} $$ (2) 将均匀分布的随机数
$ y $ 通过式(2)转换成衰变时刻$ t $ 。当$ \lambda = \frac{1}{\tau} $ ($ ^{94m}{\rm{Ru}}^{44+} $ 在相对论坐标系下的寿命$ \tau $ = 178 μs)时,可得到衰变时间$ t $ 的分布如图1(b)所示。 -
当衰变时刻的统计量只有几十个甚至更少时,就无法形成合适的直方图,用直接拟合法对数据进行拟合可能会得到不合理的结果。为了能有效地计算统计量较少情况下的半衰期,人们提出了两种计算方法。一种是对衰变时刻取对数,使得衰变时刻分布更为集中,从而有效计算出寿命。另外一种是基于极大似然法严格推导出的计算公式。两种方法在文献[2, 6]中有详细的描述,本小节主要是用模拟数据来确定两种方法的适用范围。尤其是当衰变时刻的个数不同时,以本次实验中的情况为例,研究两种方法可以适用的范围。根据两篇文献中寿命计算公式的特点,我们通过“压缩”衰变时间的方法研究了当衰变事件的个数为10,20,···,100,200,300时,
$ ^{94m} {\rm{Ru}} ^{44+} $ 的寿命及其标准偏差($ \sigma $ )随衰变事件个数的变化情况,在本文中将此方法称为对数时间法。以及通过基于极大似然法得到的计算公式,研究了当衰变事件的个数为3,4,···,20时,$ ^{94m} {\rm{Ru}} ^{44+} $ 的寿命及其标准偏差($ \sigma $ )随衰变事件个数的变化情况,在本文中将此方法称为极大似然法。基于指数衰减公式
$ N = N_0 {\rm{e}}^{-\lambda t} (0\leqslant t < \infty) $ ,可得当衰变时刻$ t $ 取$ {\rm e} $ 的对数时,其对应的密度分布函数为$$ \frac{{\rm{d}}N}{{\rm{d}}({\ln}t)} = \frac{{\rm{d}}N}{{\rm{d}}t}\frac{{\rm{d}}t}{{\rm{d}}({\ln}t)} = -N_0\lambda t {\rm e}^{-\lambda t}, \; \;\; 0\leqslant t < \infty{\text{。}} $$ (3) 当我们令
$ \theta = {\ln} t $ ,则式(3)变为$$ \Big|\frac{{\rm{d}}N}{{\rm{d}}\theta}\Big| = N_0\lambda {\rm{e}}^{\theta-\lambda {\rm{e}}^{\theta}},\; \; -\infty < \theta < +\infty\,, $$ (4) 其中常数
$ N_0 $ 为初始时刻的计数,$ \lambda $ 为衰变常数。由式(4) 所得的曲线,即
$ \frac{{\rm{d}}N}{{\rm{d}}\theta} $ 随$ \theta $ 的变化曲线,是一个钟形,略不对称的曲线。为了求得衰变常数$ \lambda $ ,将衰变时刻$ t $ 取对数变为$ \theta $ ($ \theta = \ln t $ ),画出自变量为$ \theta $ 的直方图。用式(4)拟合直方图即可得到常数$ N_0 $ 和衰变常数$ \lambda $ ,即可得寿命$ \tau = \frac{1}{\lambda} $ 。在实验测量中,很难测量得到非常光滑的曲线,而只可以得到直方图。在这种条件下,
$ \theta $ 的标准偏差$ \sigma_{\theta_{\exp}} $ 为$$ \sigma_{\theta_{\exp}} = \sqrt{\frac{\sum\nolimits_{i = 1}^{N}(\theta_i-\overline{\theta}_{\exp})^2}{N}}, $$ (5) 其中
$ \overline{\theta}_{exp} $ 的表达式为$$ \overline{\theta}_{\exp} = \frac{\sum\nolimits_{i = 1}^{N}\theta_i}{N}{\text{。}} $$ (6) 如果把衰变时刻
$ t $ 取10为底的对数$ \lg t $ 时,则$ \vartheta = \lg t $ 的分布函数变为$$ \Big|\frac{{\rm{d}}N}{{\rm{d}}\vartheta}\Big| = {\ln}10\cdot N_0\cdot\lambda \cdot10^{\vartheta}\cdot {\rm{e}}^{-\lambda\cdot10^{\vartheta}}, \;\; -\infty < \vartheta < +\infty{\text{。}} $$ (7) $ \vartheta $ 标准偏差同样可由式(5)和式(6)求得。我们模拟了当衰变事件的个数为 10, 11, 12,···, 20, 30,···, 100, 200, 300 时,用以上方法计算得到寿命,并得到寿命随衰变事件个数的变化情况,如图2(a)中的黑点所示。寿命的标准偏差(
$ \sigma $ )随衰变事件数目的变化如图2(b)中的黑点所示。另外当观测时间窗口分别为178, 356, 534, 712, 890及1 068 μs(即1~6倍的寿命),且固定衰变事件的个数为50时,用式(7)计算得到了寿命随观测时间窗口的变化情况。模拟中设定的寿命为
$ \tau $ 为178 μs,衰变事件的个数固定为50个,模拟结果如图3 所示。图3中的每个数据点值是50次模拟寿命的平均值。从图3中我们可知,当观测时间窗口小于3 倍的寿命时,求得的寿命比设定值小。尤其当观测时间窗口为1倍的平均寿命时(与实验情况一致),求得寿命偏离设定值较大。此外,当衰变事件的个数为个位时,由传统的极大似然法计算寿命及上下限公式[6]如下:
$$ \hat{\tau} = \overline{t} = \frac{1}{n}\sum\limits^{n}_{i = 1}t_i, $$ (8) $$ {\tau}_{\rm u} = \frac{\overline{t}}{1-z/\sqrt{n}},\; \; \;\; {\tau}_{\rm l} = \frac{\overline{t}}{1+z/\sqrt{n}}{\text{。}} $$ (9) 其中:
$ t_i $ 为衰变时间;$ n $ 为衰变事件的个数;$ {\tau}_{\rm u} $ 和$ {\tau}_{\rm l} $ 分别为寿命$ \hat{\tau} $ 的上限和下限;定量$ z $ 与选取的置信水平有关。在本文的计算中$ z $ 取1,即置信度为68.3%。本文模拟了当衰变事件的个数为3, 4, ···, 20时,寿命及其标准偏差(
$ \sigma $ )随衰变事件的个数的变化,如图2中的蓝点所示。基于图2的模拟结果可知,当观测时间窗口不受限时,衰变事件的个数不同,所采用的方法也不同。当衰变事件的个数为70以上时,即当数据可以得到一个良好的指数分布柱状图时,可用直接拟合法求得寿命,且误差也较小。当衰变事件的个数不能累成指数分布直方图且大于15 时,依然可以通过对数时间法求得寿命。当衰变事件的个数小于15时,可用极大似然法求得寿命。以上数据仅为一次数据模拟得到的结果,可以作为参考。在实际实验中,可根据实验的数据特点灵活选取 合适的寿命计算方法。
-
以上方法的适用条件是观测时间窗口不受限或足够大。显然,这样的条件实际上只能充分延长测量时间。然而,在某些偶然情况下,只能测得一个特定持续观测时间窗口内的衰变事件且所测得的个数较少,此时对数时间法就不再适用。因此,必须开发新的计算方法来准确提取寿命。本小节将介绍当观测时间受限时,用极大似然法计算寿命的方法。
对于衰变时刻
$ t\in (T_1, T_2) $ ,衰减概率函数的表达式为式(10),而当测量时刻$ t $ 超过$ T_2 $ 时,其存活概率S量化为式(11)。$$ f(t) = \frac{1}{\tau} {\rm{e}}^{-(t-T_1)/\tau},\; \; \;\; T_1< t < \infty\,, $$ (10) $$ S = \int\nolimits^{\infty}_{T_2}f(t){\rm{d}}t = {\rm{e}}^{-(T_2-T_1)/\tau}, \; \;\;\; t = T_2\,, $$ (11) 其中常数
$ \tau $ 为寿命。极大似然函数可以写成:$$ \begin{split} L(t_i|\tau) &= \prod\limits^n_{i = 1}f(t_i)\prod\limits^m_{i = 1}S \\&= \left(\frac{ {\rm{e}}^{-(t-T_1)/\tau}}{\tau}\right)^n \left( {\rm{e}}^{-(T_2-T_1)/\tau}\right)^m{\text{。}} \end{split} $$ (12) 基于极大似然估计法,可得到极大似然估计值
$ \hat{\tau} $ ,表达式为$$ \hat{\tau} = \frac{1}{n}\sum\limits^{n}_{i = 1}t_i-T_1+\frac{m(T_2-T_1)}{n}\,, $$ (13) $$ \sigma_{\hat{\tau}} = \frac{\hat{\tau}}{\sqrt{(m+n)(1- {\rm{e}}^{-(T_2-T_1)/\hat{\tau}})}}{\text{。}} $$ (14) 其中:(
$ T_1, T_2 $ )为观测时间窗口;$ n $ 为该时间窗口内观测到的衰变事件的个数;$ t_i $ 为衰变事件的衰变时间;$ m $ 为衰变时间$ t $ 大于$ T_2 $ 的事例数目,即在观测时间内未衰变事件的个数。式(14)表示的是寿命的不确定度$ \sigma_{\hat{\tau}} $ ,其推导的详细过程可参考文献[7]。我们通过模拟数据来验证该方法的适用范围,如图4所示。通过模拟结果可知,该方法有较为广泛的适用范围,不同计数和不同观测时间窗口下,都能得到一个较为准确的结果。
Simulation Study of Lifetime Calculation Methods for Radioactive Nuclides
-
摘要: 根据几种常用放射性核素的寿命计算方法,通过模拟数据研究了直接拟合法、对数时间法、极大似然法、观测时间受限时的极大似然法等四种寿命计算方法的适用范围。当观测时间不受限时,研究了在不同计数下寿命计算方法的适用范围。当观测时间受限时,研究了在不同观测时间窗口下寿命计算方法的适用范围。模拟中选用全剥离离子94mRu44+作为目标核素,得到了不同计数及不同观测时间窗口下的寿命及其误差,并给出了四种方法的适用范围。94mRu44+寿命的模拟结果与在兰州等时性质量谱仪上获得的实验结果在一倍标准偏差范围内一致,从而进一步验证了寿命计算方法的适用范围及模拟数据的可靠性。该模拟结果可为寿命测量实验设计提供理论依据和参考。
-
关键词:
- 兰州重离子加速器冷却储存环 /
- 等时性质量谱仪 /
- 全剥离离子 /
- 寿命模拟 /
- 94mRu
Abstract: According to several common lifetime calculation methods of radioactive nuclides, the scope of applicable of the four calculation methods, which are named Method of Direct Fitting, Method of Logarithmic Time, Method of Maximum Likelihood and Method of Maximum likelihood when observation time windows is limited, are studied based on the simulation data. As the observation time window is limited or not, the applicable range of the lifetime calculation methods in different observation time windows and different counts are discussed. In simulation, fully stripped ion 94mRu44+ was selected as the target nuclide, the lifetime and error in different counts and different observation time windows are obtained, and the applicable range of the four methods is given. The experimental data of 94mRu44+ was obtained from the lifetime measurement experiment which is performed by using the Isochronous Mass Spectrometry (IMS) at the HIRFL-CSR facility in Lanzhou. The simulation results are consistent with the experimental results within one error bar, thereby it is further verified the applicable range of the calculation method and the reliability of the simulation data. The simulation results provide theoretical basis and reference for the design of the future lifetime experiments.-
Key words:
- HIRFL-CSR /
- isochronous mass spectrometry /
- fully stripped ions /
- lifetime simulation /
- 94mRu
-
-
[1] LITVINOV Y A, BOSCH F. Rep Prog Phys, 2011, 74: 016301. doi: 10.1088/0034-4885/74/1/016301 [2] SCHMIDT K H. Eur Phys J A, 2000, 8: 141. doi: 10.1007/s100500070129 [3] OGANESSIAN Y T, YEREMIN A V, POPEKO A G, et al. Nature, 1999, 400: 242. doi: 10.1038/22281 [4] OGANESSIAN Y T, ABDULLIN F S, BAILEY P D, et al. Phys Rev Lett, 2010, 104: 142502. doi: 10.1103/PhysRevLett.104.142502 [5] KHUYAGBAATAR J, YAKUSHEV A, DÜLLMANN C E, et al. Phys Rev Lett, 2014, 112: 172501. doi: 10.1103/PhysRevLett.112.172501 [6] SCHMIDT K H, SAHM C C, PIELENZ K, et al. Z Phys A, 1984, 316: 19. doi: 10.1007/BF01415656 [7] CHEN X C, ZENG Q, LITVINOV Y A, et al. Phys Rev C, 2017, 96: 034302. doi: 10.1103/PhysRevC.96.034302 [8] HAUSMANN M, ATTALLAH F, BECKERT K, et al. Nucl Instr and Meth A, 2000, 446: 569. doi: 10.1016/S0168-9002(99)01192-4 [9] ZENG Q, WANG M, ZHOU X H, et al. Phys Rev C, 2017, 96: 031303(R). doi: 10.1103/PhysRevC.96.031303