Peak Fit With Gnuplot Part3 複数ピーク＋核種データ編
Peak Fit With Gnuplot 最初の一歩とＫ４０編、
Peak Fit With Gnuplot Part2 任意の単一ピーク編の続編です。
他に、
Peak Fit With Gnuplot Part4 検出器の特性調査編、
Peak Fit With Gnuplot Part5 検出器の自動特性調査編、
Peak Fit With Gnuplot Part6 色々な関数編があります。
複数ピークを扱えないと、セシウム１３４や１３７のピークもろくに調べられないので、
単一ピークで徐々に実験したことを基にして、幾つかの要素を加えます。
また、複数ピークを扱う場合、色々な核種の分岐率というか放出率というか、そういう情報とエネルギーとかも
取り扱えると嬉しいので、その辺も工夫します。
さらに、分解能は直線で誤魔化せそうですが、効率の方は、さすがに直線ではマズイだろう、ということで、
効率曲線をどうするのか、考えます。
現在実験中のスクリプトがこれ。
ＰｙｔｈｏｎとかからＧｎｕｐｌｏｔを呼ぶのでもなく、ＧｎｕｐｌｏｔからＡｗｋやらＳｈｅｌｌｓｃｒｉｐｔを呼ぶのでもなく、
「Ｇｎｕｐｌｏｔの機能だけ」で、何でもやってしまう、変態的スクリプトです。
まず、f にスペクトルのファイル名を入れたら、パート２で作った準備スクリプトを走らせて、変数とかを色々と設定します。
例：
f='ThereminoMCA_2015-04-01-Uso.txt'
load 'ft0.plt'
そしたら、続けて、このスクリプトを走らせるだけ。
load 'nl.plt'
ただし、まだ実験中なので、どの範囲でフィットする、とかは、スクリプトの中の
fmin と fmax に適当な値を入れて試すことになります。
また、初期値とか、あるいは、スペクトルによっては、Ｇｎｕｐｌｏｔが計算を終えるまでに
１００回以上もぐるぐると回って、時間がかかったり、結果がまとまらない可能性もあります。
とにかく実験中のスクリプトですので、あんまり期待したり、怒ったりなさっても、私は全く感知致しませんので、あしからず。
#========================== # Automatic peak hight fitting test # 自動ピーク高フィッティングのテスト # エネルギー位置は、データからの固定、分解能は、標準的なＮａＩ３インチの式をそのまま # で、残ったピーク高のみのフィッティングをしてみます。これは、下手にピークの中央や # 分解能もフィットすると、とんでもない位置になったり、変な分解能のピークが乱立したり、という # そういう可能性もありますし、動かす変数が少ないほうが、フィットしてくれるＧｎｕｐｌｏｔさんも # 楽なんじゃないだろうか、という心遣いからです。
# List of Radio Isotopes # 核種データ。 # ＳＷ８３Ａ様のページからコピペしたものとか、Ｉｄｅｎｔｉｆｙ．ｅｘｅの核種データからのものとか # 混じっていまして、微妙に値も違っているのですが、まあ、テスト用としても、素人用としても # 特に問題はないだろうと思います。 # まだ、少ししか含めていません。
# 最初は、ピーク高を調べる核種のリスト。これに、好きな核種をいれます。 # 一番最後のものが、現在やっているＥｕ１５２と鉛のＸ線の指定。 nl='Cs134 Cs137 Bi214 Ann' # K40' #nl='Ann Cs134 Cs137' nl='Cs134 Cs137' nl='Cs137 Bi214 Cs134' nl='Cs137 Bi214 Ann' nl='Eu152 Pbx'
# ここからが、核種データ # まず、核種の名前の変数に、核種の名前の文字列と、その核種のどのピークをフィットするかを指定。 # 単一ピークの核種は、１を入れます。複数のピークがある場合、好きなピークを選んで並べます。 # 次に、核種の名前に続いて、アンダースコア（下線）とピーク番号の変数に、それぞれのピークの # エネルギーと、放出率が入ります。「放出率」という日本語で正しいのか、ちょっとまだ確認してませんが、 # 要は、その核種がその崩壊の仕方をする割合と、その崩壊に於いて、そのエネルギーの光子が # 飛ぶ割り合いを掛け合わせたもので、Ｉｄｅｎｔｉｆｙ．ｅｘｅとかの核種情報もこの形になってるみたいです。
# ちなみに、ピーク番号に１を指定したピークを基準ピークとみなし、それを他のピークの関数は、 # このピークの放出率とお目当てのピークの放出率の割合に、効率曲線の数字を噛ませて高さを計算します。 # これは、現在のテストの仕方で、もっと自由にすることも可能です。
# Name of RI and which photopeaks to use # and Photo peaks for each RI. Energy and ratio
# 鉛のＸ線と５１１ｋｅＶの場合、放出率は、まあ、関係ないと思うので、１００を入れときました。 Pbx='Pbx 1' Pbx_1='77 100'
Cs137='Cs137 1' Cs137_1='661.66 85.1' Cs137_2='32.19 3.77' Cs137_3='31.81 2.05' Cs137_4='36.35 1.04'
Cs134='Cs134 1 2 3 4 5' Cs134_5='563.2 8.4' Cs134_3='569.3 15.4' Cs134_2='604.7 97.6' Cs134_1='795.9 85.5' Cs134_4='802.0 8.7'
K40='K40 1' K40_1='1460.83 99.53'
Bi214='Bi214 1' Bi214_1='609.31 46.09' Bi214_2='1120.27 15.04' Bi214_3='1764.49 15.92'
Ann='Ann 1' Ann_1='511 100'
Eu152='Eu152 1 2 3 4 5 6 8 9 10 11 12 13 14 15' # EC 13.33 Y Eu152_2 = '39.52 21.18' Eu152_3 = '40.11 38.39' Eu152_4 = '45.37 11.12' Eu152_5 = '121.77 28.38' Eu152_6 = '244.69 7.507' Eu152_1 = '344.28 26.59' Eu152_8 = '411.12 2.233' Eu152_9 = '443.89 2.802' Eu152_10 = '778.91 12.98' Eu152_11 = '867.38 4.210' Eu152_12 = '964.11 14.46' Eu152_13 = '1085.88 9.941' Eu152_14 = '1112.07 13.57' Eu152_15 = '1408.00 20.85'
# 効率曲線。これもまだ、幾つかテスト中で、良く分かっていません。 # Efficiency curve # ln(eff)=a1+a2*ln(E) ef(x) = exp(a1+a2*log(X))
ef(x) = e1*x+e2 e1=-0.05 e2=100
ef(x) = e1 * x**e2 e1= 1.895; e2= -0.62
# これは、使うかもしれないガウス関数の一つ。 gg(x,hn,En,sn)=hn*exp(-((x-En)/sn)**2/2)
# マクロの使用を許可 # 最初は、マクロでやろうと思ったら、何故か上手くいかなくて、結局ｅｖａｌ文を使いました。 set macro
# 全てのピークのガウス関数を含んだ、まとめ用関数を定義する為の文字列を初期化 gallcmd='gall(x)='
# フィットする時に、どの変数群を使うか、あるいは固定しておくかによって、Ｆｉｔコマンドの # Ｖｉａの項目が変わってくるので、それに対応するための下準備 vall='' vsall='' vhall='' veall=''
# その核種に複数のピークをフィットするかどうかのフラグ用変数 e1e2=0
# 外側のループ。核種ごとの作業がこのレベルで行われます。 do for [i=1:words(nl)]{
# ri に核種の名前が入ります。 ri=word(nl,i)
# それをプリント print ri
# その核種用のガウス関数の定義をする為の文字列を初期化 gcmd='g'.ri.'(x)='
# その核種用のＶｉａの項目を初期化 vcmd='' vscmd='' vhcmd='' vecmd=''
# 全ての核種のまとめ関数に、この核種のガウス関数を追加 gallcmd=gallcmd.'+g'.ri.'(x)'
# ｗｒｉ という変数に、その核種のピークの数を入れる為のコマンド文字列を作成し、ｅｖａｌ文で実行。 cmd='wri=words('.ri.')' eval cmd
# この変数なんだったか、早くも覚えていない・・・ rphr=1
# 内側のループ。この中で、ピークごとの作業をやっています。 do for [rn=2:wri]{
# ｒｎｎ という変数に、ピーク番号をセットするコマンドの文字列を作って、それを実行。 cmd='rnn=word('.ri.',rn)' eval cmd
# ｐｈという変数に、そのピークの名前の文字列（核種の名前＋下線＋ピーク番号）をセット ph=ri.'_'.rnn
# 色んな変数をセットするコマンドの文字列を作って、それを実行。 # まず、eCs137_1 という風に、ピークの名前の前にｅを付けたものに、そのピークのエネルギーを入れ、 ｅｐｈにもその値をコピー。 # rCs137_1 という風に、ピークの名前の前にｒを付けたものが放出率で、ｒｐｈにもコピー cmd='e'.ph.'=word('.ph.',1);eph=e'.ph.';r'.ph.'=word('.ph.',2);rph=r'.ph eval cmd
# ｓｐｈにそのピークのシグマ（ピークの幅を指定する値）をエネルギー位置から計算して入れる。 sph=SigE(eph)
# sCs137_1 という風に、ピークの名前の前にｓを付けた変数にもその値を入れる為のコマンドを作って、実行。 cmd='s'.ph.'=sph' eval cmd
# ピーク番号が１の場合は、そのピークの放出率とピーク名を記憶 if(rn==2){ rph1 = rph ph1 = ph # ピーク番号が２以上の時は、ピークの高さを決定する為の割り合いを、番号が１のピークと該当ピークの放出率から計算 # また、複数ピークの核種だというフラグをセット。 }else{ rphr = sprintf("%f",rph/rph1) e1e2=1 }
# Ｓｔａｔコマンドで、そのエネルギーの左右１シグマの範囲での最高地点を見つけ、ｈｐｈという変数にセット stat [eph-sph:eph+sph] f every ::HL::BGL u 1:2 nooutput hph=STATS_max_y
# hCs137_1 という風に、ピークの名前の前にhを付けた変数に、ピークの高さの初期値をセット # この変数は、複数ピークの場合、２番目以降のピークについては、現在使われていませんが、後で使う予定。 cmd='h'.ph.'=hph;print ph," ",eph," ",hph' eval cmd
# フィットする核種が一つだけか、ピーク番号が１なら、以下の様なガウス関数を定義して、ＶｉａのマクロにhCs137_1を追加。 # gCs137_1(x)=hCs137_1*exp(-((x-eCs137_1)/sCs137_1)**2/2) if(wri==2 || rn==2){ cmd='g'.ph.'(x)=h'.ph.'*ef(e'.ph.')*exp(-((x-e'.ph.')/s'.ph.')**2/2) vhall=vhall.',h'.ph
# 複数ピークの２番目からはフィットする核種が一つだけか、ピーク番号が１なら、以下の様なガウス関数を定義。 # 高さの部分は、一番目のピークに放出率をかけたもの。 # gCs132_1(x)=hCs134_1*3.4566444*exp(-((x-eCs134_2)/sCs134_2)**2/2) }else{ cmd='g'.ph.'(x)=h'.ph1.'*'.rphr.'*ef(e'.ph.')*exp(-((x-e'.ph.')/s'.ph.')**2/2) } eval cmd
# 核種の全部の関数をまとめた関数の定義の為の文字列に追加。 gcmd=gcmd.'+g'.ph.'(x)' } # 核種の全部の関数をまとめた関数の定義を実行 eval gcmd } # 全ての核種の関数をまとめた関数の定義を実行 eval gallcmd
# 複数ピークがあったら、効率曲線の変数もフィットの時に決めてもらう。 if(e1e2>0){ vhall=vhall.',e1,e2' }
# 実験の残りのゴミですが、ちょっと保存 # fit [E-n*s:E+n*s] gall(x) f every ::HL using 1:2 noerror via h,h2,a,b
# バックグラウンドの曲線の定義とその初期値の設定。これをもっと工夫する予定。 b(x)=ba*x**2+bb*x+bc ba=5;bb=-10;bc=9600
#b(x)=ba*x**2+bb*x+bc+bd*x**3 #bd=1
# フィットする範囲の指定。ここは、まだ手動で実験中。これも、自動化する予定。 fmin=540 fmin=20 fmax=1600 fmax=550
# xrange も、まだ手動設定でテスト。 set xrange [0:1800]
# 色々なＦｉｔコマンドの例。Ｖｉａのところに変数を追加することで、その項目が変動します。なければ、その項目は（初期値で）固定。 # @vhall で、選択された核種のピーク高の変数や、必要なら効率曲線の変数が変動するようになってます。 #fit [fmin:fmax] gall(x)+b(x) f every ::HL::BGL using 1:2 noerror via ba,bb,bc @vhall ,eCs137_1,eCs134_1 fit [fmin:fmax] gall(x)+b(x) f every ::HL::BGL using 1:2 noerror via ba,bb,bc @vhall #,eCs137_1 #fit [fmin:fmax] gall(x)+b(x) f every ::HL::BGL using 1:2 noerror via ba,bb,bc @vhall ,eCs137_1
# 色々なＰｌｏｔコマンド。 #set xrange [250:1300] #plot f every ::HL::BGL using 1:2 with dot, [fmin:fmax] gall(x)+b(x),[eAnn_1-3*sAnn_1:eAnn_1+3*sAnn_1] gAnn_1(x), [fmin:fmax] b(x) #plot f every ::HL::BGL using 1:2 with dot, [fmin:fmax] gall(x)+b(x),gAnn_1(x),gCs137_1(x),gBi214_1(x), [fmin:fmax] b(x) plot f every ::HL::BGL using 1:2 with dot, [fmin:fmax] gall(x)+b(x),[fmin:fmax] b(x)
# ｑ は、Ｑｕｉｔの略で、この後は実行されません。昔のフィットコマンドとかを置いてあるだけです。 q fit gall(x)+bb*x+bc f every ::HL using 1:2 noerror via ba,bb,bc @vhall fit gall(x) f every ::HL using 1:2 noerror via ba,bb,bc @vhall
#====================
上記のスクリプトで、Ｅｕ１５２線源のスペクトルの低い方をフィットしてみた例。
まだまだではありますが、少し近くなってきました。
２０ＢｑくらいのＣｓ１３７の汚染がある日本の土、２０ｇのスペクトルで再テスト。
After 6 iterations the fit converged. final sum of squares of residuals : 5.41401e+006 rel. change during last iteration : -5.68351e-007
degrees of freedom (FIT_NDF) : 220 rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 156.873 variance of residuals (reduced chisquare) = WSSR/ndf : 24609.1
Final set of parameters Asymptotic Standard Error ======================= ========================== ba = 0.00256898 +/- 0.0007405 (28.82%) bb = -8.22867 +/- 1.141 (13.86%) bc = 7048.43 +/- 428.7 (6.083%) hCs137_1 = 8276.16 +/- 49.13 (0.5936%) hBi214_1 = 1138.19 +/- 60.57 (5.322%) hCs134_1 = 1222.67 +/- 40.69 (3.328%)
correlation matrix of the fit parameters: ba bb bc hCs137 hBi214 hCs134 ba 1.000 bb -0.996 1.000 bc 0.978 -0.993 1.000 hCs137_1 0.272 -0.230 0.161 1.000 hBi214_1 -0.397 0.419 -0.437 -0.163 1.000 hCs134_1 0.367 -0.338 0.281 0.456 -0.569 1.000
まだ、ちょっとと言うか、色々と問題があり、スクリプトも随時改造中です。
フィットしたピークをそれぞれベースラインの上と、画面の下に全部出す改造と、
横軸がズレちゃっている時でも、再調整しながらフィットする改造を施しました。
つづきはこちら。検出器の特性調査編