xraylib の使い方

ここでは蛍光X線分析のためのデータベース xraylib の使い方について紹介します。

主に python をベースの解説です。X線の基礎については知っていることを前提とはせずに、できる限り背景やツールの使い方から不定性も含めて解説してみます。(訂正コメントなど歓迎です。)

xraylib を用いて自己吸収について計算したい方は、

を参考にしてください。

X線と xraylib に関する基本事項

まずは xraylib がどんなものか雰囲気を掴むために、web インターフェースがあるので、

を使ってみよう。雰囲気としては、このようにX線に関する様々な定数や関数が取り揃えてられて、好きな言語で扱う場合のコードの例も表示してくれる優れものです。

エネルギーや相互作用断面積などの源泉情報はどこにあるのか？は大切なことで、それらは

に書かれています。

例えば、

xraylib/doc/xraydoc.txt

kissel_pe.dat
=============

Partial and Total Photoelectric cross sections.

 This data set contains data calculated by Lynn Kissel
 (lkissel@llnl.gov). The original data is available from
 ftp://www-phys.llnl.gov/rayleigh
 and an interactive WWW interface is in 
 http://www-phys.llnl.gov/V_Div/scattering/elastic.html

のように、Partial and Total Photoelectric cross sections については、この文献を使ってます、などと明記されていますので、必要に応じて源泉情報も確認しましょう。

X線と物質の相互作用や既存のデータベースなど

X線と物質は、「光電吸収」「コンプトン散乱」「(電子陽電子)対生成」の主に３つの反応があります。もっと高いエネルギーのガンマ線帯域まで行くと、「光核反応」が起こります。(注:「(電子陽電子)対生成」と「光核反応」は1MeV以上で起こる現象なので、主にX線分光向けのソフトウェアである xraylib にはガンマ線帯域の効果は入っていません。)

xraylib は、何かというと、主に地上での中性物質の蛍光X線ためのデータベースを統合したソフトウェア群という感じです。もっと厳密に知りたい人は、

などを参照ください。

xraylib の前に、昔からあるデータベースを簡単に紹介します。手計算レベルであれば今でもこれらで十分な場合も多いです。

X-ray Data Booklet
- 全元素の結合エネルギー、蛍光X線エネルギーなど
- 昔はX線の研究者なら手元に一冊、と言われた本

X-ray Transition Energies Database
- 蛍光X線の吸収端のエネルギーをGUIでささっと出してくれる。
- X-ray Data Booklet の GUI 版に近い。(同じ機関ではなく、NISTとバークレーが作ってるので微妙に違う?のかもしれません。)

NIST XCOM
- エネルギー帯域 : 主に硬X線からガンマ線帯域 (1 keV - 100 GeV)
- 光電吸収、コンプトン散乱(coherent, incoherent で弁別)、対生成(Nucler field, Electric filed で弁別)の反応断面積

webelements
- 電子配置、binding energy などを調べるのに便利な化学者

Filter Transmission
- エネルギー帯域 : 主に軟X線帯域 ( 0.01 keV - 30 GeV)
- X線フィルターの透過率の計算に便利

XrayDB
- xraylib とかなり似ていて、比較的最近開発されているもの。主にpython+webインターフェース対応。
- https://millenia.cars.aps.anl.gov/xraydb/ が web インターフェースで、周期表からパパッと全体像が見渡せて、かなり使いやすい。
- xraylib と XrayDB の違いは、xraylib は発光分光、XrayDBは結晶分光から進化したのではないかと。xraylib は、蛍光X線発光率が簡単に出せて、XrayDBはダーウィン角やアナライザ(アナライザ結晶のジャーゴン)の角度が充実してそうである(xraylibにも含まれている)。
- 周期表はすごく良い！ https://xraypy.github.io/XrayDB/periodictable.html#periodic-tables

おまけ
- 荷電粒子と物質の相互作用
  - PSTAR : 陽子(proton)と物質のdE/dx (ベーテブロッホの公式)
  - ASTAR : ヘリウム原子核(alpha particle)と dE/dx
  - ESTAR : 電子の物質の dE/dx
    - 電子は軽いので、陽子やアルファ線に比べて軽いので、制動放射も効くことに注意。
- ASTRO-H CookBook : 宇宙X線分光
  - 宇宙の精密X線分光向けにまとめられたもの(英語かつX線天文屋さん向けです)。

