Python による放射線データの解析

tenoto/lecture_fitting @ GitHub

tags: lecture python　

高エネルギー天文学で使われる HEASoft や素粒子・原子核実験で使われる ROOT などの専用解析ツールを使わず、データ解析を行いたい場合もあります。Python だけで放射線データの基本的な解析をする手順を記載します。

練習: 乱数生成した標準正規分布をフィット

まず、実際の測定データを扱う前に、疑似データを乱数から生成して遊んでみましょう。

この一連のページでは、Python での解析と実例を紹介します。コードは Google Colaboratory に一部置いています。Google Colaboratory では標準で多くのライブラリをもっていますが、足りないものは !pip install xxx で使えるようにできます。

以下のすべてを Google Colab には置いていないので、注意してください。後半では、GitHub にコードを置いているで、Google Colab でも GitHub でも使いやすい方を試してください

乱数

乱数はいくつかのライブラリで扱えますが、様々な目的で使える numpy で呼ぶのが使いやすいことが多いです。

import numpy as np
np.random.seed(0)
evt = np.random.randn(1000)

ここで、seed(0) は乱数のシード設定、randn(1000) は標準正規分布(平均 0, 標準偏差 1)に従う乱数を 1000 個、生成して numpy の array に収容しています。この array は、素粒子物理学やＸ線天文学などで、放射線イベントの何らかの物理量の列に対応する、というイメージで evt という変数名にしてあります。

ヒストグラムに詰める

生成された数値データをヒストグラムに詰めてプロットするには、numpy.histogram が使えます。ヒストグラムの上限、下限とその間のビン数を指定するのは、ROOT (もっと古い古代語だと display45, pow?)などと同じにできます。

hy, xedges = np.histogram(evt,bins=40,range=(-4.,4.))
hx = 0.5*(xedges[1:]+xedges[:-1]) 
hyerr = np.sqrt(hy)

hy はヒストグラムに詰められた値で、xedges はヒストグラムの境界の値のため、ヒストグラムの中心値は hx = 0.5*(xedges[1:]+xedges[:-1]) で計算して hx に詰めておく必要があります。統計誤差は、各ヒストグラムの平方根で hyerr に詰めています。

(注意) 実は、天文学で扱うことがある小統計の統計誤差について、平方根を使うのには注意が必要です。Gehrels "N.Confidence Limits for Small Numbers of Events in Astrophysical Data", ApJ, 303, 336 (1986) や、HEASoft の mathpha のコード中での取り扱い Ian "The MATHPHA User’s Guide" (1995)　の "3 CALCULATION & PROGATION OF ERRORS" などが参考になります。

これらを用いて、

import matplotlib.pylab as plt 
fig = plt.figure(figsize=(6,4)) 
plt.errorbar(hx, hy, yerr=hyerr, marker='', drawstyle='steps-mid') 

とすれば、

のようにプロットできます。marker='', drawstyle='steps-mid' のオプションにすることで、標準的な物理での誤差付きヒストグラムの表現にしています。

フィッティング

上のヒストグラムは、標準正規分布から生成しているので、十分に統計を溜めれば、平均 0 で、標準偏差 1 のフィット結果に収束するはずです。

素粒子物理学で蓄積されてきたフィット（より正確には、評価関数の最小化）のアルゴリズムは、Minuit のライブラリになっており、python では、iminuit という名前で、一般的な関数最小化法として使用できます。これは、モデルのデータへの尤度フィットや、尤度プロファイル分析から、モデルのパラメータ誤差推定値を得るためにも使用できます。iminuit チュートリアルが参考になります。

さらに、probfit で iminuit で使いたい標準的な関数や、尤度関数、Chi-square 計算などを呼べるようにします。

Chi-square フィッティング

from iminuit import Minuit
from probfit import BinnedChi2, Extended, gaussian

ext_gauss = Extended(gaussian)
bc2 = BinnedChi2(ext_gauss, evt)

m = Minuit(bc2, mean=-1, sigma=0.5, N=100)

plt.figure(figsize=(12, 4))
plt.subplot(121)
bc2.draw(m)
plt.title('Before')

m.migrad() # fit

plt.subplot(122)
bc2.draw(m)
plt.title('After')

gaussian 関数を、Extended することで規格化因子を付けています。その上で、BinnedChi2 でガウス関数を生成したデータに適用したときの評価関数をつくり、Minuit で最小化できるように準備します。この際、初期値を入れています。m.migrad()がこの評価関数の最小値を探してくる本体です。

