X線自己吸収について
X線による自己吸収、という言葉は分野によって定義が様々であるが、ここでは次の図に示したような、単純な一層の媒質にX線が照射した場合に、入射後に光電吸収して出射する単純な場合を考えます。
下記の記事の内容は、全て、google colab 上で動作するように、https://colab.research.google.com/drive/1tHsspwdMx_4vdbtKZ_gEeSDg89R7w7XX?usp=sharing を作成してますので、コードは この記事のcolabの を参照ください。
xraylib を使うことが前提なので、xraylib の動作方法などは、
などを参考にしてください。
自己吸収の式の導出
細かい定義などは、Broadband high-energy resolution hard x-ray spectroscopy using transition edge sensors at SPring-8, Review of Scientific Instruments 92, 013103 (2021); https://doi.org/10.1063/5.0020642
の式(6)(7)の notation に従うことにします。
下記では、透過する距離は、x/sin(θ)、x/sin(φ)、のようにいくらか補正がかかるが、一般的には45度程度で入射と出射すると考えて、最初は簡単のために角度は露わには書かないことにします。
ρ は全体の物質の密度で、ρi は i番元素のみの密度、Ein は、入射X線のエネルギー、Ea は出射X線のエネルギーです。σatt は、全質量吸収係数で、σia は、i番元素のaという輝線の発光断面積です。
\begin{align}
I_i^{\alpha}(E_{\text {in }}, E_{\alpha})/I_0 &= \int_0^d \left\{e^{-\sigma^{\text {att }}\left(E_{\text {in }}\right) \rho x}\right\} \left\{1 - e^{-\sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i dx} \right\} \left\{e^{-\sigma^{\text {att }}\left(E_{\text {in }}\right) \rho dx} \right\} \left\{e^{-\sigma^{\text {att }}\left(E_{\alpha}\right) \rho x}\right\} \\
&= \int_0^d \left\{入射光子透過率 \right\} \left\{i番元素吸収率 \right\} \left\{微小透過率(ほぼ1) \right\} \left\{出射光子透過率 \right\} \\
&= \int_0^d \left\{e^{-\sigma^{\text {att }}\left(E_{\text {in }}\right) \rho x}\right\} \left\{\sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i\right\} \left\{e^{-\sigma^{\text {att }}\left(E_{\alpha}\right) \rho x}\right\} dx \\
&= \sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i d \times \frac{
\left[ 1 - e^{- ( \sigma^{\text {att }}\left(E_{\text {in}}\right) + \sigma^{\text {att }}\left(E_{\alpha}\right)) \rho d) } \right]}{(\sigma^{\text {att }}\left(E_{\text {in}}\right) + \sigma^{\text {att }}\left(E_{\alpha}\right)) \rho d} \\
&= (i番目元素による吸収) \times (自己吸収)
\end{align}
上式を全部まとめて自己吸収補正という時があるが、その言葉は必ずしも最適な表現ではなく、自己吸収による効果と、該当元素による効果の積に分けた方が考えやすい、というのがこの記事の一番のポイントです。
入射X線の吸収と、出射X線の吸収を合わせて一つの変数 $x$ (平均的な optical depth に相当する量)として、
x \equiv (\sigma^{\text {att }}\left(E_{\text {in}}\right) + \sigma^{\text {att }}\left(E_{\alpha}\right)) \rho d
と定義し、
f(x) \equiv \frac{1-e^{-x}}{x}
という関数を用いて書き直すと、
\begin{align}
I_i^{\alpha}(E_{\text {in }}, E_{\alpha})/I_0 &= \sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i d \times \frac{
\left[ 1 - e^{- ( \sigma^{\text {att }}\left(E_{\text {in}}\right) + \sigma^{\text {att }}\left(E_{\alpha}\right)) \rho d) } \right]}{(\sigma^{\text {att }}\left(E_{\text {in}}\right) + \sigma^{\text {att }}\left(E_{\alpha}\right)) \rho d} \\
&= \sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i d \times f(x) \\
\end{align}
となる。この x は、X線の入射エネルギー Ein と、X線の出射エネルギー Ea の２つの変数に依存するのですが、光電吸収が卓越する帯域では、E^-3.5 の依存性があるので、Ea < Ein の場合は σ(Ea) >> σ(Ein) なので、Ein にはあまり依存しないです。
f(x) のグラフの概形は、
#!/usr/env/python
import numpy as np
import matplotlib.pyplot as plt
params = {'xtick.labelsize': 13, # x ticks
'ytick.labelsize': 13, # y ticks
'legend.fontsize': 13, # legend
}
plt.rcParams['font.family'] = 'serif'
plt.rcParams["font.size"] = 16
plt.rcParams.update(params)
x = np.linspace(0,100,100000)
y = (1 - np.exp(-x))/x
F = plt.figure(figsize=(12,8.))
ax = plt.subplot(1,1,1)
plt.xscale('linear')
plt.yscale('linear')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'X')
plt.ylabel(r'Y')
plt.plot(x,y)
outfile = "plot_func.png"
plt.savefig(outfile)
などでお絵かきをして確認すると、
となるので、x が 0 で 1 に収束し、xが大きくなると 0 に収束する。密度ρは、現実的には有限なことがあるので、厚み d が十分に小さいと f(x)はほぼ1となり、自己吸収が効かない、ということになる。
入射角 θ と出射角度 φ を省略せずに書いてみると次のようになります。
\begin{align}
I_i^{\alpha}(E_{\text {in }}, E_{\alpha})/I_0 &= \int_0^d \left\{e^{-\sigma^{\text {att }}\left(E_{\text {in }}\right) \rho \frac{x}{\sin \theta}}\right\} \left\{1 - e^{-\sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i dx} \right\} \left\{e^{-\sigma^{\text {att }}\left(E_{\text {in }}\right) \rho dx} \right\} \left\{e^{-\sigma^{\text {att }}\left(E_{\alpha}\right) \rho \frac{x}{\sin \phi}}\right\} \\
&= \int_0^d \left\{入射光子透過率 \right\} \left\{i番元素吸収率 \right\} \left\{微小透過率(ほぼ1) \right\} \left\{出射光子透過率 \right\} \\
&= \int_0^d \left\{e^{-\sigma^{\text {att }}\left(E_{\text {in }}\right) \rho \frac{x}{\sin \theta}}\right\} \left\{\sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i\right\} \left\{e^{-\sigma^{\text {att }}\left(E_{\alpha}\right) \rho \frac{x}{\sin \phi}}\right\} dx \\
&= \sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i d \times \frac{
\left[ 1 - e^{- ( \sigma^{\text {att }}\left(E_{\text {in}}\right) \frac{d}{\sin \theta} + \sigma^{\text {att }}\left(E_{\alpha}\right) \frac{d}{\sin \phi} ) \rho) } \right]}{(\sigma^{\text {att }}\left(E_{\text {in}}\right) \frac{d}{\sin \theta} + \sigma^{\text {att }}\left(E_{\alpha}\right) \frac{d}{\sin \phi} ) \rho} \\
&= (i番目元素による吸収) \times (自己吸収)
\end{align}
ここで、先ほどの同様の入射X線の吸収と、出射X線の吸収を合わせて一つの変数 $x$ (平均的な optical depth に相当する量)を
x \equiv (\sigma^{\text {att }}\left(E_{\text {in}}\right) \frac{d}{\sin \theta} + \sigma^{\text {att }}\left(E_{\alpha}\right) \frac{d}{\sin \phi} ) \rho
と定義する。上で定義した f(x)を用いると、
\begin{align}
I_i^{\alpha}(E_{\text {in }}, E_{\alpha})/I_0
&= \sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right) \rho_i d \times f(x) \\
&= (i番目元素による吸収) \times (自己吸収)
\end{align}
となり、自己吸収がない場合の発光量と自己吸収を分けることができ、入射角 θ と出射角度 φ の依存性は x に押し込めることができる。
ここで、$\sigma_{\mathrm{i}}^\alpha\left(E_{\text {in }}\right)$ は、「蛍光X線生成断面積 (flourescence cross section)」 という定義で、
に詳しく書いたのでそちらを参考にしつつ、具体的には、
xraylib.CS_FluorLine_Kissel_Cascade(i, xraylib.KL3_LINE, inenergy)
のように、原子番号(i)、輝線の名前(xraylib.KL3_LINE)、入射X線エネルギー (inenergy) keV、の３つのパラメータを与えると、蛍光X線生成断面積を g/cm^2 の単位で返してくれる。
自己吸収の計算方法と具体例
ここでは、「太陽組成のリュウグウ試料の自己吸収」「太陽組成の降着円盤の自己吸収」「NIST610の自己吸収」のプロットの例を紹介します。データは私の作成したcsvファイルをダウンロードして使う仕様になってます。全て、google colab 上で動作するように、
https://colab.research.google.com/drive/1tHsspwdMx_4vdbtKZ_gEeSDg89R7w7XX?usp=sharing
を作成してますので、コードはこちらを参照ください。
太陽組成のリュウグウ試料の自己吸収
太陽組成については、
を参考ください。ここでは、Lodders (2003) の abundance の物質が、密度は 1.27 g/cm^3 で存在した時に、厚みが 10um, 100um, 1000um の３パターンの場合に、自己吸収 f(x) および x がどうなるかを計算した例です。
(入射エネルギーは 18 keV を勝仮定して計算してますが、30keVでも50keVでも、出射するエネルギー(の方が普通は小さい）の方で x が決まるので、入射エネルギーにはあまり依存しません。)
コードの一部はこのような感じで、全体は https://colab.research.google.com/drive/1tHsspwdMx_4vdbtKZ_gEeSDg89R7w7XX?usp=sharing をみてください。
def CalcSelfAbs_Lo2003_Ryugu(xmin = 1e2, xmax = 18e3, xbins = 1e3, theta = np.pi/4, phi = np.pi/4, pltflag = True):
# e.g., Ryugu : density = 1.27 g/cm^3, thickness = 1um
density = 1.27 # 1.27 g/cm^3
ev2kev = 1e-3
list_ene_ev = np.linspace(xmin,xmax,int(xbins))
input_ene_ev = list_ene_ev[-1]
reduce_flag = True
thickness = np.array([1e-3, 1e-2, 1e-1]) # 10um, 100um, 1000um [cm]
Avogadro_constant = 6.02e23
H_1mol = 1.0 # g/mol
number_density = Avogadro_constant * (density / H_1mol) # 1/cm^3
nh = number_density * thickness # [1/cm^2]
sfx = []
slist_x = []
for one_t,one_nh in zip(thickness,nh):
print("density = ", density, " [g/cm^3] thickness = ", one_t, " [cm] Nh = ", one_nh, " [1/cm^2]")
list_x, fx = _Get_x_fx(Get_CS_Total_Lo2003, input_ene_ev, list_ene_ev, density, one_t, reduce_flag, theta, phi)
sfx.append(fx)
slist_x.append(list_x)
# plot
F = plt.figure(figsize=(8,10))
ax = plt.subplot(2,1,1)
plt.figtext(0.1,0.96,"Ryugu using Lodders (2003) : density = " + str("%3.2e" % density) + r" [g/cm$^3$]")
plt.xscale('log')
plt.yscale('log')
plt.grid(linestyle='dotted',alpha=0.5)
plt.ylabel(r'Self Absorption : f(x)')
for one_fx, one_x, one_nh, one_t in zip(sfx, slist_x, nh, thickness):
plt.errorbar(list_ene_ev, one_fx, label = "Ein = " + str("%3.1f" % (input_ene_ev * ev2kev)) + " (keV) Nh = " + str("%3.1e" % one_nh) + r" [cm$^{-2}$] " + " d = " + str("%3.3f" % one_t) + r" [cm]")
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=1, borderaxespad=0.,fontsize=7)
ax = plt.subplot(2,1,2)
plt.xscale('log')
plt.yscale('log')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (eV)')
plt.ylabel(r'mean optical depth : x')
for one_fx, one_x, one_nh in zip(sfx, slist_x, nh):
plt.errorbar(list_ene_ev, one_x, label = "Ein = " + str("%3.1f" % (input_ene_ev * ev2kev)) + " (keV) Nh = " + str("%3.1e" % one_nh) + r" [cm$^{-2}$]")
outfile = "selfabs_ryugu_ein" + str("%3d" % input_ene_ev) + "eV.png"
plt.savefig(outfile)
print("saved as ", outfile)
if pltflag: plt.show()
太陽組成の降着円盤の自己吸収
同様に、
- density = 1.66-07 [g/cm^3] thickness = 100.0 [m] Nh = 1e+21 [1/cm^2]
- density = 1.66e-07 [g/cm^3] thickness = 1000.0 [m] Nh = 1e+22 [1/cm^2]
- density = 1.66e-07 [g/cm^3] thickness = 10000.0 [m] Nh = 1e+23 [1/cm^2]
というパラメータを持つ降着円盤の場合について計算した結果がこちらです。
地上にある物質はオーダーで密度が 1g/cc 程度で、それが数100ミクロン程度で、柱密度が 10^22 [1/cm^2] 程度になります。ブラックホールの降着円盤の場合は、密度がその７桁落ちくらいで 1e-7g/cc の場合で、厚みが７桁分厚い 1000 m と同じ程度の柱密度となって、X線の見え方としては同じようになります。
(こう考えると、地球上の砂つぶの分析と、ブラックホールの降着円盤のX線観測は、一見違うようで、似たようなことに思えてきますね。)
NIST610 の自己吸収
最後は、ちょっとアバンダンスが太陽組成ではなくて、重元素が多く含まれた標準試料(NIST610)について、比較してみましょう。
厚みはだいたい1mmの試料ですが、ここでは思考実験として、薄い場合の計算をしています。
- density = 2.65 [g/cm^3] thickness = 0.001 [cm] Nh = 3.7e+19 [1/cm^2]
- density = 2.65 [g/cm^3] thickness = 0.01 [cm] Nh = 3.7e+20 [1/cm^2]
- density = 2.65 [g/cm^3] thickness = 0.1 [cm] Nh = 3.7e+21 [1/cm^2]
柱密度 Nh の計算は、単純に密度を平均分子量で割っただけなので、雑な計算ですが、重元素が多い分だけ、太陽組成に比べて小さい柱密度でも結構吸収が効いてることがわかります。
蛍光X線の比の活用
このように、ある媒質にX線が照射されて、そこからの発光を見る、ということは、その媒質の中での吸収という情報がエンコードされていることに相当しています。それはつまり、観測された発光の情報から媒質の吸収という情報をデコードすることも原理的には可能なはずです。
鉄の Ka と Kb の活用
具体的に鉄のKaとKbについて考えてみます。
x(E) \equiv (\sigma^{\text {att }}\left(E_{\text {in}}\right) \frac{d}{\sin \theta} + \sigma^{\text {att }}\left(E \right) \frac{d}{\sin \phi} ) \rho
とすると、
\begin{align}
I_{Fe}^{K \alpha}(E_{\text {in }}, E_{Fe K\alpha})/I_0
&= \sigma_{\mathrm{Fe}}^{K\alpha}\left(E_{\text {in }}\right) \rho_{Fe} d \times f(x(E_{Fe K\alpha})) \\
\end{align}
\begin{align}
I_{Fe}^{K \beta}(E_{\text {in }}, E_{Fe K\beta})/I_0
&= \sigma_{\mathrm{Fe}}^{K \beta}\left(E_{\text {in }}\right) \rho_{Fe} d \times f(x(E_{Fe K \beta})) \\
\end{align}
となりこれらの比をとると、
\begin{align}
\frac{I_{Fe}^{K \beta}(E_{\text {in }}, E_{Fe K\beta})}{I_{Fe}^{K \alpha}(E_{\text {in }}, E_{Fe K\alpha})} &= \frac{\sigma_{\mathrm{Fe}}^{K \beta}\left(E_{\text {in }}\right) \times f(x(E_{Fe K \beta}))}{\sigma_{\mathrm{Fe}}^{K \alpha}\left(E_{\text {in }}\right) \times f(x(E_{Fe K \alpha}))} \\
&\simeq 0.1 \times \frac{f(x(E_{Fe K \beta}))}{f(x(E_{Fe K \alpha}))}
\end{align}
最後は、KaとKbの強度比がだいたい10:1であることを使ったが、量子力学的な遷移確率で決まるものと思えば、強度比が理論的に決まることがポイントで、その場合は、実測された輝線の強度比から、f(x)の２点が決まる。f(x)には、厚みd, 入射角度θ、出射角度φ、密度ρ、などの未知数があるが、dρ や、φなど、部分的に制約することができる。FeのKa, Kb以外に、Cr, Ni, などのKa,Kbも使うことができるとより制約できるはず。
その他、FeのLラインとKラインの比や、電離した輝線など、理論的にある程度の精度で分岐比がわかっているラインの比を使うと、吸収量の診断ができはずです。
自己吸収の有無による蛍光X線の違い
「蛍光X線生成断面積 (flourescence cross section)」 の計算方法
Lnames という配列に、KとLラインの名前のリストを作って、xraylib_lname = eval("xraylib."+lname) という eval を用いて、文字列から xraylib.ライン名 の変数を取得し、 cs = xraylib.CS_FluorLine_Kissel_Cascade(an, xraylib_lname, iene_kev) * density * length の部分で、 「蛍光X線生成断面積 (flourescence cross section)」 を計算している。
Lnames = \
["KL2_LINE", # Ka2
"KL3_LINE", # Ka1
"KM2_LINE", # Kb3
"KM3_LINE", # Kb1
"L3M4_LINE", # La2
"L3M5_LINE", # La1
"L2M4_LINE", # Lb1
"L2N4_LINE", # Lg1
"L1N2_LINE", # Lg2
"L1N3_LINE"] # Lg3
def Get_LineEnergy_CSF(an, iene_kev = 20.0, density = 1.0, length = 1.0, debug = False):
list_lene_kev = []
list_csf = []
list_lname = []
elem = xraylib.AtomicNumberToSymbol(an)
for lname in Lnames:
# print(lname)
xraylib_lname = eval("xraylib."+lname)
try:
le = xraylib.LineEnergy(an, xraylib_lname) # energy
# CS_FlourLine [cm^2/g] * density [g/cm^3] * length [cm]
cs = xraylib.CS_FluorLine_Kissel_Cascade(an, xraylib_lname, iene_kev) * density * length
list_lene_kev.append(le)
list_csf.append(cs)
list_lname.append(elem + " " + lname)
except:
if debug: print("not found.", an, lname)
continue
list_lene_kev = np.array(list_lene_kev)
list_csf = np.array(list_csf)
if not len(list_lene_kev) == len(list_csf):
print("[ERROR] length mismatch")
return -1
return list_lene_kev, list_csf, list_lname
Get_LineEnergy_CSF(26, iene_kev = 20.0, density = 1.0, length = 1.0)
これを用いて、計算してプロットすると、
import matplotlib.colors as colors
import matplotlib.cm as cm
def Plot_CS_FluorLine(list_an, iene_kev = 20.0, pltflag = True):
list_an = np.array(list_an)
F = plt.figure(figsize=(12,8.))
ax = plt.subplot(1,1,1)
plt.xscale('linear')
plt.yscale('log')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency')
usercmap = plt.get_cmap('jet')
cNorm = colors.Normalize(vmin=0, vmax=len(list_an))
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=usercmap)
for i, an in enumerate(list_an):
elem = xraylib.AtomicNumberToSymbol(an)
c = scalarMap.to_rgba(i)
list_lene_kev, list_csf, list_lname = Get_LineEnergy_CSF(an, iene_kev)
for j, (lene, csf, lname) in enumerate(zip(list_lene_kev, list_csf, list_lname)):
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c, label = elem)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=10, borderaxespad=0.,fontsize=9)
else:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c)
outfile = "fig_cs_flour.png"
plt.savefig(outfile)
このように、各元素のラインに色をつけて表示される。縦軸は、xraylib.CS_FluorLine_Kissel_Cascade(an, xraylib_lname, iene_kev) * density * length という量で、無次元なので optical depth と同じです。
太陽組成のリュウグウ試料の発光スペクトルの計算
リュウグウの発光スペクトルについて計算してみましょう。計算方法は先ほどと同様で、Lodders 2003 のアバンダンスに従って、質量比に従って重みを変えて計算します。入射エネルギーは可変ですが、ここでは20keV程度の例を示します。
import matplotlib.colors as colors
import matplotlib.cm as cm
def Plot_CS_FluorLine_Ryugu(iene_kev = 20.0, pltflag = True, debug = False, outcsv = "fig_cs_flour_ryugu.csv"):
# e.g., Ryugu : density = 1.27 g/cm^3, thickness = 1um
density = 1.27 # 1.27 g/cm^3
length = 1e-2 # 100um
ev2kev = 1e-3
list_atomicnumber, list_El, list_mass, list_adundance = _get_solar()
# calculate CS_Total for several elements
total = np.sum(list_mass * list_adundance)
F = plt.figure(figsize=(12,8.))
ax = plt.subplot(1,1,1)
plt.xscale('linear')
plt.yscale('log')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency')
usercmap = plt.get_cmap('jet')
cNorm = colors.Normalize(vmin=0, vmax=len(list_atomicnumber))
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=usercmap)
f = open(outcsv,"w")
outstr ="e_keV,csf, lname\n"
f.write(outstr)
for i, (an, ma, abund) in enumerate(zip(list_atomicnumber, list_mass, list_adundance)):
elem = xraylib.AtomicNumberToSymbol(an)
c = scalarMap.to_rgba(i)
mass_fraction = (ma * abund) / total
partial_density = mass_fraction * density
list_lene_kev, list_csf, list_lname = Get_LineEnergy_CSF(an, iene_kev = iene_kev, density = partial_density, length = length)
for j, (lene, csf, lname) in enumerate(zip(list_lene_kev, list_csf, list_lname)):
outstr = str(lene) + "," + str(csf) + "," + str(lname) + "\n"
if debug: print(outstr)
f.write(outstr)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c, label = elem)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=10, borderaxespad=0.,fontsize=9)
else:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c)
f.close()
outfile = "fig_cs_flour_ryugu.png"
plt.savefig(outfile)
この関数を、次のように実行すると、
Plot_CS_FluorLine_Ryugu(iene_kev = 20.0, pltflag = True, outcsv = "fig_cs_flour_ryugu.csv")
次の図が生成されます。
これが、吸収が全くない場合の発光スペクトルになります。
次に、自己吸収がある場合について計算してみます。
from scipy.interpolate import interp1d
func_selfabs = interp1d(list_ene_ev * ev2kev, fx)
で、自己吸収の関数を作成している点だけが先ほどと異なります。
import matplotlib.colors as colors
import matplotlib.cm as cm
def Plot_CS_FluorLine_Ryugu_withSelfAbs(iene_kev = 20.0, pltflag = True, debug = False, outcsv = "fig_cs_flour_ryugu_selfabs.csv", ymin=1e-18, ymax=1e-2):
# e.g., Ryugu : density = 1.27 g/cm^3, thickness = 1um
density = 1.27 # 1.27 g/cm^3
length = 1 # 1cm, too thick just to visualize the absorption effect
ev2kev = 1e-3
list_atomicnumber, list_El, list_mass, list_adundance = _get_solar()
# calculate CS_Total for several elements
total = np.sum(list_mass * list_adundance)
# calculate self absorption function f(x)
xmin = 100 # eV
xmax = iene_kev * 1e3
xbins = 2e3
list_ene_ev = np.linspace(xmin,xmax,int(xbins))
input_ene_ev = list_ene_ev[-1]
reduce_flag = True
theta = np.pi/4
phi = np.pi/4
print("density = ", density, " [g/cm^3] length = ", length, " [cm]")
list_x, fx = _Get_x_fx(Get_CS_Total_Lo2003, input_ene_ev, list_ene_ev, density, length, reduce_flag, theta, phi)
from scipy.interpolate import interp1d
func_selfabs = interp1d(list_ene_ev * ev2kev, fx)
F = plt.figure(figsize=(14,8.))
ax = plt.subplot(1,2,1)
plt.xscale('linear')
plt.yscale('log')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency w/o SelfAbs')
plt.ylim(ymin, ymax)
ax = plt.subplot(1,2,2)
plt.xscale('linear')
plt.yscale('log')
plt.ylim(ymin, ymax)
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency w/ SelfAbs')
usercmap = plt.get_cmap('jet')
cNorm = colors.Normalize(vmin=0, vmax=len(list_atomicnumber))
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=usercmap)
f = open(outcsv,"w")
outstr ="e_keV,csf, absfactor, abscsf, lname\n"
f.write(outstr)
for i, (an, ma, abund) in enumerate(zip(list_atomicnumber, list_mass, list_adundance)):
elem = xraylib.AtomicNumberToSymbol(an)
c = scalarMap.to_rgba(i)
mass_fraction = (ma * abund) / total
partial_density = mass_fraction * density
list_lene_kev, list_csf, list_lname = Get_LineEnergy_CSF(an, iene_kev = iene_kev, density = partial_density, length = length)
for j, (lene, csf, lname) in enumerate(zip(list_lene_kev, list_csf, list_lname)):
selfabs_factor = func_selfabs(lene)
selfabs_csf = selfabs_factor * csf
outstr = str(lene) + "," + str(csf) + "," + str(selfabs_factor) + "," + str(selfabs_csf) + "," + str(lname) + "\n"
if debug: print(outstr)
f.write(outstr)
ax = plt.subplot(1,2,1)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c, label = elem)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=10, borderaxespad=0.,fontsize=9)
else:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c)
ax = plt.subplot(1,2,2)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=selfabs_csf, color=c, label = elem)
else:
plt.vlines(x=lene, ymin=0, ymax=selfabs_csf, color=c)
f.close()
outfile = "fig_cs_flour_ryugu_selfabs.png"
plt.savefig(outfile)
この関数を、
Plot_CS_FluorLine_Ryugu_withSelfAbs(iene_kev = 20.0, pltflag = True, debug = False, outcsv = "fig_cs_flour_ryugu_selfabs.csv")
で実行すると、
のように、自己吸収なし(左)と自己吸収あり(右)になります。ここでは、あえて吸収が大きくなるように厚みを大きめにしたパラメータの例で計算しています。
太陽組成の降着円盤の発光スペクトルの計算
降着円盤の蛍光X線のスペクトルは、リュウグウ試料と密度と厚みが違うだけで計算できるので、これも見てみよう。(厳密には、降着円盤は温度によって電離度が変わることや、回転運動や一般相対論的な効果も考慮する必要がある場合があります。)
import matplotlib.colors as colors
import matplotlib.cm as cm
def Plot_CS_FluorLine_AccretionDisk_withSelfAbs(iene_kev = 20.0, pltflag = True, debug = False, outcsv = "fig_cs_flour_adisk_selfabs.csv", ymin=1e-18, ymax=1e-2):
Avogadro_constant = 6.02e23
number_density = 1e17 # [1/cm^3] fig8 of https://arxiv.org/pdf/1908.07272.pdf
H_1mol = 1.0 # g/mol
density = H_1mol * (number_density / Avogadro_constant) # g/cm^3
ev2kev = 1e-3
length = 1e6 # [cm]
nh = number_density * length # [1/cm^2]
list_atomicnumber, list_El, list_mass, list_adundance = _get_solar()
# calculate CS_Total for several elements
total = np.sum(list_mass * list_adundance)
# calculate self absorption function f(x)
xmin = 100 # eV
xmax = iene_kev * 1e3
xbins = 2e3
list_ene_ev = np.linspace(xmin,xmax,int(xbins))
input_ene_ev = list_ene_ev[-1]
reduce_flag = True
theta = np.pi/4
phi = np.pi/4
print("density = ", density, " [g/cm^3] length = ", length, " [cm]", " Nh = ", nh, " [/cm^2]")
list_x, fx = _Get_x_fx(Get_CS_Total_Lo2003, input_ene_ev, list_ene_ev, density, length, reduce_flag, theta, phi)
from scipy.interpolate import interp1d
func_selfabs = interp1d(list_ene_ev * ev2kev, fx)
F = plt.figure(figsize=(14,8.))
ax = plt.subplot(1,2,1)
plt.xscale('linear')
plt.yscale('log')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency w/o SelfAbs')
plt.ylim(ymin, ymax)
ax = plt.subplot(1,2,2)
plt.xscale('linear')
plt.yscale('log')
plt.ylim(ymin, ymax)
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency w/ SelfAbs')
usercmap = plt.get_cmap('jet')
cNorm = colors.Normalize(vmin=0, vmax=len(list_atomicnumber))
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=usercmap)
f = open(outcsv,"w")
outstr ="e_keV,csf, absfactor, abscsf, lname\n"
f.write(outstr)
for i, (an, ma, abund) in enumerate(zip(list_atomicnumber, list_mass, list_adundance)):
elem = xraylib.AtomicNumberToSymbol(an)
c = scalarMap.to_rgba(i)
mass_fraction = (ma * abund) / total
partial_density = mass_fraction * density
list_lene_kev, list_csf, list_lname = Get_LineEnergy_CSF(an, iene_kev = iene_kev, density = partial_density, length = length)
for j, (lene, csf, lname) in enumerate(zip(list_lene_kev, list_csf, list_lname)):
selfabs_factor = func_selfabs(lene)
selfabs_csf = selfabs_factor * csf
outstr = str(lene) + "," + str(csf) + "," + str(selfabs_factor) + "," + str(selfabs_csf) + "," + str(lname) + "\n"
if debug: print(outstr)
f.write(outstr)
ax = plt.subplot(1,2,1)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c, label = elem)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=10, borderaxespad=0.,fontsize=9)
else:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c)
ax = plt.subplot(1,2,2)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=selfabs_csf, color=c, label = elem)
else:
plt.vlines(x=lene, ymin=0, ymax=selfabs_csf, color=c)
f.close()
outfile = "fig_cs_flour_adisk_selfabs.png"
plt.savefig(outfile)
これを、次のパラメータで実行すると、下記の図が生成されます。
Plot_CS_FluorLine_AccretionDisk_withSelfAbs(iene_kev = 20.0, pltflag = True, debug = False, outcsv = "fig_cs_flour_adisk_selfabs.csv")
X線入射窓とX線吸収体の効果
次に、より現実的な効果として、X線の入射窓やX線検出器の吸収体の効果についても計算方法を確認しておきます。
np.exp(-np.array([xraylib.CS_Total(an,E) for E in list_ene_ev * ev2kev])densitythickness) の部分で、xraylib.CS_Total(原子番号,エネルギー)により、そのエネルギーのその原子の吸収の全断面積を計算します。
X線の吸収率の計算の準備
# 透過率の計算
def Get_CS_Total_from_One(list_ene_ev, an, density, thickness):
# xraylib.CS_Total is Cross sections (cm2/g)
ev2kev = 1e-3
trans = np.exp(-np.array([xraylib.CS_Total(an,E) for E in list_ene_ev * ev2kev])*density*thickness)
return trans
# 単元素の透過率の計算
def Filter_One(list_ene_ev, el = "Al", thickness=330e-9 * 1e2):
list_ene_ev = np.array(list_ene_ev)
an = xraylib.SymbolToAtomicNumber(el)
afilter = {"an":an,"density":xraylib.ElementDensity(an),"thickness":thickness}
trans = Get_CS_Total_from_One(list_ene_ev, afilter["an"], afilter["density"], afilter["thickness"])
return trans
# 複数元素の透過率の計算
def Filter_Elements(list_ene_ev, list_el = ["Cr","Co","Cu"], list_massfrac=[1,1,1], thickness=0.1):
list_ene_ev = np.array(list_ene_ev)
list_massfrac = np.array(list_massfrac)
if not len(list_el) == len(list_massfrac):
print("ERROR length mismatch", list_el, list_massfrac)
return -1
list_fracdensity = []
for el, massfrac in zip(list_el, list_massfrac):
an = xraylib.SymbolToAtomicNumber(el)
fracdensity = xraylib.ElementDensity(an) * massfrac / np.sum(list_massfrac)
list_fracdensity.append(fracdensity)
list_trans = []
for el, fracdensity in zip(list_el, list_fracdensity):
an = xraylib.SymbolToAtomicNumber(el)
afilter = {"an":an,"density":fracdensity,"thickness":thickness}
onetrans = Get_CS_Total_from_One(list_ene_ev, afilter["an"], afilter["density"], afilter["thickness"])
list_trans.append(onetrans)
list_trans = np.array(list_trans)
trans = np.ones(len(list_ene_ev))
for one in list_trans:
trans = trans * one
return trans
の関数を定義して、動作確認します。
# 単元素の吸収の計算例 (A)
list_ene_ev = [100,1000,6000,1e4, 2e4] # eV
Filter_One(list_ene_ev, el = "Al", thickness=0.01)
# array([0.00000000e+00, 1.28200738e-14, 4.45281212e-02, 4.92896624e-01, 9.11291329e-01])
# 単元素の吸収の計算例 (Aと同じになることを確認)
list_ene_ev = [100,1000,6000,1e4, 2e4] # eV
Filter_Elements(list_ene_ev, list_el = ["Al"], list_massfrac=[1], thickness=0.01)
# array([0.00000000e+00, 1.28200738e-14, 4.45281212e-02, 4.92896624e-01, 9.11291329e-01])
# 複数元素の吸収の計算例 (Aと同じになることを確認)
list_ene_ev = [100,1000,6000,1e4, 2e4] # eV
Filter_Elements(list_ene_ev, list_el = ["Al","Al","Al"], list_massfrac=[1,1,1], thickness=0.01)
# array([0.00000000e+00, 1.28200738e-14, 4.45281212e-02, 4.92896624e-01,9.11291329e-01])
このように、Filter_One という関数で単元素の計算、Filter_Elements という関数で、複数元素の指定した質量比で、透過率を計算する。
X線吸収率のプロットの例 (空気 + Alフィルタ + BiのX線吸収体)
まずは、簡単なフィルタの透過率、Biの吸収体の吸収率を計算してみます。
xmin = 1e3 # eV
xmax = 20e3 # eV
xbins = 1e2
ev2kev = 1e-3
list_ene_ev = np.linspace(xmin,xmax,int(xbins))
F = plt.figure(figsize=(8,7))
ax = plt.subplot(1,1,1)
plt.xscale('linear')
plt.yscale('linear')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (eV)')
plt.ylabel(r'Efficiency')
# Air 5cm
thickness = 5 # cm
filter_Air = Filter_Elements(list_ene_ev, list_el = ["N","O","Ar"], list_massfrac=[75.524,23.139,1.288], thickness=thickness)
plt.errorbar(list_ene_ev, filter_Air, label = "Air 5cm")
# Al 300nm
thickness = 300e-9 * 1e2 # 300nm = 3e-5 cm
filter_Al = Filter_Elements(list_ene_ev, list_el = ["Al"], list_massfrac=[1], thickness=thickness)
plt.errorbar(list_ene_ev, filter_Al, label = "Al 110nm x 3 (100mK, 3K, 50K)")
# Bi 4um
thickness = 4e-6 * 1e2 # 4e-6m x 100 = 4e-4 cm
abs_Bi = 1.0 - Filter_Elements(list_ene_ev, list_el = ["Bi"], list_massfrac=[1], thickness=thickness)
plt.errorbar(list_ene_ev, abs_Bi, label = "Bi 4um (100mK)")
# all
eff_all = filter_Air * filter_Al * abs_Bi
plt.errorbar(list_ene_ev, eff_all, label = "Total")
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=10, borderaxespad=0.,fontsize=10)
outfile = "fig_filter.png"
plt.savefig(outfile)
X線吸収率のプロットの例(Heガス + Alフィルタ + LUXEL フィルタ + Bi吸収体)
次に、もう少し高度なフィルターを計算してみます。
# download data
fileURL="https://drive.google.com/file/d/1x1Z9R4IYR4tlKP7B746c-nQqMRKuKcfs/view?usp=sharing" # download LUXEL_filter_HT_large.csv
luxel_file="local_LUXEL_filter_HT_large.csv"
gdown.download(fileURL, luxel_file, quiet=False, fuzzy=True)
def _trans_luxel_ht_window(ene_ev, csvfile="./local_LUXEL_filter_HT_large.csv"):
from scipy.interpolate import interp1d
df = pd.read_csv(csvfile)
ene = np.array(df['ev'].values,dtype=float)
tra = np.array(df['trans'].values,dtype=float)
f = interp1d(ene, tra)
tras = f(ene_ev)/100.
return tras
xmin = 1e3 # eV
xmax = 20e3 # eV
xbins = 1e2
ev2kev = 1e-3
list_ene_ev = np.linspace(xmin,xmax,int(xbins))
F = plt.figure(figsize=(8,7))
ax = plt.subplot(1,1,1)
plt.xscale('linear')
plt.yscale('linear')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (eV)')
plt.ylabel(r'Efficiency')
# He 5cm
thickness = 5 # cm
filter_He = Filter_Elements(list_ene_ev, list_el = ["He"], list_massfrac=[1.0], thickness=thickness)
plt.errorbar(list_ene_ev, filter_He, label = "He 5cm")
# LUXEL Filter
filter_luxel = _trans_luxel_ht_window(list_ene_ev)
plt.errorbar(list_ene_ev, filter_luxel, label = "LUXEL HT Filter")
# Al 300nm
thickness = 300e-9 * 1e2 # 300nm = 3e-5 cm
filter_Al = Filter_Elements(list_ene_ev, list_el = ["Al"], list_massfrac=[1], thickness=thickness)
plt.errorbar(list_ene_ev, filter_Al, label = "Al 110nm x 3 (100mK, 3K, 50K)")
# Bi 4um
thickness = 4e-6 * 1e2 # 4e-6m x 100 = 4e-4 cm
abs_Bi = 1.0 - Filter_Elements(list_ene_ev, list_el = ["Bi"], list_massfrac=[1], thickness=thickness)
plt.errorbar(list_ene_ev, abs_Bi, label = "Bi 4um (100mK)")
# all
eff_all = filter_He * filter_luxel * filter_Al * abs_Bi
plt.errorbar(list_ene_ev, eff_all, label = "Total")
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=2, borderaxespad=0.,fontsize=10)
outfile = "fig_filter.png"
plt.savefig(outfile)
現実的なスペクトルの例
X線入射窓付きでリュウグウ試料を観測した場合のスペクトル
最後に、これまでのを駆使して、自己吸収とX線入射窓付きで、BiのX線吸収体で観測した場合の例を計算してみます。
import matplotlib.colors as colors
import matplotlib.cm as cm
def Plot_CS_FluorLine_Ryugu_withSelfAbsFilter(iene_kev = 20.0, pltflag = True, debug = False, outcsv = "fig_cs_flour_ryugu_selfabsfilter.csv", ymin=1e-18, ymax=1e-2):
# e.g., Ryugu : density = 1.27 g/cm^3, thickness = 1um
density = 1.27 # 1.27 g/cm^3
length = 1 # 1cm, too thick just to visualize the absorption effect
ev2kev = 1e-3
list_atomicnumber, list_El, list_mass, list_adundance = _get_solar()
# calculate CS_Total for several elements
total = np.sum(list_mass * list_adundance)
# calculate self absorption function f(x)
xmin = 100 # eV
xmax = iene_kev * 1e3
xbins = 2e3
list_ene_ev = np.linspace(xmin,xmax,int(xbins))
input_ene_ev = list_ene_ev[-1]
reduce_flag = True
theta = np.pi/4
phi = np.pi/4
print("density = ", density, " [g/cm^3] length = ", length, " [cm]")
list_x, fx = _Get_x_fx(Get_CS_Total_Lo2003, input_ene_ev, list_ene_ev, density, length, reduce_flag, theta, phi)
from scipy.interpolate import interp1d
func_selfabs = interp1d(list_ene_ev * ev2kev, fx)
# Air 5cm
thickness = 5 # cm
filter_Air = Filter_Elements(list_ene_ev, list_el = ["N","O","Ar"], list_massfrac=[75.524,23.139,1.288], thickness=thickness)
# # He 5cm
# thickness = 5 # cm
# filter_He = Filter_Elements(list_ene_ev, list_el = ["He"], list_massfrac=[1.0], thickness=thickness)
# LUXEL Filter
filter_luxel = _trans_luxel_ht_window(list_ene_ev)
# Al 300nm
thickness = 300e-9 * 1e2 # 300nm = 3e-5 cm
filter_Al = Filter_Elements(list_ene_ev, list_el = ["Al"], list_massfrac=[1], thickness=thickness)
# Bi 4um
thickness = 4e-6 * 1e2 # 4e-6m x 100 = 4e-4 cm
abs_Bi = 1.0 - Filter_Elements(list_ene_ev, list_el = ["Bi"], list_massfrac=[1], thickness=thickness)
# all
eff_all = filter_Air * filter_luxel * filter_Al * abs_Bi
func_eff = interp1d(list_ene_ev * ev2kev, eff_all)
F = plt.figure(figsize=(14,8.))
ax = plt.subplot(1,2,1)
plt.xscale('linear')
plt.yscale('log')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency w/o SelfAbs+Filter')
plt.ylim(ymin, ymax)
ax = plt.subplot(1,2,2)
plt.xscale('linear')
plt.yscale('log')
plt.ylim(ymin, ymax)
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency w/ SelfAbs+Filter')
usercmap = plt.get_cmap('jet')
cNorm = colors.Normalize(vmin=0, vmax=len(list_atomicnumber))
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=usercmap)
f = open(outcsv,"w")
outstr ="e_keV,csf, eff, effcsf, lname\n"
f.write(outstr)
for i, (an, ma, abund) in enumerate(zip(list_atomicnumber, list_mass, list_adundance)):
elem = xraylib.AtomicNumberToSymbol(an)
c = scalarMap.to_rgba(i)
mass_fraction = (ma * abund) / total
partial_density = mass_fraction * density
list_lene_kev, list_csf, list_lname = Get_LineEnergy_CSF(an, iene_kev = iene_kev, density = partial_density, length = length)
for j, (lene, csf, lname) in enumerate(zip(list_lene_kev, list_csf, list_lname)):
selfabs_factor = func_selfabs(lene)
eff_factor = func_eff(lene)
eff_csf = eff_factor * selfabs_factor * csf
outstr = str(lene) + "," + str(csf) + "," + str(eff_factor) + "," + str(eff_csf) + "," + str(lname) + "\n"
if debug: print(outstr)
f.write(outstr)
ax = plt.subplot(1,2,1)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c, label = elem)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=10, borderaxespad=0.,fontsize=9)
else:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c)
ax = plt.subplot(1,2,2)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=eff_csf, color=c, label = elem)
else:
plt.vlines(x=lene, ymin=0, ymax=eff_csf, color=c)
f.close()
outfile = "fig_cs_flour_ryugu_selfabsfilter.png"
plt.savefig(outfile)
図(左)が全く吸収がないスペクトルで、図(右)がフィルターによる吸収とBiによる光電吸収の効果を考えたスペクトルになってます。
X線入射窓付でNIST610を見た場合のスペクトル
import matplotlib.colors as colors
import matplotlib.cm as cm
def Plot_CS_FluorLine_NIST610_withSelfAbsFilter(iene_kev = 20.0, pltflag = True, debug = False, outcsv = "fig_cs_flour_nist610_selfabsfilter.csv", ymin=1e-12, ymax=1e1):
# e.g., NIST610 density = 2.65 # g/cm^3, thickness = 0.1 # 0.1 cm
density = 2.65 # g/cm^3,
length = 0.1 # 0.1 cm
list_atomicnumber, El, list_mass, list_mass_fraction, list_adundance = _get_nist610()
ev2kev = 1e-3
total = np.sum(list_mass * list_adundance)
# calculate self absorption function f(x)
xmin = 100 # eV
xmax = iene_kev * 1e3
xbins = 2e3
list_ene_ev = np.linspace(xmin,xmax,int(xbins))
input_ene_ev = list_ene_ev[-1]
reduce_flag = True
theta = np.pi/4
phi = np.pi/4
print("density = ", density, " [g/cm^3] length = ", length, " [cm]")
list_x, fx = _Get_x_fx(Get_CS_Total_NIST610, input_ene_ev, list_ene_ev, density, length, reduce_flag, theta, phi)
from scipy.interpolate import interp1d
func_selfabs = interp1d(list_ene_ev * ev2kev, fx)
# Air 5cm
thickness = 5 # cm
filter_Air = Filter_Elements(list_ene_ev, list_el = ["N","O","Ar"], list_massfrac=[75.524,23.139,1.288], thickness=thickness)
# # He 5cm
# thickness = 5 # cm
# filter_He = Filter_Elements(list_ene_ev, list_el = ["He"], list_massfrac=[1.0], thickness=thickness)
# LUXEL Filter
filter_luxel = _trans_luxel_ht_window(list_ene_ev)
# Al 300nm
thickness = 300e-9 * 1e2 # 300nm = 3e-5 cm
filter_Al = Filter_Elements(list_ene_ev, list_el = ["Al"], list_massfrac=[1], thickness=thickness)
# Bi 4um
thickness = 4e-6 * 1e2 # 4e-6m x 100 = 4e-4 cm
abs_Bi = 1.0 - Filter_Elements(list_ene_ev, list_el = ["Bi"], list_massfrac=[1], thickness=thickness)
# all
eff_all = filter_Air * filter_luxel * filter_Al * abs_Bi
func_eff = interp1d(list_ene_ev * ev2kev, eff_all)
F = plt.figure(figsize=(14,8.))
ax = plt.subplot(1,2,1)
plt.xscale('linear')
plt.yscale('log')
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency w/o SelfAbs+Filter')
plt.ylim(ymin, ymax)
ax = plt.subplot(1,2,2)
plt.xscale('linear')
plt.yscale('log')
plt.ylim(ymin, ymax)
plt.grid(linestyle='dotted',alpha=0.5)
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'Efficiency w/ SelfAbs+Filter')
usercmap = plt.get_cmap('jet')
cNorm = colors.Normalize(vmin=0, vmax=len(list_atomicnumber))
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=usercmap)
f = open(outcsv,"w")
outstr ="e_keV,csf, eff, effcsf, lname\n"
f.write(outstr)
for i, (an, ma, abund) in enumerate(zip(list_atomicnumber, list_mass, list_adundance)):
elem = xraylib.AtomicNumberToSymbol(an)
c = scalarMap.to_rgba(i)
mass_fraction = (ma * abund) / total
partial_density = mass_fraction * density
list_lene_kev, list_csf, list_lname = Get_LineEnergy_CSF(an, iene_kev = iene_kev, density = partial_density, length = length)
for j, (lene, csf, lname) in enumerate(zip(list_lene_kev, list_csf, list_lname)):
selfabs_factor = func_selfabs(lene)
eff_factor = func_eff(lene)
eff_csf = eff_factor * selfabs_factor * csf
outstr = str(lene) + "," + str(csf) + "," + str(eff_factor) + "," + str(eff_csf) + "," + str(lname) + "\n"
if debug: print(outstr)
f.write(outstr)
ax = plt.subplot(1,2,1)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c, label = elem)
plt.legend(bbox_to_anchor=(0., 1.01, 1., 0.01), loc='lower left',ncol=10, borderaxespad=0.,fontsize=9)
else:
plt.vlines(x=lene, ymin=0, ymax=csf, color=c)
ax = plt.subplot(1,2,2)
if j == 0:
plt.vlines(x=lene, ymin=0, ymax=eff_csf, color=c, label = elem)
else:
plt.vlines(x=lene, ymin=0, ymax=eff_csf, color=c)
f.close()
outfile = "fig_cs_flour_nist610_selfabsfilter.png"
plt.savefig(outfile)
Plot_CS_FluorLine_NIST610_withSelfAbsFilter(iene_kev = 20.0, pltflag = True, debug = False, outcsv = "fig_cs_flour_nist610_selfabsfilter.csv")
NIST610についての、図(左)が全く吸収がないスペクトルで、図(右)がフィルターによる吸収とBiによる光電吸収の効果を考えたスペクトルになってます。
まとめ
X線の自己吸収や、蛍光X線の発光率などは、一見すると複雑な計算に見えるが、分解して考えると一つ一つは簡単な計算です。
- 自己吸収が効くかどうかは、f(x)を計算して調べる
- 蛍光X線の発光については、 xraylib.CS_FluorLine_Kissel_Cascade(an, xraylib_lname, iene_kev) * density * length の部分で、 「蛍光X線生成断面積 (flourescence cross section)」 を計算する
- 化合物や混合物は、各元素の個数密度(質量密度)を計算して、断面積をその密度の重みで計算する
- フィルターやX線吸収体の効率については、xraylib.CS_Total を用いて全吸収断面積を計算する。
あたりを考えて計算すればよいです。
最後、断面積の計算は一般的に、
exp(-質量吸収/発光断面積[cm^2/g] * density [g/cm^3] * length [cm])
となるので、density [g/cm^3] * length [cm] はセットで未知数として入ってくるので、密度と奥行きを独立に決めるには、吸収だけでは難しく、別の独立した観測量(e.g., 輝線)が必要になります。