xraylib は、X線帯域の相互作用、特に蛍光X線に関して、C、Perl、Fortran 2003、IDL、Python、C++、Java、C#、Lua、Ruby、PHP、、Python-Numpy、Pascal の言語から使えるようにに整備されたものです。(Cが御本尊で他がバインディングかな)

現在(@2022.8.5)、xraylib を組み込んだソフトウェアは、

X-ray yield simulation (xys)
- 竜野氏製のX線の発光率の計算プログラム (TES向けに作られたもの)
- 断面積の説明

があります。SPring-8やJ-PARCでのTES応用の中で作成されたもので、GUIで、X線のフィルタ、標的、自己吸収なども計算してくれます。

X線の呼び方と flourescene yield について

ツールの使い方の前に、X線の呼び方だけは知っておく必要があります(ご存知の方は読み飛ばしてください)。

「シーグバーン記法」という 1924年にノーベル化学賞を受賞した Manne Siegbahn 氏の作成したものと、より整合性が取れる命名方として国際純正・応用化学連合 IUPAC 記法の２つが現在は共存している状態です。「蛍光X線スペクトルの読み方について」河合潤@京大など、X線分光の資料には初めに注意が必ず書いてある。
XrayDBの周期表はそのあたりも考慮してよく出来ていると思う。それを拝借して簡単にまとめると、

という感じです。IUPAC記法はエネルギーの低い順に数字をつけて、蛍光X線を表す場合は、遷移前後の軌道がわかるように明記する、というルールで、この表現の方が明確ではあります。

xraylibは、include/xraylib.h の中で定義していて、IUPAC の Table VIII.2 from Nomenclature system for X-ray spectroscopy に従った記法になっています。

に xraylib の中で使われているマクロの説明があります。

光電吸収で、K殻の電子が自由電子として原子の外に叩き出されると、K殻に一つ空の軌道ができるので、そこへ上の軌道の電子が落ちるときにX線が発生します。全確率を足すと、0.29(Ka2) + 0.54(Ka1) + 0.05(Kb3) + 0.09(Kb1) + 0.03(Kb2) = 1.0 となります。K線については大事な関係で、ライン同定の際によく使われるだいたいの目安として、

Ka2 : Ka1 = 1 : 2
Ka (=Ka1+Ka2) : Kb (=Kb1+Kb2+Kb3) = 10 : 1

も覚えておくとよいです。厳密には、分岐比は、量子力学の摂動論から <終状態の波動関数|相互作用ハミルトニアン|終状態の波動関数> で計算されるものです。より高次の効果を考えると他の遷移もありえますが、強度は小さいので普通の解析では寄与が小さいことが多いはずです。

もっと大事な効果として、
X-Ray Data Booklet Section 1.3 FLUORESCENCE YIELDS FOR K and L SHELLS
に記載がありますが、蛍光X線発生率 (flourescence yields) は大切です。定義は、空孔が出来たときに蛍光X線がでる確率です。

例えば、原子番号20のカルシウムのK殻の電子が外に叩き出されたすると、約20%が蛍光X線で、残りの80%はオージェ電子となります。オージェ電子のエネルギーも元素固有の値となるため、それを利用したオージェ電子分光（AES）という方法もあります。

ラインプロファイルや宇宙のX線プラズマとの違いなど

ラインプロファイルまで詳細に知りたい人は、少し注意が必要です。あくまで、データベースとして入っているのは、現実というよりは、理想的な中性物質の元素からの蛍光X線のエネルギーは強度です。実際は、化学的な効果(chemical shift)というのが、ラインの中心エネルギーや強度比に影響しますし、メタルの効果を踏まえた金属元素蛍光X線のラインプロファイルに書きましたが、MnのKalpha線などは、Voigt関数を7本や8本入れて現象論的に合わせ込むようなことをすることがあります。

宇宙のX線のプラズマ分光の場合は、冷たい降着円盤からの蛍光X線の場合はそのまま使えることも多いです。高階に電離したプラズマになると、AtomDBや、SPEXなど別のデータベースが必要になります。宇宙向けの詳細は ASTRO-H CookBook などご覧ください。

インストール方法

conda を用いる方法

に従えば基本的にはOKです。個人的には、conda が一番安定してる気がしますが、２年くらい前(2020年?)は、xraylibが既存の環境を破壊することがあったので、僕は今でも仮想環境を xraylib 専用に作ることにしていますが、最近は普通に conda install -c conda-forge xraylib で一発OKだと思います。

一応レベルで昔の例を挙げておきます。

install.sh