最小値を求めた後、以下のコマンドで非対称な誤差を計算します。

m.minos() 

これらの誤差は、以下で表示できます。

m.print_param()

フィットした後のパラメータや誤差を取り出すには、

m.values
m.errors
m.merrors

などとすればよいです。上記のフィットの結果から、中心値は

$-0.046\pm 0.032$ なので、

$2\sigma$の範囲で最初の乱数生成の仮定と一致し、幅については

$0.987\pm 0.026$ と

$1\sigma$の範囲で一致していることがわかります。

なお、今回は BinnedChi2 を使いましたが、他にも色々あります（詳細は別途）。

測定セットアップ

では、ここから実際に、測定データで試してみます。

• 放射線測定は、2020年11月9日、19日に、理化学研究所の榎戸研実験室(研究本館 610号室)で、茨城高専の村上悦基君、理研・榎戸研の沼澤正樹さんにやっていただきました。
• 放射線源として、ここでは 137Cs, 60Co の 2つの測定データと、放射線源を置かない場合の環境放射線バックグラウンドの測定データを扱う。これらの線源から放出される代表的なガンマ線輝線のエネルギーは物理定数表にも記載している。
• 検出器は「雷雲プロジェクト」用に開発してきた、コガモ検出器(Compact Gamma-ray Monitor; CoGaMo)を用いた。5x5x15 cm の CsI(Tl) 結晶シンチレータと MPPC 光検出器を用いている。

測定データと解析スクリプト

• 解析スクリプトは、tenoto/lecture_fitting@GitHub からダウンロードできる。
• 測定データは(今回の実習用に加工したものを) Dropbox からダウンロードでき、GitHub からプルした lecture_fitting の直下に置いてください。

準備

まず、必要なライブラリを読み込んでおきます。

#!/usr/bin/env python

import os 
import numpy as np
import pandas as pd
import matplotlib.pylab as plt 

from iminuit import Minuit
from probfit import BinnedChi2, Chi2Regression, Extended, gaussian

最後の iminuit と probfit が今回の放射線計測のために重要で、これらのライブラリは、pip install xxx (xxx はライブラリ名)などで準備してください。

CSV 形式のデータ読み込み

もし Google Colaboratory でデータアップロードする場合、

from google.colab import files
files = files.upload()

でアップロードできます。

data ディレクトリ下に置いてあるdata/201109_id30_137cs_15min.csv という CSV ファイルは、Time (sec) と ADC channel (pha) の列になっています。

python のライブラリ pandas には csv 形式を読める read_csv が用意されています。詳しくは、仕様を読んでもらえれば良いですが、コラム名はファイル内には含まれていないので、index_col=False, names=['sec','pha']で指定しています。

また、dtype でコラムのフォーマットを指定しています。sec の方は、誤ったフォーマットにすると数値がずれてしまうので注意。

読み込んだものは、data frame (df) に取り込まれます。

df_137cs = pd.read_csv('data/201109_id30_137cs_15min.csv',
	index_col=False, names=['sec','pha'],
	dtype={'sec':np.float,'pha':np.uint16})
print(df_137cs)

df_bgd = pd.read_csv('data/201109_id30_bkg_15min.csv',
	index_col=False, names=['sec','pha'],
	dtype={'sec':np.float,'pha':np.uint16})
print(df_bgd)

ライトカーブ（ヒストグラムを表示する）

matplotlib を使う方法

os.system('rm -rf out; mkdir -p out;')

nbins = 2**9 
tmax = 2**9 

