ピークの情報のフィッティングによる導出

このスペクトルの x 軸はエネルギーを示していますが、 ここでは既知のエネルギーの$\gamma$ 線がどのチャンネルに来ているかが わかります。 $^{137}$Cs のデータと $^{60}$Co のデータをほぼ同じ条件で得た場合、 顕著に見られるピークがそれぞれのスペクトルで2本、合計4本見られます。 これらのエネルギーとチャンネルの関係からx軸の較正を行うことができます。

gnuplot は version 3.7 から、マルカール法を使った 非線型フィッティング が使えるようになりました。 詳細については文献 ¹⁰⁴に譲ります。

コマンドファイルとして、 text_test04.gp を使います。 まずはピークの位置がほかの成分から別れている $^{137}$Cs のスペクトルを見てみましょう。 230 ch. 辺りのピークはほとんど正規分布していますので、

$$f_1(x) = a_{11} \exp \left(-\frac{(x-a_{12})^2}{2a_{13}^2}\right)$$ (1)

という関数を定義し、このうち変数 $a_{11}, a_{12}, a_{13}$ を 最適化しましょう。 非線形フィットを行う場合には初期値が必要となります。 これは実際にスペクトルを表示して適当な値を最初に代入する必要があります。 また、スペクトルの全体をフィットするわけではなく、 ある範囲に限ってフィットするわけですから、 フィットのコマンドに際して範囲を指定する必要があります。 さらに、誤差としてデータの平方根を使うことにします。 gnuplot のコマンドファイルは次のようになります。

f1(x) = a11*exp(-((x-a12)/a13)**2/2)
a11 = 2000 ; a12 = 230 ; a13 = 10
fit [210:240] f1(x) "text_test04_cs137.dat" using 1:2:($2**0.5) via a11,a12,a13

この部分のフィッティングの結果は次のようになります。

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a11             = 1989.14          +/- 21.14        (1.063%)
a12             = 225.805          +/- 0.06535      (0.02894%)
a13             = 7.25809          +/- 0.06076      (0.8371%)

このようにして、ピークの位置と高さ、幅の情報を得ることができます。 右に表示されているのは、それぞれのパラメータの決定精度（誤差）です。 フィッティングの結果は、

plot [0:300] f1(x), "text_test04_cs137.dat"

とplot コマンドを用いて表示させます。

ところが、データによっては初期値を工夫してもうまくフィッテングできない場合があります。 これは、関数がデータと合っていない場合もありますが、パラメータの大きさが非常に違う場合にも起こります。 上の例では、a11とa13の比は300程度でした。この比が大きくなるとgnuplotは一方を0と判断してエラーを返すことがあります。 また、幅a13がピーク位置a12に比べて非常に小さい場合にも、フィッテングがうまくできない場合があります。 このようなときは、

f1(x) = a11*100*exp(-((x-a12-200)/a13)**2/2)
a11 = 20 ; a12 = 30 ; a13 = 10
fit [210:240] f1(x) "text_test04_cs137.dat" using 1:2:($2**0.5) via a11,a12,a13

のように、定数倍や定数値のオフセットをつけてa11,a12,a13の値を同じ桁にするとうまくいきます。

次にK-X 線のピークのフィットを試みましょう。 K-X 線のピークの部分には、 ほぼ一定の値を持つピーク成分以外の成分があります。 先ほどのような単純なGauss 関数ではなく、

$$f_2(x) = a_{21} \exp \left(-\frac{(x-a_{22})^2}{2a_{23}^2}\right) + a_{24}$$ (2)

という、ピーク成分以外を示す一定値を加えた関数にしてみましょう。 先ほどと同様におおよその値を初期値にして、範囲を指定して fit を 実行してみましょう。

f2(x) = a21*exp(-((x-a22)/a23)**2/2)+a24
a21 = 1400 ; a22 = 15 ; a23 = 3 ; a24 = 200
fit [10:16] f2(x) "text_test04_cs137.dat" using 1:2:($2**0.5) \
via a21,a22,a23,a24

図 15: $^{137}$Cs のフィッティングの結果の描画

