-
寿命是放射性核素的一个基本性质。半衰期的长短反映了放射性核素衰变的快慢,提供了原子核内部结构的信息。定量地计算寿命不仅有助于揭示精确的核素合成路径,同时还有助于研究人们一直讨论的问题:核衰变常数是如何变化的[1]。目前用于核素寿命测量的实验技术多种多样,用数据“重现”密度函数并提取寿命是一种最为直接、使用最为广泛的测量方法。换句话说,通过实验测得的衰变时间得到直方图,并用指数函数最小二乘法拟合直方图获得寿命。尽管该方法在各方面研究中都有较为广泛的应用,但当衰变事件较少且无法形成合适的直方图时,最小二乘法拟合可能会得到不合理的结果。本文通过模拟产生不同数量的衰变事件,研究了该方法的适用范围。为了方便区分,在本文中我们把该方法称为直接拟合法。
为了有效地从稀疏数据(衰变事件较少)中准确提提取半衰期,Schmidt等[2] 提出了一种巧妙的方法来处理稀疏数据。该方法是按对数坐标“压缩”衰变时间,即对衰变时间取对数,从而使得衰变时间分布更为集中。随后用最小二乘法指数拟合这些“压缩”衰变时间的直方图,得到寿命。该方法被核谱学学家广泛采用,尤其是超重核数据分析[3-5]。在本文中,该方法称为对数时间法。此外,当计数只有个位数时,可用极大似然法求得寿命及误差[6]。
当观测时间不受限制时,以上方法可以有效地应用于寿命计算中。但当观测时间受限时,Schmidt的方法就有一定的局限性,或会得到一些不合理的结果。正如其文章中所提到的该方法的适用条件是“...观察时间不受限制...”。为了求得当观测时间受限时的寿命,理论推导出了一套有效法的计算方法。该方法是基于指数函数的概率密度函数,用极大似然法推导得到的。该方法已经成功地应用于计算有限观测时间窗口、低计数率条件下全剥离
$ ^{^{94m}}{\rm{Ru}} $ 的寿命,且在核谱学衰变研究中有广泛的潜在应用。另外,在某些偶然情况下,最多只能计算一段持续时间内衰减事件的总数,而以上方法都不适用。因此,必须开发更适用的数据分析方法来处理这些实验数据和提取寿命信息,即可以基于极大似然估计推导出计算寿命的方法。在本文中,将以上方法称为观测时间受限时的极大似然法。以上几种方法是目前较为有代表性的寿命的计算方法,更多的计算方法研究可以参考文献[7]。目前还没就这些方法的适用范围进行较为详细且具体的限定或给出可参考的适用范围,尤其是不同方法计算标准偏差的估计。本文基于模拟数据,研究了在不同计数或不同的测量时间窗口条件下,这四种方法的适用条件和范围。数据模拟中选取全剥离
$ ^{{94m}}{\rm{Ru}} $ 作 为目标核素,其寿命为137(7) μs,当洛伦兹因子$ \gamma $ 为1.302 时(与实验设定条件一致),得到它在实验室坐标系下的寿命为178(10) μs。 -
产生随机数的方法有很多种,本文中使用了一种新的产生方法。该方法不仅适用于产生指数衰减的随机数,还适用于产生大部分已知概率密度函数的随机数。本方法的基本思路是:由于概率密度函数的累积分布(应变量
$ 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)所示。 -
当衰变事件在统计量较高、观测时间窗口不受限的情况下,把测量得到的衰变时刻画成直方图。用指数函数通过最小二乘法直接对直方图进行拟合,拟合结果如图1(b)所示。
$ ^{94m} {\rm{Ru}} ^{44+} $ 的寿命$ \tau $ 及其误差$ \sigma $ 可通过拟合得到。图1(b)中的红色实线为拟合曲线,拟合得到的寿命为178.2(4) μs,与设定值178.0 μs 在误差范围内符合。同时我们还模拟产生了50,60,···,100,150,200,···,10 000个衰变时刻,分别用指数函数进行拟合,计算得到的结果如图2中的红点所示。从图2中我们可知,当衰变事件的个数大于70 时,通过拟合也能得到较为准确的结果,且寿命的标准偏差$ \sigma $ 约为20 μs。当衰变事件的个数增多时,寿命的准确度及精度都有所提升,尤其是精度(寿命的标准偏差$ \sigma $ )有较为明显的提升。图2(b)中的虚线为寿命设置值的标准偏差$ \sigma $ = 10 μs),当衰变事件的个数在400 左右时,寿命的标准偏差将小于10 μs。根据此模拟结果,可以估计出在不同精度下所需的衰变事件的个数。当个数小于70 时,用此方法得到标准偏差($ \sigma $ )较大。接下来我们介绍,当衰变事件的个数少于70 时计算半衰期的方法。 -
当衰变时刻的统计量只有几十个甚至更少时,就无法形成合适的直方图,用直接拟合法对数据进行拟合可能会得到不合理的结果。为了能有效地计算统计量较少情况下的半衰期,人们提出了两种计算方法。一种是对衰变时刻取对数,使得衰变时刻分布更为集中,从而有效计算出寿命。另外一种是基于极大似然法严格推导出的计算公式。两种方法在文献[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所示。通过模拟结果可知,该方法有较为广泛的适用范围,不同计数和不同观测时间窗口下,都能得到一个较为准确的结果。
-
在兰州重离子加速器冷却储存环上开展了全剥离同核异能态
$ ^{94m} {\rm{Ru}} ^{44+} $ 的寿命测量实验。实验中储存环设置在等时性模式下,光学设置的转换点$ \gamma_{\rm{t}} $ [8]为1.302,等于目标离子($ ^{94m} {\rm{Ru}} ^{44+} $ )的洛伦兹因子$ \gamma $ ,即$ \gamma = \gamma_{\rm{t}} = 1.302$ 。通过分析实验数据分析,在观测时间窗口(20 μs,180 μs)内观测到了49个衰变事例。这49的衰变事例用对数时间法计算得到$ ^{94m} {\rm{Ru}}^{44+} $ 的寿命为88(19) μs,与模拟结果94(13) μs 接近,但与理论值178(10) μs 相差较大,如图5所示。此外,通过磁场修正,在观测时间窗口内鉴别出了$ m(m = 29)$ 个未衰变事件,以及剩余了$ n (n = 37)$ 个衰变事件(原来49个)。实验中的观测时间窗口大小$ T $ =160 μs,37个衰变事件的平均衰变时间$ \bar{t} $ =66.9 μs,$ n $ =37,$ m $ =29,代入观测时间受限时的极大似然法中的式(13)、式(14)可求得寿命为192(32) μs,与模拟结果169(34) μs及理论预期值178(10) μs在误差范围内符合,如图5所示。实验及数据分析的详细过程可参考文献[9]。 -
本文通过模拟数据,研究了四种计算寿命的方法。前三种计算方法的观测时间窗口不受限,根据模拟结果,非常直观地展示了使用不同方法得到的寿命及误差随衰变事件个数的变化情况,并给出了这三种方法的适用范围。第四种方法的观测时间窗口受限,同样通过模拟结果发现了观测时间受限时的极大似然法有较为广泛的适用性。通过
$ ^{94m} {\rm{Ru}} ^{44+} $ 的实验数据来进一步验证两种方法的适用范围及模拟结果的可靠性。本文的模拟结果,不仅为将来寿命测量实验的方案设计、寿命计算方法的选择提供参考,还为将来开展短寿命核的寿命测量(且不局限于短寿命核素的寿命测量)提供了实验参考和数据处理经验。
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 -