conda create -n xraylib # xraylibという名前のcondaの仮想環境を作成
conda activate xraylib # activate する。
conda install -c anaconda pyqt # qtを使わないなら入れなくても良い
conda install pandas # 必要なら
conda install -c conda-forge xraylib # xraylib の install 
conda install matplotlib # 必要なら
pip install scipy # should work # “conda install scipy” fails due to some mismatch

で入れていました。scipy も昔は mismatch がありましたが、今は大丈夫ではないでしょうか。

google colab で xraylib を使う方法

google colab で全サンプルを動かせるようにしてあります！!

コードを見れば分かる方は、xraylibをgoogle colabで動かした例をご覧ください。下記に要所の説明していきます。

xraylib を使ってみよう

例題を動かしてみる

ソースーコード :

に一通りの関数の動かし方の例があるので、まずはこれを動かしてみましょう。必要な関数や、使い方はこれをみれば大体はわかるかと思います。

サンプルは長いので、冒頭の一部だけ紹介します。

import xraylib
import math
import numpy as np

xraylib.XRayInit()
print("Example of python program using xraylib")
print("xraylib version: {}".format(xraylib.__version__))
print("Density of pure Al : {} g/cm3".format(xraylib.ElementDensity(13)))
print("Ca K-alpha Fluorescence Line Energy: {}".format(xraylib.LineEnergy(20, xraylib.KA_LINE)))
print("Fe partial photoionization cs of L3 at 6.0 keV: {}".format(xraylib.CS_Photo_Partial(26, xraylib.L3_SHELL,6.0)))
print("Zr L1 edge energy: {}".format(xraylib.EdgeEnergy(40, xraylib.L1_SHELL)))
print("Pb Lalpha XRF production cs at 20.0 keV (jump approx): {}".format(xraylib.CS_FluorLine(82, xraylib.LA_LINE,20.0)))
print("Pb Lalpha XRF production cs at 20.0 keV (Kissel): {}".format(xraylib.CS_FluorLine_Kissel(82, xraylib.LA_LINE,20.0)))
print("Bi M1N2 radiative rate: {}".format(xraylib.RadRate(83, xraylib.M1N2_LINE)))
print("U M3O3 Fluorescence Line Energy: {}".format(xraylib.LineEnergy(92, xraylib.M3O3_LINE)))
print("Ca(HCO3)2 Rayleigh cs at 10.0 keV: {}".format(xraylib.CS_Rayl_CP("Ca(HCO3)2",10.0)))

とすると、下記のようなメッセージが帰ってくれば正常です。

Example of python program using xraylib
xraylib version: 4.1.2
Density of pure Al : 2.6989 g/cm3
Ca K-alpha Fluorescence Line Energy: 3.6904905714961105
Fe partial photoionization cs of L3 at 6.0 keV: 22.444635592574553
Zr L1 edge energy: 2.5316
Pb Lalpha XRF production cs at 20.0 keV (jump approx): 11.414356266865491
Pb Lalpha XRF production cs at 20.0 keV (Kissel): 11.181532140402364
Bi M1N2 radiative rate: 0.53915
U M3O3 Fluorescence Line Energy: 4.11
Ca(HCO3)2 Rayleigh cs at 10.0 keV: 0.39754478321999853

分光についての略語を知らないと、意味がわからないことが多々あります。

XRF : X-Ray Fluorescence 蛍光X線
M1N2 : M殻の1番目の順位から、N殻の2番目の遷移 (IUPACの記法)
Kissel : Kissel Lynn さんというデータベース e.g., RTAB を作成されてる方のお名前であって、物理とは無関係。xraylib では、Kissel を明記した関数名が用意されていることがある。

原子番号とモル質量を csv ファイルに書き出して読み込む

最も簡単な例として、原子番号と、モル質量を csv ファイルに書き出して、pandas で読み込む例を示します。

# convert data into csv file, and open it with pandas
amin=1
amax=104 # this is a mix value of atomic number @2022.8.7
outfname = "xraylib_basic.csv" 
f = open(outfname, 'w')
outstr = "Z,El,atomicMass\n"
f.write(outstr)
for i in np.arange(amin,amax):
    atomicWeight = xraylib.AtomicWeight(int(i))
    atomicSymbol = xraylib.AtomicNumberToSymbol(int(i))
    outstr = str(i) + ","+str(atomicSymbol)+","+str(atomicWeight)+"\n"
#    print(outstr)
    f.write(outstr)
f.close()
data = pd.read_csv(outfname)

とすると、

