Pico Tech • Peak Fit With Gnuplot Part6

Radiation Detection / Software Projects ... Peak Fit With R ... Peak Fit With Octave ... Peak Fit With Python / 日本語ページインデックス ... English Documents / 掲示板 / Peak Fit With Gnuplot

Peak Fit With Gnuplot Part6 色々な関数編

Peak Fit With Gnuplot 最初の一歩とＫ４０編、
Peak Fit With Gnuplot Part2 任意の単一ピーク編、
Peak Fit With Gnuplot Part3 複数ピーク＋核種データ編、
Peak Fit With Gnuplot Part4 検出器の特性調査編、
Peak Fit With Gnuplot Part5 検出器の自動特性調査編の続編です。

次は、
Peak Fit With Gnuplot Part7 実用的なスクリプトの実験。

Ｇｎｕｐｌｏｔの色々な使い方は、以下のページに書いてあります。
Gnuplot Memo
Gnuplot Memo3D
Gnuplot Memo2

私の目標は、まず、私が使っているＮａＩ３インチの測定器の応答関数をでっち上げることで、
それを使って、色々な核種のピークをフィットすることです。

で、その手始めに、Ｋ４０のピークを使って色々やったり、Ｒａ２２６やＥｕ１５２などで、
測定器の各種特性を調べたりしているわけですが、
肝心の応答関数の方は、コンプトンたんは、一応良さそうで、コンプトン丘も三次関数で良さそうなのですが、
４３０ｋｅＶ以下のエネルギーの低い方や、バックスキャッターの取り扱い、
そして、さらには２００ｋｅＶより下の方の扱いなどで悩んでいるところです。

で、悩んでいる理由は、応答関数を作るのに、どういう曲線を使うと
　　メインピークに連動させるのが楽で、
　　Ｇｎｕｐｌｏｔが楽に計算できて、
　　フィットが収束しやすいか、
という、基礎知識が圧倒的に欠けているので、
まず、色々な関数と曲線の例をこのところ（というか、かなり前から）眺めたり、
プロットしたり、フィットしたりして試しながら学習しているところです。

そういうわけで、このページは、そういう、色々な関数のメモのページです。

このページや、
http://www.lightstone.co.jp/origin/flist6.html

このページなどが、非常に参考になります。
https://www.hulinks.co.jp/software/peakfit/functions.html

ガウス関数のおさらい

ガウス関数と言っても、ピークの高さに着目した形、面積に着目した形、
統計の様に、面積が常に１になるように正規化したもの、などなどがあります。

で、私が今のところ良く使っているのは、高さに着目した形です。
これは、ｙ＝０の線の上の、エネルギーがＥの位置に、高さがｈで、幅（というか太り方）がｓ（シグマ）の
ガウス曲線を出す関数です。＊は、掛け算で、＊＊２はＧｎｕｐｌｏｔやＰｙｔｈｏｎでの二乗の書き方です。

g(x)=h*exp(-(x-E)**2/(2*s**2))

計算機に、もう少し優しい形にすると、こうでしょうか。
g1(x)=h*exp(-0.5*((x-E)/s)**2)

変数でなく、引数で、パラメーターを指定するなら、こうすればよいだけ。
gg(x,h,E,s)=h*exp(-0.5*((x-E)/s)**2)

高さを１、中心のエネルギーがゼロ、シグマが１でプロットすると、こうなります。
h=1;E=0;s=1
reset;unset xr;unset yr
set title 'h=1;E=0;s=1'
plot g1(x) t 'g1(x)=h*exp(-0.5*((x-E)/s))**2)'

ガウス関数では、ＦＷＨＭ（半値幅）が良く出てきて、それは、シグマと定数で簡単に出せます。
fwhm=s*2*sqrt(2*log(2))

この部分は、定数なので、一度計算して、変数に入れておくと、それ以降、Ｇｎｕｐｌｏｔは計算しないで済みます。
s2fwhm=2*sqrt(2*log(2))
fwhm=s*s2fwhm