次に$^{60}$Co のスペクトルですが、2本のピークはお互いが 少なからず重なっており、さらに コンプトン端と重なる部分もあり、 $^{137}$Cs のスペクトルの時に比べると 初期値を与えることも簡単ではありません。 gnuplot では、定義した変数の値が保持されるので、 この特徴を活かして、順次初期値を変更しながら フィッティングを実行してみましょう。

まず、 ピークを2本のガウス分布で表現し、 コンプトン散乱等の成分を部分的に2次関数で表現することにしましょう。 $f_3$ は 1.332 MeV のピーク、 $f_4$ は 1.173 MeV のピーク、 $g_1$ は コンプトン散乱等の成分として、

$$f_3(x) = a_{31} \exp \left(-\frac{(x-a_{32})^2}{2a_{33}^2}\right)$$ (3)

$$f_4(x) = a_{41} \exp \left(-\frac{(x-a_{42})^2}{2a_{43}^2}\right)$$ (4)

$$g_1(x) = b_{1}(x-b_{2})^2+b_{3}$$ (5)

とおきます。 最終的な目標はこれらの総和である、

$$h(x) = f_3(x) + f_4(x) + g_1(x)$$ (6)

を9変数( $a_{31}, a_{32}, a_{33}, a_{41}, a_{42}, a_{43}, b_{1}, b_{2}, b_{3}$) を使ってスペクトルに最適化するわけです。 実際にやってみればわかりますが、 適当な初期値を与えないと簡単にはフィッティングは収束しません。

まず、手順ですが、次のように考えてみました。

(1)

425 ch. から470 ch. のデータに対して $f_3(x)$ をフィッティングし、 $a_{31}, a_{32}, a_{33}$ の値をより初期値として適当な値にする。

(2)

380 ch. から410 ch. のデータに対して $f_4(x)$ をフィッティングし、 $a_{41}, a_{42}, a_{43}$ の値をより初期値として適当な値にする。

(3)

300 ch. から360 ch. のデータに対して $g_1(x)$ をフィッティングし、 $b_{1}, b_{2}, b_{3}$ の値をより初期値として適当な値にする。

(4)

得られた変数値を初期値にして$h(x)$を9変数でフィッティングする。

(1-3) で得た初期値をプロットしてみたところ図16 のようになってしまい、必ずしも良い初期値とは言えません。 この原因は(3) の初期値に大きく依存していると考えられます。 $b_1, b_2, b_3$の部分は(3)の手順で自動的に取得すると、 思わぬ値を入れてしまう可能性があります。 ここでは(3) の手順に頼らず、

b1=0.01 ; b2 = 500 ; b3=1
#fit [320.0:360.0] g1(x) "text_test04_co60.dat" via b1, b2, b3

のようにフィッティングの部分をコメントアウトして、 手作業での入力が初期値になるようにしましょう。

図 16: 初期値をきめるプロセス、 左の図は $f_3(x), f_4(x), g_1(x)$とも個別にフィットして得た初期値、 右の図は $f_3(x), f_4(x)$はフィットして得た初期値、 $g_1(x)$については手作業で入力したもの。

個別にフィッティングを行ったため、1.173 MeV のピークの部分が 過大に評価されている様子がわかります。

さて、これらの初期値を使って最終目標である全体のフィットを 行ってみましょう。 コマンドはかなり長くなりますが、次の通りです。

fit [300.0:500.0] h(x) "text_test04_co60.dat" using 1:2:(sqrt($2)) \
via a31,a32,a33,a41,a42,a43,b1,b2,b3

plot [200:550] h(x), f3(x),f4(x),g1l(x),\
"text_test04_co60.dat" using 1:2:($2**(0.5)) with yerrorbar

結果は全体を示す関数 $h(x)$ と 個別の関数とを、データとともに示すように指示されています。

結果は図17の通りで、2つのピークの部分は 良く再現されています。 念のため書き加えますが、 フィッティングの結果はあくまで gnuplot の持つフィッティングルーチンの出した答えであって、 その妥当性は解析する人が考える必要があります。

図 17: 全体のフィッティングの結果