のように、pandas.read_csv が、csv を解釈して、pandas の Dataframe オブジェクトにしてくれるので扱いがとても楽になります。この記事では csvファイルに保存して pandas で読み込むスタイルで書いておきました。

X線と物質の相互作用断面積

xraylib は、基本的に100keV以下の光電吸収が効く帯域で使うことを想定しているので、対生成の断面積は入っていない。下記では、「光電吸収」「コンプトン散乱」「レーリー散乱」の３つの散乱断面積と、それらの総和をプロットする例を示します。

# convert data into csv file, and open it with pandas
el_list=["He","O","Fe","Bi"] # element name list
for el in el_list:
    emin=0.1 # keV
    emax=100.0 # keV
    de=0.05; # keV
    outfname = "xraylib_cs_" + el + ".csv"; f = open(outfname, 'w')
    outstr = "ene,total,photo,rayl,compt\n"
    f.write(outstr)
    for ene in np.arange(emin,emax,de):
        an = xraylib.SymbolToAtomicNumber(el)
        cs_total = xraylib.CS_Total(an, ene)
        cs_photo = xraylib.CS_Photo(an, ene)
        cs_rayl = xraylib.CS_Rayl(an, ene)
        cs_compt = xraylib.CS_Compt(an, ene)
        outstr = str(ene) + ","+str(cs_total)+","+str(cs_photo) + ","+str(cs_rayl)+","+str(cs_compt)+"\n"
    #    print(outstr)
        f.write(outstr)
    f.close()
    csdata = pd.read_csv(outfname)

    fig = plt.figure(figsize =(6, 6))
    plt.title(el)
    plt.plot(csdata['ene'], csdata['total'],"k-",label='total')
    plt.plot(csdata['ene'], csdata['photo'],"r--",label='Photoionization')
    plt.plot(csdata['ene'], csdata['rayl'],"g--",label='Rayleigh scattering')
    plt.plot(csdata['ene'], csdata['compt'],"b--",label='Compton scattering')
    plt.grid(alpha=0.2, ls="dotted")
    plt.yscale("log")
    plt.xscale("log")
    plt.xlim(emin,emax)
    plt.ylim(1e-5,1e5)
    plt.ylabel(r" cross section ($cm^2/g$)")
    plt.xlabel("Energy (keV)")
    plt.legend(loc="best")
    plt.tight_layout()
    plt.show()

次のような図が作成されれば正常に動いている。

ヘリウムは約10keVで、光電吸収よりもコンプトン散乱の方が卓越する。原子番号が大きくなるにつれて、光電吸収の cross section が大きくなる。この傾向は、モーズリーの法則として 1913年に発見されていますが、精密分光の場合はこの法則の精度でよいことは少ないので、xraylib or データベースの値が必要になります。(光電吸収の反応断面積の厳密な計算は、高度な量子力学の計算が必要で、水素様原子ではない原子の場合には電子多体系とか難しい話が絡んできます。)

コンプトン散乱についても、xraylib.CS_Compt(an, ene) で計算されているように、エネルギーと原子を与えて計算されているのは、自由電子ではなくて束縛電子とのコンプトン散乱の反応断面積なので、元素の個性が入っています。(逆に束縛電子が初期状態で運動量がゼロではないことを利用して、物質中の電子の運度量をコンプトンプロファイルから調べるという物性研究もされています。)

蛍光X線のエネルギーと結合エネルギー

Ka1の蛍光X線のエネルギーと、K shell の binding energy をプロットしてみよう。この時に必要なのは、"IUPACで決められている名前"ですので、名前について、

から名前を確認しましょう。xraylib を使う時にはIUPAC記法をつかいましょう。ここでは、Ka1, Kedge, La1, Kedge のエネルギーをプロットする方法を紹介します。

amin=1; amax=104 # 104 is a mix value of atomic number @2022.8.7

# Ka1, Kedge
koutfname = "xraylib_Ka.csv" 
f = open(koutfname, 'w')
outstr = "Z,El,Kedge,Kalpha1\n"
f.write(outstr)
for i in range(amin,amax):
    atomicSymbol = xraylib.AtomicNumberToSymbol(i)    
    try:
        kedge = xraylib.EdgeEnergy(i, xraylib.K_SHELL) # see https://github.com/tschoonj/xraylib/wiki/Appendix-xraylib-macros#shell-macros
        ka1 = xraylib.LineEnergy(i, xraylib.KL3_LINE)  # Ka1, see https://github.com/tschoonj/xraylib/wiki/Appendix-xraylib-macros#siegbahn-x-ray-fluorescence-line-macros
    except:
        kedge = -1; ka1 = -1 # fill negative if NaN