このＦＷＨＭは、ｋｅＶの値で、よくある％の値にしたい時は、ピークのエネルギーＥで割って、１００をかけます。
Pfwhm=fwhm/E*100

また、面積（Ｎｐ）は、ｓ（シグマ）に√２πと高さをかければ良いので、簡単。
Np=h*sqrt(2*pi)*s

で、√２πも定数なので、こうしておくとＧｎｕｐｌｏｔは少し楽。
rt2pi=sqrt(2*pi)
Np=h*rt2pi*s

後は、中央から、両側にｓの分だけの面積は、全体の約６８％で、２ｓだと約９５％で、３ｓだと約９９．７％というのは、
統計で使うガウス曲線（正規分布）と同じなので、おぼえておくと、便利な場合があります。

このガウス曲線をｙ＝ｂ、傾きがａの直線の上に乗せたい時には、こんな風にするだけ。
要は、混ぜたい線を足してやれば良いだけです。
g2(x)=h*exp(-0.5*((x-E)/s)**2)+a*(x-E)+b
a=0.5;b=2
set title 'h=1;E=0;s=1;a=0.5;b=2'
plot g2(x) t 'g2(x)=h*exp(-0.5*((x-E)/s))**2)+a*(x-E)+b',g1(x),a*(x-E)+b

二次曲線や三次曲線に乗せても良いですし、まあ、なんでも好きに試してみると良いでしょう。

参考ページの式やグラフィック（これら二つは、ｙ０があるかどうかを除くと、同じものですが、変数の名前が違っています。）：
GaussAmp：ガウスピーク関数(振幅パラメータ) http://www.lightstone.co.jp/origin/flist6.html

Gaussian (Amplitude) https://www.hulinks.co.jp/software/peakfit/functions.html

エクストリーム関数（極値分布）

reset;unset xr;unset yr
set title 'h=1;E=0;w=1'
h=1;E=0;w=1
exf(x)=h*exp(-exp(-(x-E)/w)-(x-E)/w+1)
plot exf(x) t 'exf(x)=h*exp(-exp(-(x-E)/w)-(x-E)/w+1)'

引数指定形式
exfg(x,h,E,w)=h*exp(-exp(-(x-E)/w)-(x-E)/w+1)

バックグラウンドと一緒に出して、様子を見ています。
h=0.16
E=175
w=200
set title sprintf('h=%f, E=%d, w=%f',h,E,w)
plot [0:2600] exf(x), f every ::HL::BGL u 1:($2/sec) with dot

Extreme：統計分野で使われるExtreme関数 http://www.lightstone.co.jp/origin/flist6.html

Extreme Value (Amplitude)　https://www.hulinks.co.jp/software/peakfit/functions.html

ＬｏｇＮｏｒｍ

reset;unset xr;unset yr
set title 'h=1;E=100;w=1'
h=1;E=100;w=1
logn(x)=h*exp(-0.5*(log(x/E)/s)**2 )
plot [0:2000] logn(x) t 'logn(x)=h*exp(-0.5*(log(x/E)/s)**2 )'

引数指定形式
logng(x,h,E,s)=h*exp(-0.5*(log(x/E)/s)**2 )

h=0.17;E=160;s=0.85
set title 'h=0.17;E=160;s=0.85'
plot [0:2000] logn(x) t 'logn(x)=h*exp(-0.5*(log(x/E)/s)**2 )',f every ::HL::BGL u 1:($
2/sec) with dot

エクストリーム関数も（引数指定形式で）追加して、比較。
plot [0:2000] logn(x) t 'logn(x)=h*exp(-0.5*(log(x/E)/s)**2 )',f every ::HL::BGL u 1:($
2/sec) with dot,exfg(x,0.17,160,120)

二つを足して２で割って、平均化
plot [0:2000] logn(x) t 'logn(x)=h*exp(-0.5*(log(x/E)/s)**2 )',f every ::HL::BGL u 1:($
2/sec) with dot,exfg(x,0.17,160,120),0.5*(exfg(x,0.17,160,120)+logn(x))