Peak Fit With Gnuplot Part7 実用的なスクリプト
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 Part6 色々な関数編。
Ｇｎｕｐｌｏｔの色々な使い方は、以下のページに書いてあります。
Gnuplot Memo
Gnuplot Memo3D
Gnuplot Memo2
実用的なスクリプト
応答関数をでっち上げて、色々やろうと思ったのですが、
その前にもう少し簡単なやり方で、実用的なスクリプトを作ってみることにしました。
ここで最終的に目指すのは、汎用のスクリプトで、５００ｋｅＶくらいから１１００ｋｅＶくらいの範囲で
ベースラインを三次式で近似し、その上にＣｓ１３７，Ｃｓ１３４，ウラン系、トリウム系などの光電ピークのみをフィットし、
Ｋ４０については、その後で他の核種の線量を固定して、Ｋ４０の値だけフィットする、という方法でやってみます。
これは、Ｋ４０をまずフィットして、その値でベースラインを推測しようとすると、ウラン系やトリウム系の影響、
そして、Ｃｓ１３４やそのサムピークの影響を受けたり、応答関数の精度によってセシウム汚染の推定への
影響が大き過ぎたりするためです。
使用する式のおさらい
# 放射能（Ｂｑ）とピークの計数率（ｃｐｓ）の関係。
# 実測値から割り出すピーク効率の式に、補正などの分を含んでいる為、簡略化されてます。
# 放射能（Ｂｑ） = 計数率（ｃｐｓ） / （ピーク効率 ｘ 放出率）
# A = n / (ef * r)
# n = A * ef * r
#
# 「ピーク効率」は、ピークのエネルギー位置と効率曲線の式から求めます。ef(E)
# 「放出率」と「ピークのエネルギー」は、核種データ（nlist.plt）から読み込みます。核種データは、現在は仮の形式および数値です。
# 「効率曲線」は、自動特性調査スクリプト（ft3.plt）で割り出すか、他の方法で推測します。 Peak Fit With Gnuplot Part5
# 「計数率」は、ピークと思われるデータに、ガウス曲線をフィットして求めます。
#
# ピークフィットに使用するガウス曲線の式
# ｈ：ピーク高、Ｅ：ピークのエネルギー位置、ｓ：ガウス曲線の幅（シグマ）で分解能に比例、ａ，ｂ：ベースラインを直線とした場合のパラメータ
# g(x)=h*exp(-((x-E)/s)**2/2)+a*(x-E)+b
#
# ピーク高に、ルート２パイ＝√(２π) とシグマ(s)をかけると、計数率(n)が出ます。
# rt2pi=sqrt(2 * pi)
# n = h * rt2pi * s # 計数率ｎ（ｃｐｓ） ＝ ピーク高ｈ（ｃｐｓ） ｘ √(２π) ｘ シグマｓ
# h = n /(rt2pi * s) # ピーク高ｈ（ｃｐｓ） = 計数率ｎ（ｃｐｓ） / （√(２π) ｘ シグマｓ）
# h = (A * ef * r) / (rt2pi * s)
#
#
# この他に、特定のエネルギーの分解能を得る、分解能曲線の式、エネルギーからチャンネルを求める校正式、
# などを使います。これらは、自動特性調査スクリプトで推測、生成されるか、適当に設定する必要があります。
ピークフィットの実験が少し前進
まだ、不確かさはＧｎｕｐｌｏｔのフィットのエラーの数字をそのまま出しているだけで、
検出限界とかも計算していません。
効率曲線もしっかり確認してないので、その辺はテキトーです。
（おそらくこのせいもあって、Ｃｓ１３７，Ｃｓ１３４の比が正確ではありません）
ガウス曲線の式に色々詰め込んで、フィットすると検体中のそれぞれの核種の放射能がそのまま出ます。
濃度にしたい場合は、それを重量で割るだけ。
なので、ＧｎｕｐｌｏｔのＦｉｔ出力を見ると、それで直ぐに検体中の放射能（Ｂｑ）が分かります。
計数率を知りたい場合は、逆に計算をする必要がありますが、それも帳票形式の出力を
自動でする様にしてあります。（まだ、追加してない項目などがありますが）
凡例に、核種のエネルギー、放出率、効率、ピーク面積（計数率）、放射能と不確かさ、放射能濃度（Ｂｑ/ｋｇ）を
出すようにしてみました。
Ｃｓ１３７が１０００Ｂｑ/ｋｇ程度の土を小分けにして２０ｇくらいにしたもの。
Ｃｓ１３７が１０００Ｂｑ/ｋｇ程度の土を小分けにして１０ｇくらいにしたもの。
線香の灰（原発大災害以前のもの） ６３ｇ （重量を６２ｇと入れてました・・・）
ＫＣＬ ５０ｇ
粒状炭（原発大災害以前のもの） ３２７ｇ これは、５時間程度の測定で、
もっと長時間のデータもあったと思いますが、短時間だとどうなるかのテストで
まだ、検出限界とかを計算するようにしていないので、正確にはわかりませんが、
否検出になるかもしれないケース。
Ｇｎｕｐｌｏｔ画面の出力の例：
最初の段階のＫ４０の値は、出鱈目です。
次のＫ４０専用のフィットで値が出て、ＢＧのＫ４０の値を引いて実際の値（と思われるもの）を計算しています。
Final set of parameters Asymptotic Standard Error ======================= ========================== bb0 = 0.0368674 +/- 0.007157 (19.41%) bb1 = -7.29122e-005 +/- 3.042e-005 (41.72%) bb2 = 6.42779e-008 +/- 4.111e-008 (63.96%) bb3 = -2.08948e-011 +/- 1.785e-011 (85.43%) BqCs137 = 0.327326 +/- 0.1337 (40.84%) BqCs134 = 0.0952125 +/- 0.1284 (134.9%) BqAnn = 0.642645 +/- 0.06478 (10.08%) BqK40 = 267.871 +/- 124.9 (46.61%) BqRa226 = 0.442247 +/- 0.2872 (64.95%) BqTh232 = 0.249158 +/- 0.3075 (123.4%)
Final set of parameters Asymptotic Standard Error ======================= ========================== bK40a = -5.59262e-006 +/- 4.912e-007 (8.784%) bK40b = 0.0110037 +/- 0.00072 (6.544%) BqK40 = 111.169 +/- 1.739 (1.564%)
ftitle = t 'ThereminoMCA_2015_03_22_05_49_19-RyujouTan-327g3_3-8_TPA_LC4-N3.txt 18719 sec' Mass (weight) = 0.327kg ======================================================================================= Isotope Energy Ratio Efficiency Rate Activity Uncertainty Act(Bq/kg) --------------------------------------------------------------------------------------- Cs137 661.7keV, 85.1%, 0.073cps/Bq, 0.02023cps 0.3Bq +-0.13Bq 1.0Bq/kg Cs134 795.9keV, 85.5%, 0.059cps/Bq, 0.00478cps 0.1Bq +-0.13Bq 0.3Bq/kg K40 1460.8keV, 10.66%, 0.032cps/Bq, 0.37344cps 111.2Bq +-1.74Bq 340.0Bq/kg Ra226 609.3keV, 46.09%, 0.080cps/Bq, 0.01636cps 0.4Bq +-0.29Bq 1.4Bq/kg Th232 238.6keV, 43.65%, 0.379cps/Bq, 0.04127cps 0.2Bq +-0.31Bq 0.8Bq/kg ======================================================================================= K40 Activity = 12.838331838516Bq/kg Cs137:Cs134 = 1:0.29
スクリプトは、もう少し色々と手を加えてからアップロードします。
まだ、実験中のゴミもあるし、コメントとかも書いてない部分があったり、
解決したい課題も残っていますので。
どうしても直ぐに見たい、という方は、nkom あっと pico.dreamhosters.com までご連絡ください。
工事中