#        print("warning ", i, atomicSymbol, sys.exc_info())                        
    outstr = str(i) + ","+str(atomicSymbol)+","+str(kedge)+","+str(ka1)+"\n"
    f.write(outstr)            
f.close()
kdata = pd.read_csv(koutfname)

# La1, L3edge
loutfname = "xraylib_La.csv" 
f = open(loutfname, 'w')
outstr = "Z,El,L3edge,Lalpha1\n"
f.write(outstr)
for i in range(amin,amax):
    atomicSymbol = xraylib.AtomicNumberToSymbol(i)    
    try:
        l3edge = xraylib.EdgeEnergy(i, xraylib.L3_SHELL) # see https://github.com/tschoonj/xraylib/wiki/Appendix-xraylib-macros#shell-macros
        la1 = xraylib.LineEnergy(i, xraylib.L3M5_LINE)  # Ka1, see https://github.com/tschoonj/xraylib/wiki/Appendix-xraylib-macros#siegbahn-x-ray-fluorescence-line-macros
    except:
        l3edge = -1; la1 = -1 # fill negative if NaN
#        print("warning ", i, atomicSymbol, sys.exc_info())                        
    outstr = str(i) + ","+str(atomicSymbol)+","+str(l3edge)+","+str(la1)+"\n"
    f.write(outstr)                        
f.close()
ldata = pd.read_csv(loutfname)

fig = plt.figure(figsize =(8, 6))
plt.plot(kdata['Z'], kdata['Kalpha1'],"r.",label=r'K$\alpha_1$')
plt.plot(kdata['Z'], kdata['Kedge'],"m.",label='K edge')
plt.plot(ldata['Z'], ldata['Lalpha1'],"g.",label=r'L$\alpha_1$')
plt.plot(ldata['Z'], ldata['L3edge'],"b.",label='L3 edge')
plt.grid(alpha=0.2, ls="dotted")
plt.yscale("log")
#plt.xscale("log")
#plt.xlim(emin,emax)
plt.ylim(0.1,150)
plt.ylabel(r"Energy (keV)")
plt.xlabel("atomic number")
plt.legend(loc="best")
plt.tight_layout()
plt.savefig(koutfname.replace(".csv",".png"))
plt.show()

実行して、

このような図が生成されれば正常です。xraylib の癖として、値が存在しない場合に、プログラムがコケるので、自分で try except を使って、エラーをキャッチする必要があります。この場合は、原子番号が小さい時に、K or L のエネルギーがX線帯域よりも低くなる場合に、xraylib では定義がされてないので、それを回避する必要があります。try except を使うことで、バグに気がつかないことも増えるので、try-except文の多用はお勧めはしませんが。

蛍光X線発生率 (flourescence yields)

次に蛍光X線発生率 (flourescence yields) をプロットしてみよう。

# Flourescence yield 
amin=1; amax=104 # 104 is a mix value of atomic number @2022.8.7

# Ka1, Kedge
fyoutfname = "xraylib_fy.csv" 
f = open(fyoutfname, 'w')
outstr = "Z,El,Kfy,L3fy\n"
f.write(outstr)

for i in range(amin,amax):
    atomicSymbol = xraylib.AtomicNumberToSymbol(i)    
    try: # Kshell
        kfy = xraylib.FluorYield(i, xraylib.K_SHELL)
    except:
        kfy = -1
        print("warning ", i, atomicSymbol, sys.exc_info())                        
    try: # L3 shell
        l3fy = xraylib.FluorYield(i, xraylib.L3_SHELL)
    except:
        l3fy = -1 # fill negative if NaN
        print("warning ", i, atomicSymbol, sys.exc_info())                        
                
    outstr = str(i) + ","+str(atomicSymbol)+","+str(kfy)+","+str(l3fy)+"\n"
    f.write(outstr)            
f.close()
fydata = pd.read_csv(fyoutfname)

fig = plt.figure(figsize =(7, 5))
plt.plot(fydata['Z'], fydata['Kfy'],"r.",label='K shell')
plt.plot(fydata['Z'], fydata['L3fy'],"b.",label='L3 shell')

plt.grid(alpha=0.2, ls="dotted")
plt.ylim(0,1)
plt.ylabel(r"Flourescence yield")
plt.xlabel("atomic number")
plt.legend(loc="best",fontsize=15)
plt.tight_layout()
plt.savefig(outfname.replace(".csv",".png"))
plt.show()

これを実行すると次の図が生成される。

