Gnuplotを使った解析
Kromek K102 MCA
MCAについて
原子核散乱では、Kromek社のMCAを使用しています。
https://www.kromek.com/product/k102-multichannel-analyser/
この手のMCAとしては使い勝手が良くコンパクトで設定も簡単なので、簡単な放射線計測にはとても便利です。コントロールするソフトは同社のKSpectと言うソフトを使います。
https://www.kromek.com/product/free-gamma-spectroscopy-software/
Windows用の無料ソフトなので自分のパソコンにインストールして使うこともできます。ピークフィットもしてくれたりしますが、結果に誤差をつけてくれないので、測定したデータをレポートにそのまま使うのには不十分で、自分でプログラムを組んだり他のソフトを使って解析をする必要があります。
データ構造
Kspectで測定したヒストグラムはspeまたはtxt型式で保存できます。どちらでもテキストデータなので、適当なテキストエディタなどで読み込みことができます。Gnuplotで解析する場合は、これをGnuplotでも読み込める型式に変換してあげる必要があります。以下のいずれの方法を使っても構いません。
ファイルを見ればわかると思いますが、speファイルは始めの10行はヘッダーで、ここにファイルの取得時刻や測定時間などの情報が書かれています。ヒストグラムの情報は11行目から4106行目までの4096行 (=12bit) です。各行がヒストグラムの各binの高さ情報になっています。txtファイルは、この4096行のヒストグラムの情報のみが保存されています。
txtで保存しておけば、そのままGnuplotで読み込むことができます。speで保存してある場合は、適当なプログラムを使って、この11行めから4106行目を抜き出してあげれば良いです。せっかくプログラムを使う場合は、Gnuplotで読みやすいように変換するしておくとあとの解析が楽になります。
サンプルファイル
サンプルファイルを置いておきます。下で解説するプログラムや、標準線源の137Csを使ってGe検出器で取得したヒストグラムを置いておきます。
データの変換
得られたデータファイルを137Cs_sample.spe、137Cs_sample.txtとします。まずはこれをGnuplotで扱いやすい形に書き換えて、137Cs_sampe.datという名前で保存することにします。
Gnuplotでそのまま読み込む
txt型式で保存したデータはGnuplotで直接読み込むことができます。
gnuplot> plot '137Cs_sample.txt' using ($0):1:(sqrt($1)) with yerr
gnuplotでは、行番号をusing 中の $0 表示することができます。またそれぞれのヒストグラムのbinの誤差はその平方根とすればbinに対して誤差棒も付けることができます。
AWKで変換する
ファイルの中身を一行ずつ処理するようなプログラムを書くときは、awkというプログラム言語を使うと楽です。awkは大抵のUNIX系のOS（MacOSやLinux）には初めから入っていますのでインストール不要です。Windowsには入れたことがないので解説できないです。他の方法を使ってください。
以下の一行コマンドをコピーしてターミナルに入力してください。文法はC言語に似ているので、Cを知っていればなんとなく理解できるのではないかと思います。ここでは、1列目をbin番号、2列目がヒストグラムの高さ、3列目をその平方根として出力します。あとはリダイレクト > でファイルに入力します。
> awk '{if(NR>10 && NR<4107) printf("%d %f %f\n",NR-10,$1,sqrt($1));}' 137Cs_sample.spe > 137Cs_sample.dat
あとは普通にGnuplotでヒストグラムを表示することができます。
gnupolot> plot '137Cs_sample.dat' using 1:2:3 with yerrorbar
Cで変換する
Cでファイルを読み込むのはAWKのよりは少し大変です。下記のプログラムを使えば読み込むことができます。
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char *argv[])
{
FILE *f;
char s[256];
int bin = 1;
if(argc!=2){
printf("Usage: ./conv <FILENAME>.txt\n\n");
exit(-1);
}
if((f = fopen(argv[1],"r")) == NULL){
printf("file open error!\n");
exit(-1);
}
while (fgets(s, 256, f) != NULL){
printf("%d %s",bin,s);
bin++;
}
fclose(f);
return 0;
}
これをconv.cなどの名前で保存し、コンパイル、実行します。
> gcc -o conv conv.c
> chmod +x conv
> ./conv 137Cs_sample.spe > 137Cs_sample.dat
出力されるファイルの内容はAWKの場合と同じなので、Gnuplotでの表示方法も同じです。
gnupolot> plot '137Cs_sample.dat' using 1:2:3 with yerrorbar
Pythonで変換する
Pythonを使っても同じことができます。
import math
import sys
fname_in = sys.argv[1]
fname_out = sys.argv[1][:-4]+"+dat"
with open(fname_in) as fin:
lines = fin.read().split("\n")
with open(fname_out,mode="w") as fout:
for i in range(10,4106):
fout.write("{} {} {}\n".format(i-9,lines[i],math.sqrt(float(lines[i]))))
これをconv.pyなどの名前で保存し、実行します。ファイル名が同じで拡張子が .dat と言うファイルができます。
> python3 conv.py 137Cs_sample.spe
出力されるファイルの内容はAWKの場合と同じなので、Gnuplotでの表示方法も同じです。
gnupolot> plot '137Cs_sample.dat' using 1:2:3 with yerrorbar
Gnuplotでフィッティングをする
Gnuplotで任意関数でフィットをする方法を簡単に書いておきます。例として、前述の137Csの662 keVのピークをガウシアンでフィットしてみます。
まずはフィットする関数f(x)を定義します。
gnuplot> f(x) = a/(sqrt(2*3.14)*s) * exp(-(x-m)**2/(2*s**2))
パラメータは3つで、aがガウシアンの面積、sがsigma、mは中心値です。適当な初期値を入れてあげます。初期値があまりにずれているとうまくフィットできないので、初期値が正しいかまずはプロットして確認してみます。
gnuplot> a=300000
gnuplot> s=3
gnuplot> m=1565
gnuplot> plot [1520:1600] '137Cs_sample.dat' using 1:2:3 w yerr, f(x)
なんとなくそれっぽい初期値を見つけたら、あとはフィットします。
gnuplot> fit [1540:1580] f(x) '137Cs_sample.dat' u 1:2:3 via a,s,m
... (中略) ...
Final set of parameters Asymptotic Standard Error
======================= ==========================
a = 162194 +/- 8611 (5.309%)
s = 2.04345 +/- 0.1043 (5.103%)
m = 1565.24 +/- 0.1114 (0.007118%)
... (中略) ...
gnuplot> set samples 500
gnuplot> plot [1520:1600] '137Cs_sample.dat' using 1:2:3 w yerr, f(x)
関数が滑らかに見えるように set samples 500 としています。
実はgnuplotのfitで表示される "Asymptotic Standard Error" にはバグがあって、fitするデータに誤差がある場合に正しい値が計算されないことが知られています。詳細は P. Yang (2014) のAppendix D に書かれています。一部を引用すると、
to get correct error bars on fit parameters from gnuplot when there are error bars on the points, you have to divide gnuplot’s asymptotic standard errors by the square root of the chi-squared per degree of freedom (which gnuplot calls FIT STDFIT and, fortunately, computes correctly).
となります。gnuplotでは、fitをすると自動的に FIT_STDFIT という変数にこの "the square root of the chi-squared per degree of freedom" が入っているので、それで割れば良いです。表示の仕方は set fit errorvariables というコマンドでエラーを変数名_errに入れて、それを print というコマンドで表示されるのが楽です。
gnuplot> set fit errorvariables
gnuplot> fit .... (中略)
gnuplot> print a, a_err/FIT_STDFIT
これで、積分値が1.62 ± 0.09 x 10^5、sigmaが2.04 ± 0.10 、中心値は 1565.24 ± 0.11 という数値を得ることができました。下のような表示がされると思います。どうもこのデータはピークの形はガウシアンに比べて少しエネルギーの低い側にテールがあるみたいです。
Gnuplotのデフォルトの見栄えはあまりよくないので整形する
解析をしている時は気にしなくて良いですが、プレゼンで使用したりレポートや論文として使う場合にgnuplotの出力をそのまま使うのは可読性があまりよくないです。Gnuplot自体は投稿論文でも使われるソフトですが、そこそこ自分で設定しないと初期設定ではあまり見栄えがよくないのが難点です。
例えば上のグラフを少し見栄えよくしてみます。見栄えをよくしようとすると、流石にコマンドをいちいち打っていると大変なので、マクロファイルとして保存して、それをロードすることにします。
# =========== FIT ===========
set fit errorvariables
f(x) = a/(sqrt(2*3.14)*s) * exp(-(x-m)**2/(2*s**2))
a=300000
s=3
m=1565
# plot and check init parameters
plot [1520:1600] '137Cs_sample.dat' u 1:2:3 w yerr, f(x)
pause -1
# fit and plot
fit [1540:1580] f(x) '137Cs_sample.dat' u 1:2:3 via a,s,m
plot [1520:1600] '137Cs_sample.dat' u 1:2:3 w yerr, f(x)
pause -1
# =========== MAKE FIGURE BETTER ===========
# x,y label
set xlabel "ADC channel [ch]"
set ylabel "Counts"
set lmargin 15
set bmargin 5
set xlabel offset 0,-0.5 font "Arial,18"
set ylabel offset -3.5,0 font "Arial,18"
# tics
set rmargin 5
set tics font "Arial,14"
# legend
set key font "Ariel,18"
set key left top
# line styles
set bars fullwidth 0
set style line 1 lc "#0080FF" lw 1 pt 0 ps 0
set style line 2 lc "#0080FF" lw 2 pt 0 ps 0
set style line 3 lc "#FF3333" lw 1.5 pt 0 ps 0
# smooth function
set samples 500
# write fit results
set label 1 at first set label 1 at first 1530,20000 sprintf("Fit results:\nCounts= %8.0f +- %5.0f\nsigma= %4.2f +- %3.2f\nmean= %6.2f +- %3.2f",a,a_err/FIT_STDFIT,s,s_err/FIT_STDFIT,m,m_err/FIT_STDFIT) font "Arial,16"
# plot
plot [1520:1600] '137Cs_sample.dat' u 1:2:3 w yerr ls 2 title "Histogram",\
'' u ($1-0.5):2 w step ls 1 notitle,\
f(x) ls 3 title "Fit function (gaussian)"
これで下のようなグラフを作ることができます。ここまで見栄えにこだわる必要はないですが、人に見せる図には、最低限軸のラベルをつけて、線を太くしてみやすくする程度の気遣いは必要です。
文責：新倉潤 (niikura@phys.s.u-tokyo.ac.jp)