fig = plt.figure(figsize=(8,5)) 
plt.hist(df_137cs['sec'],range=(0,tmax),bins=nbins,histtype='step')
plt.hist(df_bgd['sec'],range=(0,tmax),bins=nbins,histtype='step')
plt.xlim(0.0,tmax)
plt.xlabel('Time (sec)', fontsize=10)
plt.ylabel('Counts / 1-sec bin', fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/lc_137cs_bgd_1sec_hist.gif')

matplotlib の hist を使って、ヒストグラムを描画できます。ヒストグラムは上限、下限とビン数を指定し、表示方法を指定できます。今回、ビン数は 2^9 にしています。これは、ビンまとめするときに、2の乗数だと便利なことが多いからです。別に整数なら何でも構いません。

この図で、上の青は 137Cs、下のオレンジは環境放射線（何も放射線源を当てていない）の場合です。放射線を当てた方が、150 cps ほど増えていることがわかります。

numpy を使う方法

lc_src_y, edge = np.histogram(df_137cs['sec'],range=(0,tmax),bins=nbins)
lc_src_x = (edge[:-1] + edge[1:]) / 2.

lc_bgd_y, edge = np.histogram(df_bgd['sec'],range=(0,tmax),bins=nbins)
lc_bgd_x = (edge[:-1] + edge[1:]) / 2.

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(lc_src_x,lc_src_y,yerr=np.sqrt(lc_src_y),marker='',drawstyle='steps-mid')
plt.errorbar(lc_bgd_x,lc_bgd_y,yerr=np.sqrt(lc_bgd_y),marker='',drawstyle='steps-mid')
plt.xlim(0.0,tmax)
plt.xlabel('Time (sec)', fontsize=10)
plt.ylabel('Counts / 1-sec bin', fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/lc_137cs_bgd_1sec.png')

同じことを numpy でもできます。この際、yerr で誤差をつけてみました。

ビンの調整

今回の測定では、1秒でもそれなりの統計がありますが、例えば 8 秒でライトカーブが描きたい場合は以下のようになります。

tbin = 2**3 # 8 sec
nbins = int(nbins / tbin)

lc_src_y, edge = np.histogram(df_137cs['sec'],range=(0,tmax),bins=nbins)
lc_src_x = (edge[:-1] + edge[1:]) / 2.

lc_bgd_y, edge = np.histogram(df_bgd['sec'],range=(0,tmax),bins=nbins)
lc_bgd_x = (edge[:-1] + edge[1:]) / 2.

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(lc_src_x,lc_src_y,yerr=np.sqrt(lc_src_y),marker='',drawstyle='steps-mid')
plt.errorbar(lc_bgd_x,lc_bgd_y,yerr=np.sqrt(lc_bgd_y),marker='',drawstyle='steps-mid')
plt.xlim(0.0,tmax)
plt.xlabel('Time (sec)', fontsize=10)
plt.ylabel('Counts / (%d-sec bin)' % tbin, fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/lc_137cs_bgd_%dsec.png' % tbin)

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(lc_src_x,lc_src_y/float(tbin),yerr=np.sqrt(lc_src_y)/float(tbin),marker='',drawstyle='steps-mid')
plt.errorbar(lc_bgd_x,lc_bgd_y/float(tbin),yerr=np.sqrt(lc_bgd_y)/float(tbin),marker='',drawstyle='steps-mid')
plt.xlim(0.0,tmax)
plt.xlabel('Time (sec)', fontsize=10)
plt.ylabel('Counts / sec (%d-sec bin)' % tbin, fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/lc_137cs_bgd_%dsec_scale.png' % tbin)

最初のプロットでは、8 秒ビンごとの個数が表示されているので、1秒あたりに換算しています

ライトカーブを解析する

ライトカーブのヒストグラムを表示する

今回の測定は時間的に一定で、平均値は 426.5 counts/s です。np.mean(lc_src_y) とかやると計算できます。では、１秒値のライトカーブのレートがどの程度、分散があるのか調べてみます。これは、1秒値のレートの array をヒストグラムに詰めることになります。

lc_src_y, edge = np.histogram(df_137cs['sec'],range=(0,tmax),bins=tmax)

xmin = 300
xwid = 2**8
xmax = xmin + xwid 
nbins = 2**6

hist_y, hist_edge = np.histogram(lc_src_y,range=(xmin,xmax),bins=nbins)
hist_x = (hist_edge[:-1] + hist_edge[1:]) / 2.

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(hist_x,hist_y,yerr=np.sqrt(hist_y),marker='',drawstyle='steps-mid')
plt.xlim(xmin,xmax)
plt.xlabel('Rate (counts/sec)', fontsize=10)
plt.ylabel('Number of bins', fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/hist_lc_137cs_1sec.png')

となります。では、これをガウス関数でフィットしてみましょう。

def mygauss(x, mu, sigma, area):
	return area * np.exp(-0.5*(x-mu)**2/sigma**2)/(np.sqrt(2*np.pi)*sigma)

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(hist_x,hist_y,yerr=np.sqrt(hist_y),marker='',drawstyle='steps-mid')
plt.plot(hist_x,[mygauss(x,430,20,1000) for x in hist_x])
plt.xlim(xmin,xmax)
plt.xlabel('Rate (counts/sec)', fontsize=10)
plt.ylabel('Number of bins', fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/hist_lc_137cs_1sec_tmp1.png')

mygaus という関数を定義し、フィットの初期値を探すために、測定データを再現できるようなモデルを表示してみます。

chi2reg = Chi2Regression(mygauss,hist_x,hist_y)
fit = Minuit(chi2reg,mu=430,sigma=20,area=1000)
fit.migrad()
fit.minos() 

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(hist_x,hist_y,yerr=np.sqrt(hist_y),marker='',drawstyle='steps-mid')
model_mygauss = [mygauss(x,fit.values[0],fit.values[1],fit.values[2]) for x in hist_x]
plt.plot(hist_x,model_mygauss)
plt.xlim(xmin,xmax)
plt.xlabel('Rate (counts/sec)', fontsize=10)
plt.ylabel('Number of bins', fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/hist_lc_137cs_1sec_fit.png')

そこそこ良さそうな初期値を見つけたら、Chi-square で最適化してみます。

みたいになります。

print(fit.print_param())
print(fit.values)
print(fit.errors)

とすると、フィットで得られたパラメータとその誤差を表示できます。

┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│   │ Name  │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │
├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ mu    │  426.13   │   0.24    │   -0.24    │    0.24    │         │         │       │
│ 1 │ sigma │   20.76   │   0.23    │   -0.23    │    0.23    │         │         │       │
│ 2 │ area  │  2.073e3  │  0.021e3  │  -0.021e3  │  0.021e3   │         │         │       │
└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
None
<ValueView of Minuit at 7f92e8749060>
  mu: 426.134218801236
  sigma: 20.76215123588142
  area: 2073.1367177167863
<ErrorView of Minuit at 7f92e8749060>
  mu: 0.2432269180199041
  sigma: 0.23350886009026822
  area: 20.74372784277695

ガウス関数でフィットした中心値は 426.1+/-0.2 counts/sec ということで、最初に出した平均値 426.5 counts/s と

$2\sigma$の範囲で一致しています。

また、ガウス関数の幅 20+/-0.2 counts/sec ですが、sqrt(426)~20.6 counts/sec に（当たり前ですが）一致します。もし、この測定が統計的なバラツキ以外に変動をもつなら、レート平均値の平方根よりも幅が広がります. 今回の測定では、統計的な変動だけであることがわかります。

時間差の分布

時系列データはいろいろな解析ができます。連続して観測されるイベント間の到来時刻差(delta-t distribution)は、np.diff(df_137cs['sec']) で得ることができます。

これを片対数プロットに描くと、次のように直線になります。つまり、指数関数になります。137Cs の放射線の崩壊はランダムで、到来時間差は指数関数になります。

これは次のコードで描画し、指数関数でフィットできます。

xmin = 0.0005 
xmax = 0.02 + 0.0005 
nbins = 20

deltat_y, deltat_edge = np.histogram(np.diff(df_137cs['sec']),range=(xmin,xmax),bins=nbins)
deltat_x = (deltat_edge[:-1] + deltat_edge[1:]) / 2.

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(1000.0*deltat_x,deltat_y,
	yerr=np.sqrt(deltat_y),marker='',drawstyle='steps-mid')
#plt.xlim(0.0,900.0)
plt.xlabel('Delta time (msec)', fontsize=10)
plt.ylabel('Events', fontsize=10)
plt.xscale('linear')				
plt.yscale('log')
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")

def myexp(x, mu, norm):
	return norm * np.exp(-x/mu)

chi2reg = Chi2Regression(myexp,deltat_x,deltat_y)
fit = Minuit(chi2reg, mu=0.002356, norm=0.002356)
fit.migrad()
fit.minos() 
print(fit.print_param())
print(fit.values)
print(fit.errors)

model_myexp = [myexp(x,fit.values[0],fit.values[1]) for x in deltat_x]
plt.plot(1000*deltat_x,model_myexp,c='red')
plt.savefig('out/hist_deltat_137cs_fit.png')

さて、この指数関数のパラメータは、

┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│   │ Name │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │
├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ mu   │2.350325e-3│0.000034e-3│-0.000034e-3│0.000034e-3 │         │         │       │
│ 1 │ norm │164.6167e3 │ 0.0021e3  │ -0.0021e3  │  0.0021e3  │         │         │       │
└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
None
<ValueView of Minuit at 7fbab4b31780>
  mu: 0.0023503245242608357
  norm: 164616.71575401127
<ErrorView of Minuit at 7fbab4b31780>
  mu: 3.3962473970311133e-08
  norm: 2.1123118964314895
(py3.8.6env1) [enoto@chromy:l

ということで、

$\exp (-t/\mu )$ の

$\mu =0.00235$ ですが、

$1/\mu =425.5$ counts/sec とレートに対応することも確認できます。この delta-t distibution は本来ランダムな信号を検出器が正しく測定できているかなどの評価に使えます。

ちなみに、data frame を読み込むときの sec のフォーマットを間違うと、この delta-t ditribution は崩れることに注意してください。

スペクトルの解析

次に放射線のスペクトルを描いてみます。今回は、統計が少し悪く誤差が見えるデータの方が解説しやすいので、50 秒のデータでまずは描いてみましょう。

exposure = 50 # sec

flag_137cs = df_137cs['sec'] < exposure 
flag_bgd = df_bgd['sec'] < exposure 

spec_137cs_y, edge = np.histogram(df_137cs['pha'][flag_137cs], 
	bins=2**10, range=(-0.5, 2**10-0.5))
spec_137cs_x = (edge[:-1] + edge[1:]) / 2.
spec_137cs_x_err = (edge[:-1] - edge[1:]) / 2.
spec_137cs_y_err = np.sqrt(spec_137cs_y)

spec_bgd_y, edge = np.histogram(df_bgd['pha'][flag_bgd], 
	bins=2**10, range=(-0.5, 2**10-0.5))
spec_bgd_x = (edge[:-1] + edge[1:]) / 2.
spec_bgd_x_err = (edge[:-1] - edge[1:]) / 2.
spec_bgd_y_err = np.sqrt(spec_bgd_y)

fig, ax = plt.subplots(1,1, figsize=(11.69,8.27))
plt.errorbar(spec_137cs_x,spec_137cs_y,yerr=spec_137cs_y_err,
	marker='',drawstyle='steps-mid')
plt.errorbar(spec_bgd_x,spec_bgd_y,yerr=spec_bgd_y_err,
	marker='',drawstyle='steps-mid')
#plt.errorbar(x_bgd,y_bgd,yerr=np.sqrt(y_bgd),marker='',drawstyle='steps-mid')
plt.xlabel('ADC channel (pha)', fontsize=10)
plt.ylabel('Counts', fontsize=10)
#plt.xscale('log')				
plt.yscale('log')
plt.xlim(-0.5,300-0.5)
plt.tight_layout(pad=2)
#plt.tick_params(labelsize=fontsize)
#plt.rcParams["font.family"] = "serif"
#plt.rcParams["mathtext.fontset"] = "dejavuserif"	
plt.savefig('out/spec_137cs_bgd.png')

青が 137Cs を照射した場合で、オレンジ色が環境バックグラウンドの測定です。ADC channel で 60 以上の箇所は、バックグラウンドと137Cs は一致しています。

測定時間で割り算して（縦軸を counts から counts/sec　にする)、バックグラウンドを引いたデータについて、662 keV の輝線をガウス関数でフィットしてみます。

spec_137cs_sub_y = (spec_137cs_y - spec_bgd_y) / exposure 
spec_137cs_sub_yerr = np.sqrt(spec_137cs_y + spec_bgd_y) / exposure
spec_137cs_sub_x = spec_137cs_x
spec_137cs_sub_xerr = spec_137cs_x_err

chmin = 40
chmax = 2**6 

x = spec_137cs_sub_x[chmin:chmax]
y = spec_137cs_sub_y[chmin:chmax]
xe = spec_137cs_sub_xerr[chmin:chmax]
ye = spec_137cs_sub_yerr[chmin:chmax]

fig, ax = plt.subplots(1,1, figsize=(11.69,8.27))
plt.errorbar(x,y,yerr=ye,xerr=xe,marker='',drawstyle='steps-mid')
plt.xlabel('ADC channel (pha)', fontsize=10)
plt.ylabel('Counts', fontsize=10)
#plt.xscale('log')				
plt.yscale('log')
plt.xlim(-0.5,300-0.5)
#plt.tight_layout(pad=2)
#plt.tick_params(labelsize=fontsize)
#plt.rcParams["font.family"] = "serif"
#plt.rcParams["mathtext.fontset"] = "dejavuserif"	

def mygauss_linear(x, mu, sigma, area, c0=0.0, c1=0.0):
    return area * np.exp(-0.5*(x-mu)**2/sigma**2)/(np.sqrt(2*np.pi)*sigma) + c0 + c1 * x

chi2reg = Chi2Regression(mygauss_linear,x,y)
fit = Minuit(chi2reg, mu=50, sigma=2, area=300, c0=0.0, c1=0.0)
fit.migrad()
fit.minos() 
print(fit.print_param())
print(fit.values)
print(fit.errors)


fig = plt.figure(figsize=(7,6))
ax1 = fig.add_axes([0.1, 0.3, 0.85, 0.68])
ax1.errorbar(x,y,xerr=xe,yerr=ye,marker='o',ls='',color='k')
model_mygauss_linear = [mygauss_linear(i,
	fit.values[0],fit.values[1],fit.values[2],fit.values[3],fit.values[4]) for i in x]
#model_mygauss_linear = [mygauss_linear(i,
#	50,2,100,0,0) for i in x]
ax1.plot(x,model_mygauss_linear,c='red',drawstyle='steps-mid')
ax1.set_xlim(chmin,chmax)
#plt.xlabel('ADC channel', fontsize=10)
ax1.set_ylabel('Rate (counts/sec/bin)', fontsize=10)
ax1.axhline(y=0.0, color='k', linestyle='--')
ax1.get_xaxis().set_visible(False)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
# Second new axes

y_residual = y - model_mygauss_linear
ax2 = fig.add_axes([0.1, 0.1, 0.85, 0.20])
ax2.errorbar(x,y_residual,xerr=xe,yerr=ye,marker='o',ls='',color='k')
ax2.set_xlim(chmin,chmax)
ax2.axhline(y=0.0, color='r', linestyle='-')
ax2.set_ylim(-3.0,3.0)
ax2.set_xlabel('ADC channel', fontsize=10)
ax2.set_ylabel('Residual', fontsize=10)
#fig.align_ylabels()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/spec_137cs_bgd_fit.png')

今回のモデルでは、ガウス関数に加えて１次関数も加えています。図にはデータ点と、ベストフィットのモデル関数（赤）を示していて、下のパネルには残差を示しています。

最適のモデルパラメータは、

┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│   │ Name  │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │
├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ mu    │   50.46   │   0.08    │   -0.07    │    0.07    │         │         │       │
│ 1 │ sigma │   1.93    │   0.09    │   -0.08    │    0.08    │         │         │       │
│ 2 │ area  │    99     │     4     │     -4     │     4      │         │         │       │
│ 3 │ c0    │    4.5    │    1.6    │    -1.6    │    1.6     │         │         │       │
│ 4 │ c1    │  -0.075   │   0.031   │   -0.030   │   0.030    │         │         │       │
└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
None
<ValueView of Minuit at 7fe516f08b20>
  mu: 50.463643473409775
  sigma: 1.9267523240788242
  area: 99.07492309277377
  c0: 4.527929611161785
  c1: -0.07519915368978516
<ErrorView of Minuit at 7fe516f08b20>
  mu: 0.07590302508801514
  sigma: 0.08543496081129959
  area: 4.3795997490523275
  c0: 1.6199947988577938
  c1: 0.03050768910717568

ということで、ガウス関数のピークは 50.464+/-0.076 で、幅は Sigma=1.927+/-0.085、このピークのレートは 99.1+/-4.4 cps とわかります。

137Cs の 662 keV の輝線の時間変動

最初のライトカーブはエネルギーなどで選別はしていませんでしたが、137Cs の 662 keV をガウス関数で近似した場合のピークと幅が求まったので、3

$\sigma$ に入るイベントだけを抜き出した時間変動も描けるようになりました。

ch_662keV_mu = fit.values[0]
ch_662keV_sigma = fit.values[1]
print(ch_662keV_mu,ch_662keV_sigma)
chmin_662keV = ch_662keV_mu - 3.0 * ch_662keV_sigma
chmax_662keV = ch_662keV_mu + 3.0 * ch_662keV_sigma
print(chmin_662keV,chmax_662keV)

flag_137cs_662keV = np.logical_and(df_137cs["pha"] >= chmin_662keV, df_137cs["pha"] < chmax_662keV)
flag_bgd_662keV = np.logical_and(df_bgd["pha"] >= chmin_662keV, df_bgd["pha"] < chmax_662keV)

lc_137cs_y, edge = np.histogram(df_137cs['sec'][flag_137cs_662keV],range=(0,900.),bins=900)
lc_137cs_x = (edge[:-1] + edge[1:]) / 2.

bgd_lc_y, edge = np.histogram(df_bgd['sec'][flag_bgd_662keV],range=(0,900.),bins=900)
bgd_lc_x = (edge[:-1] + edge[1:]) / 2.

lc_137cs_mean = np.mean(lc_137cs_y)
lc_bgd_mean = np.mean(bgd_lc_y)

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(lc_137cs_x,lc_137cs_y-lc_bgd_mean,yerr=np.sqrt(lc_137cs_y),marker='',drawstyle='steps-mid')
#plt.errorbar(bgd_lc_time,bgd_lc_cnt,yerr=np.sqrt(bgd_lc_cnt),marker='',drawstyle='steps-mid')
plt.xlim(0.0,100.0)
plt.xlabel('Time (sec)', fontsize=10)
plt.ylabel('Counts / 1-sec bin', fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/lc_137cs_bgd_energy_cut.png')

このライトカーブでは、137Cs の 662 keV 輝線のレートから、同じエネルギー帯域（ADC channel)のバックグラウンドのレートの平均値を引いています。

• Average rate of the 662 keV line from 137 Cs = 125.6
• Average rate of the background = 18.1
• Subtracted average rate = 107.5

になります。ちなみに、スペクトルフィットでガウス関数と一次関数の和でモデル化したときのガウス関数のエリア（99 cps) と比べると 8 cps ほど多いのは、コンプトン等の成分が入っているから？？でしょうか。ここは未確認です。

60Co の 1173 keV と 1332 keV のライン

137Cs 同様に、60Co でのフィットを行うと以下のようになる。

今回は、２つのガウス関数と連続成分は２次関数を仮定している。

┌───┬────────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│   │ Name   │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │
├───┼────────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ mu1    │   89.90   │   0.05    │   -0.05    │    0.05    │         │         │       │
│ 1 │ sigma1 │   2.46    │   0.07    │   -0.07    │    0.07    │         │         │       │
│ 2 │ area1  │    230    │    10     │     -9     │     10     │         │         │       │
│ 3 │ mu2    │  102.27   │   0.06    │   -0.06    │    0.06    │         │         │       │
│ 4 │ sigma2 │   2.48    │   0.07    │   -0.07    │    0.07    │         │         │       │
│ 5 │ area2  │    205    │     9     │     -9     │     9      │         │         │       │
│ 6 │ c0     │    -30    │    26     │    -26     │     27     │         │         │       │
│ 7 │ c1     │    1.1    │    0.6    │    -0.6    │    0.6     │         │         │       │
│ 8 │ c2     │  -0.0073  │  0.0030   │  -0.0029   │   0.0030   │         │         │       │
└───┴────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
None
<ValueView of Minuit at 7facab4500c0>
  mu1: 89.89885681279526
  sigma1: 2.4635201398282667
  area1: 230.37813109658913
  mu2: 102.27387615149468
  sigma2: 2.476543896033184
  area2: 205.07381204912983
  c0: -29.742310299289525
  c1: 1.095319696384393
  c2: -0.007302664336358126
<ErrorView of Minuit at 7facab4500c0>
  mu1: 0.046191615919659215
  sigma1: 0.07260118129707424
  area1: 9.756022732605876
  mu2: 0.055135061931318875
  sigma2: 0.07349858005025595
  area2: 8.936257391186867
  c0: 26.114650409287663
  c1: 0.565827150139758
  c2: 0.0029774089088468775

ADC channel vs. Energy

さて、３本のラインでフィットができたので、Energy (keV) と ADC channel の相関を見てみる。

Ion	Energy (keV)	ADC channel	ADC error
137Cs	661.65	50.46	0.08
60Co	1173.23	89.90	0.05
60Co	1332.50	102.27	0.06

x = np.array([661.65,1173.23,1332.50],dtype='float') # energy 
y = np.array([50.46,89.90,102.27],dtype='float') # channel 
ye = np.array([0.08,0.05,0.06],dtype='float')

x2r = Chi2Regression(linear, x, y, ye)
fit = Minuit(x2r, m=1, c=2)
fit.migrad()
fit.minos() 
print(fit.print_param())
print(fit.values)
print(fit.errors)

fig = plt.figure(figsize=(8,5)) 
plt.errorbar(x,y,yerr=ye,marker='o',ls='',color='k')
model_x = range(0,15000)
model_linear = [linear(i,fit.values[0],fit.values[1]) for i in model_x]
plt.plot(model_x,model_linear,color='r',ls='--')
plt.xlim(0,15000)
plt.ylim(0,2**10)
plt.xlabel('Energy (keV)', fontsize=10)
plt.ylabel('ADC channel', fontsize=10)
plt.tight_layout()
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/energy_vs_channel.png')

┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│   │ Name │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │
├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ m    │ 77.20e-3  │  0.15e-3  │  -0.15e-3  │  0.15e-3   │         │         │       │
│ 1 │ c    │   -0.64   │   0.17    │   -0.17    │    0.17    │         │         │       │
└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
None
<ValueView of Minuit at 7fc4f78f78b0>
  m: 0.07720473893135693
  c: -0.6438362818461201
<ErrorView of Minuit at 7fc4f78f78b0>
  m: 0.00014678653797836546
  c: 0.1695069711635306

となるので、

$\frac{\text{ADC}}{\text{ch}}=0.0772\times \left(\frac{\text{Energy}}{\text{keV}}\right)-0.644$

という関係にあることがわかる。2^10=1024 channel は約 13 MeV に対応することがわかる。

環境バックグラウンド

最後に、ADC channel をエネルギーに変換できるようになったので、X軸をエネルギー、Y軸を Counts/bin ではなく、Counts/sec/keV の次元にそろえてみる。

exposure = 15 * 60  

hist_y, edge = np.histogram(df_bgd['pha'], bins=2**10,range=(-0.5, 2**10-0.5))
hist_x = (edge[:-1] + edge[1:]) / 2.
hist_xe = (edge[:-1] - edge[1:]) / 2.
hist_ye = np.sqrt(hist_y)

x = (hist_x+0.644)/0.0772/1000.0
xe = hist_xe/0.0772/1000.0
y = 1/77.2 * hist_y / exposure 
ye = 1/77.2 * hist_ye / exposure 

fig, ax = plt.subplots(1,1,figsize=(11.69,8.27))
plt.errorbar(x,y,xerr=xe,yerr=ye,marker='',drawstyle='steps-mid',color='k')
plt.xlabel('Energy (MeV)', fontsize=13)
plt.ylabel('Counts/sec/keV', fontsize=13)
plt.xscale('log')				
plt.yscale('log')
plt.xlim(0.3,13.0)
plt.tight_layout(pad=2)
plt.tick_params(labelsize=12)
plt.tick_params(axis="x",direction="in")
plt.tick_params(axis="y",direction="in")
plt.savefig('out/spec_bgd_energy.png')

などとなる。

ROOT のようにヒストグラムに詰める

今回は、順を追って理解できるように、クラス化したり、メソッド(関数)化したりはしていません。使ってくると、同じ操作はまとめたくなります。

たとえば、上記のヒストグラムに詰める操作は、以下のようなクラスを定義すると、ROOT っぽい操作をすることもできます。

class Hist1D(object):

  def __init__(self, nbins, xlow, xhigh):
    self.nbins = nbins
    self.xlow  = xlow
    self.xhigh = xhigh
    self.hist, edges = np.histogram([], bins=nbins, range=(xlow, xhigh))
    self.bins = (edges[:-1] + edges[1:]) / 2.

  def fill(self, arr):
    hist, edges = np.histogram(arr, bins=self.nbins, range=(self.xlow, self.xhigh))
    self.hist += hist

  @property
  def data(self):
    return self.bins, self.hist

を定義して、

h = Hist1D(220,0,12.0)
h.fill(hdu_clevt['EVENTS'].data['PI']/100.)
plt.step(*h.data)
plt.xlim(0,3.0)
plt.xlabel('Energy (keV)');plt.ylabel('Count/bin')
plt.show()

なお、ここでの du_clevt['EVENTS'].data['PI'] はイベントの array です。FITS 形式の説明のあとで使えるようになります。

参考文献

• probfit
• babar_python_tutorial/notebooks/04_Fitting.html
• Kai-Feng Chen, "Introduction to Numerical Analysis"
• http://www.sherrytowers.com/cowan_statistical_data_analysis.pdf
• Python modeling and fitting packages