原子番号が小さいと、オージェ電子でエネルギーを運び、原子番号が高くなると蛍光X線でエネルギーを解放するようになる。(注: これは中性の物質の話です。高階電離したプラズマや多価イオンでは、初期状態の電子配置が異なるので、オージェ電子の生成確率も変わります。)

蛍光X線の強度比

蛍光X線の強度比についても確認してみましょう。xraylib.RadRate という関数を使います。ここでの強度比とは、全部積分したら1になるように、各ラインごとに確率が与えられています。積分範囲はどこまでか？というと、空いた Shell を埋めるようなラインを全部積分したら 1 になるように定義しています。例えば、K shell の電子が光電吸収で自由空間に行ったとすると、このKshellの穴を埋めるように、KL2,KL3,KM2,KM3のいずれかが出るので、それらを全部足すと1になります。L3の場合は、L3M1、L3M2、L3M3、L3M4、L3M5、L3N1のいずれかのラインなので、それらを足すと１になります。これを確認してみましょう。

lines = ["KL2","KL3","KM2","KM3","L1L2","L1L3","L1M1","L1M2","L1M3","L1M4","L1M5","L1N1","L2L3","L2M1","L2M2","L2M3","L2M4","L2M5","L2N1","L3M1","L3M2","L3M3","L3M4","L3M5","L3N1"]
rates = []
L3_lines = [s for s in lines if s.startswith('L3')] # setect only lines from "L3"
for line in L3_lines:
    ls = "xraylib." + str(line) + "_LINE"    
    rates.append(xraylib.RadRate(26, eval(ls))) # assume Fe 
print(L3_lines) # ['L3M1', 'L3M2', 'L3M3', 'L3M4', 'L3M5', 'L3N1']
print(rates) # [0.095533, 8.5144e-05, 8.261e-05, 0.090617, 0.80735, 0.0063351]
print(np.sum(rates)) # 1.000002854

ここで、L3_lines が L3起因のラインリストで、それらの発生確率が rates で、その和がほぼ1になっていることがわかると思います。厳密には、もっと高次のサテライト線がありますので、そういう解析をやりたい場合は自分で追加する必要があります。

eval(ls)で、"xraylib.L3M1_LINE" などの文字列と、変数に変換しています。python では文字列を、eval で変換して変数にできます。詳しくは、https://note.nkmk.me/python-ast-literal-eval/ などご参照ください。

蛍光X線生成断面積 (flourescence cross section)

蛍光X線生成断面積については、 X-ray Fluorescence Cross Sections を一読することをお勧めします。

ここでは、

All CS_FluorLine* functions offer XRF cross sections for both K- and L-lines (provided the corresponding shells can be excited), but only the CS_FluorLine_Kissel* functions offer also the M-line XRF cross sections!
https://github.com/tschoonj/xraylib/wiki/The-xraylib-API-list-of-all-functions#x-ray-fluorescence-cross-sections

ということで、 M shell も計算できる CS_FluorLine_Kissel の例を示すことにします。

# Flourescence yield 
amin=1; amax=104 # 104 is a mix value of atomic number @2022.8.7
inenergies = [10.0, 20.0,40.0,60.0,80.0,100.0] # keV 

fig = plt.figure(figsize =(8, 6))

for inenergy in inenergies:
    # Ka1, Kedge
    cskoutfname = "xraylib_csk"+str(inenergy)+"keV.csv" 
    f = open(cskoutfname, 'w')
    outstr = "Z,El,csk\n"
    f.write(outstr)
    for i in range(amin,amax):
        atomicSymbol = xraylib.AtomicNumberToSymbol(i)    
        try: # Kshell
            csk = xraylib.CS_FluorLine_Kissel_Cascade(i, xraylib.KL3_LINE, inenergy) # Ka1            
        except:
            csk = -1
            print("warning ", i, atomicSymbol, sys.exc_info())                        

        outstr = str(i) + ","+str(atomicSymbol)+","+str(csk)+"\n"
        f.write(outstr)            

    f.close()

    cskdata = pd.read_csv(cskoutfname)
    plt.plot(cskdata['Z'], cskdata['csk'],".",label=r'K$\alpha_1$ at $E$ = ' + str(inenergy) + " keV")

    plt.grid(alpha=0.2, ls="dotted")
    #plt.xlim(0,inenergy)
    plt.yscale("log")
    plt.ylim(1e-5,1e2)
    plt.ylabel(r"X-ray flourescence production cross section ($cm^2/g$)")
    plt.xlabel("atomic number")
    plt.legend(loc="best",fontsize=15)
    plt.tight_layout()
plt.savefig(outfname.replace(".csv",".png"))
plt.show()

これは、入射するX線のエネルギーを 10keV, 20keV, 40keV, 60keV, 80keV, 100keV と振った時に、Kalpha1のX線の crosssection がどう変化するかを、原子番号ごとにプロットしたものです。

蛍光X線は吸収端(今の場合はKshellの吸収エネルギー)よりも入射エネルギーが大きくなるにつれて、光電吸収の相互作用断面積が E^-3.5 と下がってくるので、蛍光X線の収率も下がる。例えば、FeのKedgeは約7keVなので、10keVや20keVのX線はFeKa1の量に効くが、100keVのX線の寄与は２桁以上低いことになります。吸収端近傍のX線のフラックスが、蛍光X線の収量に効く、ということでもあります。
xraylib の CS の FluorLine、FluorShell、RadRate の関係は、

CS_{FluorLine} = CS_{FluorShell} \times RadRate

です。xraylibでは、

Care needs to be taken to use the correct shell and line macros!

と書かれておりますように、正しい、shell と line のマクロを使うように気をつけましょう。

厳選情報

データベースの厳選情報がズバッとわからないのですが、

あたりに書いてある。

そのさらに上流の情報は、

を使ってるような。

X線入射窓の透過率計算

X線の透過率を計算してみよう。xraylib.CS_Total を用いて、X線の相互作用断面積を cm^2/g の単位で取得したら、それに物質の密度と厚みを掛けて、指数関数の中にいれると透過率が計算できる。式で書くと、

        transmission = \exp (- CS\_TOTAL [\rm{cm}^2/g] \times thickness[\rm{cm}]  \times  density[\rm{g/cm}^3])

となる。吸収率は、1-透過率、で表せる。炭素とアルミニウムが 1um の時の透過率を計算する例になります。

# trans
# convert data into csv file, and open it with pandas
el_list=["C","Al"] # element name list
thickness_list = [0.0001,0.0001]     # cm 
for el,thickness in zip(el_list,thickness_list):
    emin=0.8 # keV
    emax=20.0 # keV
    de=0.05; # keV
    outfname = "xraylib_filt_" + el + ".csv"; f = open(outfname, 'w')
    outstr = "ene,trans,abs\n"
    f.write(outstr)
    an = xraylib.SymbolToAtomicNumber(el)
    density = xraylib.ElementDensity(an) # g/cm^3
    print("thickness = ", thickness * 10 * 1e3, " [um]")
    print(el, "density = ", density, " g/cm^3")
    for ene in np.arange(emin,emax,de):
        cs_total = xraylib.CS_Total(an, ene)
        trans = np.exp(-1.0 * cs_total * thickness * density)        
        abs =  1.0 - trans
        outstr = str(ene) + ","+str(trans)+","+str(abs)+"\n"
        f.write(outstr)
    f.close()
    trdata = pd.read_csv(outfname)

    fig = plt.figure(figsize =(7, 4))
    plt.title(el)
    plt.plot(trdata['ene'], trdata['trans'],"k-",label=el)
    plt.grid(alpha=0.2, ls="dotted")
#    plt.yscale("log")
    plt.xscale("log")
    plt.xlim(emin,emax)
    plt.ylim(0,1.1)
    plt.ylabel(r"Transmission")
    plt.xlabel("Energy (keV)")
    plt.legend(loc="best")
    plt.tight_layout()
    plt.savefig(outfname.replace(".csv",".png"))
    plt.show()

次のような図がでれば正しく動いています。

X線吸収体の吸収率計算

X線の吸収率 = 1 - 透過率、なので、これも簡単に計算できる。

# absorption 
# convert data into csv file, and open it with pandas
el_list=["Au","Bi"] # element name list
thickness_list = [0.0004,0.0004]     # cm 
for el,thickness in zip(el_list,thickness_list):
    emin=0.8 # keV
    emax=20.0 # keV
    de=0.05; # keV
    outfname = "xraylib_filt_" + el + ".csv"; f = open(outfname, 'w')
    outstr = "ene,trans,abs\n"
    f.write(outstr)
    an = xraylib.SymbolToAtomicNumber(el)
    density = xraylib.ElementDensity(an) # g/cm^3
    print("thickness = ", thickness * 10 * 1e3, " [um]")
    print(el, "density = ", density, " g/cm^3")
    for ene in np.arange(emin,emax,de):
        cs_total = xraylib.CS_Total(an, ene)
        trans = np.exp(-1.0 * cs_total * thickness * density)        
        abs =  1.0 - trans
        outstr = str(ene) + ","+str(trans)+","+str(abs)+"\n"
        f.write(outstr)
    f.close()
    trdata = pd.read_csv(outfname)

    fig = plt.figure(figsize =(7, 4))
    plt.title(el)
    plt.plot(trdata['ene'], trdata['abs'],"k-",label=el)
    plt.grid(alpha=0.2, ls="dotted")
#    plt.yscale("log")
    plt.xscale("log")
    plt.xlim(emin,emax)
    plt.ylim(0,1.1)
    plt.ylabel(r"Absorption")
    plt.xlabel("Energy (keV)")
    plt.legend(loc="best")
    plt.tight_layout()
    plt.savefig(outfname.replace(".csv",".png"))
    plt.show()

このように、金とビスマスの吸収率も簡単に計算できる。

クライン-仁科の公式

クライン-仁科の公式（英: Klein-Nishina's formula）は、量子電磁力学の最低次での、束縛を受けていない自由電子による光散乱の散乱断面積を与える関係式です。自由電子なので、特に束縛電子の影響など考えないので、xraylibで計算する必要は特にない(原子番号などの原子固有の依存性が入ってない)のですが、簡単なプロット方法だけ紹介しておきます。

import math
ene_list = [30, 100, 511, 662,1500] # keV
#scattering_angle_list = [0, 45, 90, 135, 180] # deg
azimuthal_angle_list = [0, 45, 90, 135, 180] # deg

for azimuthal_angle in azimuthal_angle_list:
  fig = plt.figure(figsize =(7, 4))

  plt.title("Klein-Nishina differential scattering cross section\n" + "Azimuthal angle (deg) = " +str(azimuthal_angle))

  for ene in ene_list:
    tag=str(ene) + "keV_" + str(azimuthal_angle) + "_deg"
    outfname = "xraylib_KN_" + tag + ".csv"; f = open(outfname, 'w')
    outstr = "scattering_angle,cs\n"
    f.write(outstr)

    for scattering_angle in np.arange(0,181,1):
        cs = xraylib.DCSP_KN(ene, math.radians(scattering_angle), math.radians(azimuthal_angle))
        outstr = str(scattering_angle) + "," + str(cs) + "\n"
        f.write(outstr)
    f.close()    
    data = pd.read_csv(outfname)

    plt.plot(data['scattering_angle'], data['cs'],"-",label=tag)
    plt.grid(alpha=0.2, ls="dotted")
    plt.yscale("log")
    plt.xscale("linear")
    plt.ylim(1e-4,1e-1)

    plt.ylabel(r"Cross Section [cm$^2$/g/sr]")
    plt.xlabel("Scattering Angle (deg)")
    plt.legend(loc="best")
    plt.tight_layout()

  plt.savefig(outfname.replace(".csv",".png"))
  plt.show()

この例では、

cs = xraylib.DCSP_KN(ene, math.radians(scattering_angle), math.radians(azimuthal_angle))

という部分で、実体は、

に、

/*////////////////////////////////////////////////////////////////////
//                                                                  //
//      Klein-Nishina differential scattering cross section         // 
//                for polarized beam (barn)                         //
//                                                                  //
//          E : Energy (keV)                                        //
//          theta : scattering polar angle (rad)                    //
//          phi : scattering azimuthal angle (rad)                  //
//                                                                  //
/////////////////////////////////////////////////////////////////// */
double DCSP_KN(double E, double theta, double phi, xrl_error **error)
{ 
  double k0_k, k_k0, k_k0_2, cos_th, sin_th, cos_phi;
  
  if (E <= 0.0) {
    xrl_set_error_literal(error, XRL_ERROR_INVALID_ARGUMENT, NEGATIVE_ENERGY);
    return 0.0;
  }

  cos_th = cos(theta);
  sin_th = sin(theta);
  cos_phi = cos(phi);
  
  k0_k = 1.0 + (1.0 - cos_th) * E / MEC2 ;
  k_k0 = 1.0 / k0_k;
  k_k0_2 = k_k0 * k_k0;
  
  return (RE2/2.) * k_k0_2 * (k_k0 + k0_k - 2 * sin_th * sin_th 
			      * cos_phi * cos_phi);
}

で計算されているものである。
online calculator によると、単位は、barn ではなくて、cm2/g/sr に変換されているようで、どこかで単位を換算していると思われる。