第4章: X線光電子分光法 (XPS: X-ray Photoelectron Spectroscopy)

4.2 化学シフトとピーク同定

4.2.1 化学シフトの起源

原子の化学状態（酸化状態、配位環境）が変化すると、内殻電子の結合エネルギーがシフトします。これを化学シフト（chemical shift）と呼びます。

4.2.2 代表的な元素のXPSピーク

コード例1: XPSスペクトルのシミュレーションとピーク同定

import numpy as np
import matplotlib.pyplot as plt

def gaussian_peak(x, amplitude, center, width):
    """
    ガウス型ピークを生成

    Parameters:
    -----------
    x : array
        結合エネルギー (eV)
    amplitude : float
        ピーク高さ
    center : float
        ピーク中心 (eV)
    width : float
        半値全幅（FWHM, eV）

    Returns:
    --------
    peak : array
        ガウス型ピークの強度
    """
    sigma = width / (2 * np.sqrt(2 * np.log(2)))
    peak = amplitude * np.exp(-((x - center)**2) / (2 * sigma**2))
    return peak

def simulate_xps_c1s_spectrum():
    """
    C 1s XPSスペクトルをシミュレート（複数の化学状態）
    """
    # 結合エネルギー範囲
    BE = np.linspace(280, 295, 1500)

    # C 1s ピーク（4つの化学状態）
    C_CC = gaussian_peak(BE, amplitude=1000, center=284.5, width=1.2)   # C-C, C-H
    C_CO = gaussian_peak(BE, amplitude=300, center=286.5, width=1.3)    # C-O
    C_O = gaussian_peak(BE, amplitude=150, center=288.0, width=1.4)     # C=O
    COO = gaussian_peak(BE, amplitude=80, center=289.5, width=1.5)      # O-C=O

    # 総スペクトル
    total_spectrum = C_CC + C_CO + C_O + COO

    # ノイズ
    noise = np.random.normal(0, 10, len(BE))
    observed_spectrum = total_spectrum + noise

    # プロット
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

    # スペクトル全体
    ax1.plot(BE, observed_spectrum, 'k-', linewidth=1.5, label='観測スペクトル')
    ax1.fill_between(BE, C_CC, alpha=0.3, color='blue', label='C-C, C-H (284.5 eV)')
    ax1.fill_between(BE, C_CO, alpha=0.3, color='green', label='C-O (286.5 eV)')
    ax1.fill_between(BE, C_O, alpha=0.3, color='orange', label='C=O (288.0 eV)')
    ax1.fill_between(BE, COO, alpha=0.3, color='red', label='O-C=O (289.5 eV)')
    ax1.set_xlabel('結合エネルギー (eV)', fontsize=12)
    ax1.set_ylabel('強度 (cps)', fontsize=12)
    ax1.set_title('C 1s XPSスペクトル（多成分）', fontsize=14, fontweight='bold')
    ax1.legend(loc='upper right', fontsize=9)
    ax1.grid(alpha=0.3)
    ax1.invert_xaxis()  # XPSの慣例（高エネルギー → 低エネルギー）

    # 各成分の分離
    ax2.plot(BE, C_CC, 'b-', linewidth=2, label='C-C, C-H')
    ax2.plot(BE, C_CO, 'g-', linewidth=2, label='C-O')
    ax2.plot(BE, C_O, 'orange', linewidth=2, label='C=O')
    ax2.plot(BE, COO, 'r-', linewidth=2, label='O-C=O')
    ax2.axvline(284.5, color='blue', linestyle='--', alpha=0.5)
    ax2.axvline(286.5, color='green', linestyle='--', alpha=0.5)
    ax2.axvline(288.0, color='orange', linestyle='--', alpha=0.5)
    ax2.axvline(289.5, color='red', linestyle='--', alpha=0.5)
    ax2.set_xlabel('結合エネルギー (eV)', fontsize=12)
    ax2.set_ylabel('強度 (cps)', fontsize=12)
    ax2.set_title('各化学状態の分離', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(alpha=0.3)
    ax2.invert_xaxis()

    plt.tight_layout()
    plt.show()

    print("C 1s ピーク同定:")
    print("  284.5 eV: C-C, C-H（炭化水素骨格）")
    print("  286.5 eV: C-O（エーテル、アルコール）")
    print("  288.0 eV: C=O（カルボニル基）")
    print("  289.5 eV: O-C=O（カルボキシル基）")

# 実行
simulate_xps_c1s_spectrum()

4.3 ピークフィッティングとデコンボリューション

4.3.1 ピーク形状の選択

XPSピークは、ガウス型とローレンツ型の混合であるVoigt関数（またはその近似であるGaussian-Lorentzian関数）で表されます。

4.3.2 Shirleyバックグラウンド

XPSスペクトルには、非弾性散乱された電子によるバックグラウンドが重畳します。David A. Shirley（1972）が提案したShirleyバックグラウンドが広く使用されます。

コード例2: Shirleyバックグラウンド減算とピークフィッティング

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def shirley_background(x, y, tol=1e-5, max_iter=50):
    """
    Shirleyバックグラウンドを計算

    Parameters:
    -----------
    x : array
        結合エネルギー (eV)
    y : array
        観測強度
    tol : float
        収束判定基準
    max_iter : int
        最大反復回数

    Returns:
    --------
    background : array
        Shirleyバックグラウンド
    """
    # データの昇順ソート（BE: 高 → 低）
    if x[0] < x[-1]:
        x = x[::-1]
        y = y[::-1]

    # 初期化
    background = np.zeros_like(y)
    y_max = np.max(y)
    y_min = np.min(y)

    # 反復計算
    for iteration in range(max_iter):
        # 各点でのバックグラウンド計算
        for i in range(1, len(x)):
            integral = np.trapz(y[:i] - background[:i], x[:i])
            background[i] = y_min + (y_max - y_min) * integral / np.trapz(y - background, x)

        # 収束判定
        if iteration > 0:
            change = np.max(np.abs(background - background_old))
            if change < tol:
                break
        background_old = background.copy()

    return background

def voigt_approximation(x, amplitude, center, sigma, gamma):
    """
    Voigt関数の近似（Gaussian-Lorentzian混合）

    Parameters:
    -----------
    x : array
        結合エネルギー
    amplitude : float
        ピーク高さ
    center : float
        ピーク中心
    sigma : float
        ガウス成分の幅
    gamma : float
        ローレンツ成分の幅

    Returns:
    --------
    voigt : array
        Voigt関数の値
    """
    gaussian = np.exp(-((x - center)**2) / (2 * sigma**2))
    lorentzian = gamma**2 / ((x - center)**2 + gamma**2)
    voigt = amplitude * (0.7 * gaussian + 0.3 * lorentzian)
    return voigt

def multi_peak_fit(x, y, initial_params):
    """
    複数ピークのフィッティング

    Parameters:
    -----------
    x : array
        結合エネルギー
    y : array
        強度
    initial_params : list of tuples
        各ピークの初期パラメータ [(A1, c1, s1, g1), (A2, c2, s2, g2), ...]

    Returns:
    --------
    fitted_params : array
        フィッティングされたパラメータ
    fitted_peaks : list
        各ピークの曲線
    """
    def multi_voigt(x, *params):
        """複数のVoigtピークの合計"""
        n_peaks = len(params) // 4
        result = np.zeros_like(x)
        for i in range(n_peaks):
            A, c, s, g = params[i*4:(i+1)*4]
            result += voigt_approximation(x, A, c, s, g)
        return result

    # 初期パラメータをフラット化
    p0 = [p for peak in initial_params for p in peak]

    # フィッティング
    popt, pcov = curve_fit(multi_voigt, x, y, p0=p0, maxfev=10000)

    # 各ピークを再構成
    n_peaks = len(popt) // 4
    fitted_peaks = []
    for i in range(n_peaks):
        A, c, s, g = popt[i*4:(i+1)*4]
        peak = voigt_approximation(x, A, c, s, g)
        fitted_peaks.append((peak, A, c, s, g))

    return popt, fitted_peaks

# シミュレーションデータ（Si 2p: Si + SiO2）
BE = np.linspace(96, 108, 1200)

# 真のピーク
Si_metal = voigt_approximation(BE, amplitude=800, center=99.3, sigma=0.4, gamma=0.2)
SiO2 = voigt_approximation(BE, amplitude=500, center=103.5, sigma=0.5, gamma=0.25)
true_spectrum = Si_metal + SiO2

# バックグラウンド（指数減衰型）
background_true = 100 * np.exp(-(BE - 96) / 10)

# 観測スペクトル
noise = np.random.normal(0, 15, len(BE))
observed = true_spectrum + background_true + noise

# Shirleyバックグラウンド減算
shirley_bg = shirley_background(BE, observed)
corrected = observed - shirley_bg

# ピークフィッティング
initial_params = [
    (800, 99.3, 0.4, 0.2),   # Si metal
    (500, 103.5, 0.5, 0.25)  # SiO2
]
fitted_params, fitted_peaks = multi_peak_fit(BE, corrected, initial_params)

# プロット
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 元のスペクトル
axes[0, 0].plot(BE, observed, 'k-', label='観測スペクトル', linewidth=1.5)
axes[0, 0].plot(BE, shirley_bg, 'r--', label='Shirleyバックグラウンド', linewidth=2)
axes[0, 0].set_xlabel('結合エネルギー (eV)')
axes[0, 0].set_ylabel('強度 (cps)')
axes[0, 0].set_title('Si 2p スペクトル（生データ）')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)
axes[0, 0].invert_xaxis()

# バックグラウンド減算後
axes[0, 1].plot(BE, corrected, 'o', markersize=3, alpha=0.5, label='BG減算後')
fitted_total = sum([peak[0] for peak in fitted_peaks])
axes[0, 1].plot(BE, fitted_total, 'r-', linewidth=2, label='フィッティング合計')
axes[0, 1].set_xlabel('結合エネルギー (eV)')
axes[0, 1].set_ylabel('強度 (cps)')
axes[0, 1].set_title('Shirley BG減算後とフィッティング')
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.3)
axes[0, 1].invert_xaxis()

# 各成分の分離
axes[1, 0].plot(BE, corrected, 'k-', alpha=0.3, label='観測データ')
colors = ['blue', 'green']
labels = ['Si metal (99.3 eV)', 'SiO₂ (103.5 eV)']
for i, (peak, A, c, s, g) in enumerate(fitted_peaks):
    axes[1, 0].fill_between(BE, peak, alpha=0.5, color=colors[i], label=labels[i])
    axes[1, 0].axvline(c, color=colors[i], linestyle='--', linewidth=1.5)
axes[1, 0].set_xlabel('結合エネルギー (eV)')
axes[1, 0].set_ylabel('強度 (cps)')
axes[1, 0].set_title('ピーク分離（デコンボリューション）')
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.3)
axes[1, 0].invert_xaxis()

# 残差
residual = corrected - fitted_total
axes[1, 1].plot(BE, residual, 'purple', linewidth=1)
axes[1, 1].axhline(0, color='black', linestyle='--', linewidth=1)
axes[1, 1].fill_between(BE, residual, alpha=0.3, color='purple')
axes[1, 1].set_xlabel('結合エネルギー (eV)')
axes[1, 1].set_ylabel('残差 (cps)')
axes[1, 1].set_title('フィッティング残差')
axes[1, 1].grid(alpha=0.3)
axes[1, 1].invert_xaxis()

plt.tight_layout()
plt.show()

# フィッティング結果の出力
print("ピークフィッティング結果:")
for i, (peak, A, c, s, g) in enumerate(fitted_peaks):
    area = np.trapz(peak, BE)
    print(f"  ピーク {i+1}: 中心 = {c:.2f} eV, 面積 = {area:.1f}")

4.4 定量分析と感度係数

4.4.1 XPS定量分析の原理

XPSスペクトルのピーク面積から、表面組成の原子濃度を決定できます。

4.4.2 Scofield相対感度係数

J.H. Scofield（1976）により理論計算された光イオン化断面積に基づく感度係数が広く使用されます。

コード例3: XPS定量分析と組成決定

import numpy as np
import matplotlib.pyplot as plt

def xps_quantification(peak_areas, sensitivity_factors):
    """
    XPS定量分析により原子濃度を計算

    Parameters:
    -----------
    peak_areas : dict
        各元素のピーク面積 {'C': area_C, 'O': area_O, ...}
    sensitivity_factors : dict
        各元素の相対感度係数 {'C': RSF_C, 'O': RSF_O, ...}

    Returns:
    --------
    atomic_concentrations : dict
        各元素の原子濃度（at%）
    """
    # 規格化強度の計算
    normalized_intensities = {}
    for element, area in peak_areas.items():
        RSF = sensitivity_factors[element]
        normalized_intensities[element] = area / RSF

    # 合計
    total = sum(normalized_intensities.values())

    # 原子濃度（at%）
    atomic_concentrations = {}
    for element, norm_int in normalized_intensities.items():
        atomic_concentrations[element] = (norm_int / total) * 100

    return atomic_concentrations

def plot_composition(atomic_concentrations):
    """
    組成を円グラフと棒グラフで可視化
    """
    elements = list(atomic_concentrations.keys())
    concentrations = list(atomic_concentrations.values())

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

    # 円グラフ
    colors = plt.cm.Set3(np.linspace(0, 1, len(elements)))
    wedges, texts, autotexts = ax1.pie(concentrations, labels=elements, autopct='%1.1f%%',
                                         colors=colors, startangle=90, textprops={'fontsize': 12})
    for autotext in autotexts:
        autotext.set_color('white')
        autotext.set_fontweight('bold')
    ax1.set_title('表面組成（原子濃度）', fontsize=14, fontweight='bold')

    # 棒グラフ
    ax2.bar(elements, concentrations, color=colors, edgecolor='black', linewidth=1.5)
    ax2.set_xlabel('元素', fontsize=12)
    ax2.set_ylabel('原子濃度 (at%)', fontsize=12)
    ax2.set_title('表面組成（棒グラフ）', fontsize=14, fontweight='bold')
    ax2.grid(alpha=0.3, axis='y')

    # 各元素の濃度を棒の上に表示
    for i, (elem, conc) in enumerate(zip(elements, concentrations)):
        ax2.text(i, conc + 1, f'{conc:.1f}%', ha='center', va='bottom', fontsize=11, fontweight='bold')

    plt.tight_layout()
    plt.show()

# 実行例：SiO2薄膜の分析
print("=== XPS定量分析例：SiO₂薄膜 ===\n")

# 測定されたピーク面積（任意単位）
peak_areas = {
    'Si': 12000,  # Si 2p
    'O': 28000,   # O 1s
    'C': 3000     # C 1s（炭素汚染）
}

# Scofield相対感度係数（Al Kα線）
sensitivity_factors = {
    'Si': 0.27,
    'O': 0.66,
    'C': 0.25
}

# 定量分析
atomic_conc = xps_quantification(peak_areas, sensitivity_factors)

print("ピーク面積:")
for elem, area in peak_areas.items():
    print(f"  {elem}: {area}")

print("\n相対感度係数:")
for elem, RSF in sensitivity_factors.items():
    print(f"  {elem}: {RSF}")

print("\n原子濃度:")
for elem, conc in atomic_conc.items():
    print(f"  {elem}: {conc:.2f} at%")

# 理論組成との比較（SiO2 = Si:O = 1:2 = 33.3:66.7）
Si_theory = 33.3
O_theory = 66.7
print(f"\n理論組成（SiO₂）: Si = 33.3 at%, O = 66.7 at%")
print(f"測定組成: Si = {atomic_conc['Si']:.2f} at%, O = {atomic_conc['O']:.2f} at%")
print(f"炭素汚染: C = {atomic_conc['C']:.2f} at%")

# プロット
plot_composition(atomic_conc)

# 汚染補正後の組成
Si_corrected = atomic_conc['Si'] / (atomic_conc['Si'] + atomic_conc['O']) * 100
O_corrected = atomic_conc['O'] / (atomic_conc['Si'] + atomic_conc['O']) * 100
print(f"\n炭素補正後: Si = {Si_corrected:.2f} at%, O = {O_corrected:.2f} at%")

4.5 深さ方向分析（Depth Profiling）

4.5.1 イオンスパッタリング法

Ar⁺イオンビームで試料表面をスパッタリングしながらXPS測定を繰り返すことで、深さ方向の組成プロファイルを取得できます。

4.5.2 非破壊角度分解XPS

検出角度を変化させることで、非破壊的に深さ情報を取得できます。

コード例4: 深さ方向プロファイルのシミュレーションと解析

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf

def depth_profile_simulation(depth, interface_position, interface_width, concentration_top, concentration_bottom):
    """
    界面を持つ深さ方向濃度プロファイルをシミュレート

    Parameters:
    -----------
    depth : array
        深さ (nm)
    interface_position : float
        界面位置 (nm)
    interface_width : float
        界面幅（拡散幅, nm）
    concentration_top : float
        表面層の濃度 (at%)
    concentration_bottom : float
        基板層の濃度 (at%)

    Returns:
    --------
    concentration : array
        各深さにおける濃度 (at%)
    """
    # erf関数で界面を表現
    concentration = concentration_bottom + (concentration_top - concentration_bottom) * \
                    0.5 * (1 - erf((depth - interface_position) / interface_width))
    return concentration

def simulate_sputter_depth_profiling():
    """
    スパッタリング深さ方向分析をシミュレート（SiO2/Si構造）
    """
    # 深さ範囲
    depth = np.linspace(0, 50, 200)  # 0-50 nm

    # SiO2層（0-20 nm）とSi基板（20 nm以降）
    Si_profile = depth_profile_simulation(depth, interface_position=20, interface_width=2,
                                          concentration_top=33, concentration_bottom=100)
    O_profile = depth_profile_simulation(depth, interface_position=20, interface_width=2,
                                         concentration_top=67, concentration_bottom=0)

    # スパッタリング測定点（離散的）
    sputter_times = np.array([0, 5, 10, 15, 20, 25, 30, 40, 50])  # スパッタリング時間（分）
    sputter_rate = 1.0  # nm/min
    measured_depths = sputter_times * sputter_rate

    # 測定濃度（ノイズ付き）
    Si_measured = np.interp(measured_depths, depth, Si_profile) + np.random.normal(0, 2, len(measured_depths))
    O_measured = np.interp(measured_depths, depth, O_profile) + np.random.normal(0, 2, len(measured_depths))

    # プロット
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

    # 理論プロファイル
    ax1.plot(depth, Si_profile, 'b-', linewidth=2, label='Si（理論）')
    ax1.plot(depth, O_profile, 'r-', linewidth=2, label='O（理論）')
    ax1.scatter(measured_depths, Si_measured, s=100, color='blue', marker='o',
                edgecolor='black', linewidth=1.5, label='Si（測定）', zorder=5)
    ax1.scatter(measured_depths, O_measured, s=100, color='red', marker='s',
                edgecolor='black', linewidth=1.5, label='O（測定）', zorder=5)
    ax1.axvline(20, color='green', linestyle='--', linewidth=2, label='界面位置 (20 nm)')
    ax1.set_xlabel('深さ (nm)', fontsize=12)
    ax1.set_ylabel('原子濃度 (at%)', fontsize=12)
    ax1.set_title('深さ方向プロファイル（SiO₂/Si）', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(alpha=0.3)
    ax1.set_xlim(0, 50)
    ax1.set_ylim(-5, 105)

    # 層構造の可視化
    ax2.fill_between([0, 20], 0, 100, alpha=0.3, color='red', label='SiO₂層 (20 nm)')
    ax2.fill_between([20, 50], 0, 100, alpha=0.3, color='blue', label='Si基板')
    ax2.text(10, 50, 'SiO₂', fontsize=16, fontweight='bold', ha='center')
    ax2.text(35, 50, 'Si', fontsize=16, fontweight='bold', ha='center')
    ax2.axvline(20, color='green', linestyle='--', linewidth=3, label='界面')
    ax2.set_xlabel('深さ (nm)', fontsize=12)
    ax2.set_ylabel('層構造', fontsize=12)
    ax2.set_title('試料構造（断面図）', fontsize=14, fontweight='bold')
    ax2.set_xlim(0, 50)
    ax2.set_yticks([])
    ax2.legend()
    ax2.grid(alpha=0.3, axis='x')

    plt.tight_layout()
    plt.show()

    # 界面厚さの推定
    interface_region = (measured_depths >= 15) & (measured_depths <= 25)
    if np.sum(interface_region) > 0:
        interface_width_est = np.max(measured_depths[interface_region]) - np.min(measured_depths[interface_region])
        print(f"推定界面幅: {interface_width_est:.1f} nm")

# 実行
simulate_sputter_depth_profiling()

4.6 Auger電子とピーク識別

4.6.1 Auger過程の原理

XPS測定では、光電子ピークに加えて、Auger電子ピークも観測されます。Auger電子は、内殻に空孔が生じた後の緩和過程で放出されます。

import numpy as np
import matplotlib.pyplot as plt

def simulate_xps_with_auger(x_ray_energy):
    """
    XPSスペクトルにおける光電子ピークとAuger電子ピークをシミュレート

    Parameters:
    -----------
    x_ray_energy : float
        X線エネルギー（eV）: Al Kα = 1486.6, Mg Kα = 1253.6

    Returns:
    --------
    スペクトルをプロット
    """
    # 結合エネルギー範囲（0-1500 eV）
    BE = np.linspace(0, 1500, 3000)

    # 光電子ピーク（結合エネルギーで表示）
    C_1s = gaussian_peak(BE, amplitude=800, center=284.5, width=1.2)
    O_1s = gaussian_peak(BE, amplitude=600, center=532.0, width=1.5)
    Si_2p = gaussian_peak(BE, amplitude=400, center=99.3, width=1.0)

    # Auger電子ピーク（運動エネルギーから結合エネルギーへ変換）
    # C KLL Auger: 運動エネルギー ≈ 270 eV（X線エネルギーに依存しない）
    C_KLL_kinetic = 270  # eV
    C_KLL_BE = x_ray_energy - C_KLL_kinetic
    C_KLL = gaussian_peak(BE, amplitude=150, center=C_KLL_BE, width=8)

    # O KLL Auger: 運動エネルギー ≈ 510 eV
    O_KLL_kinetic = 510
    O_KLL_BE = x_ray_energy - O_KLL_kinetic
    O_KLL = gaussian_peak(BE, amplitude=120, center=O_KLL_BE, width=10)

    # 総スペクトル
    total_spectrum = C_1s + O_1s + Si_2p + C_KLL + O_KLL

    # ノイズ
    noise = np.random.normal(0, 10, len(BE))
    observed = total_spectrum + noise

    # プロット
    fig, ax = plt.subplots(figsize=(14, 6))

    ax.plot(BE, observed, 'k-', linewidth=1.5, label='観測スペクトル', alpha=0.7)
    ax.fill_between(BE, C_1s, alpha=0.5, color='blue', label='C 1s (284.5 eV) 光電子')
    ax.fill_between(BE, O_1s, alpha=0.5, color='red', label='O 1s (532.0 eV) 光電子')
    ax.fill_between(BE, Si_2p, alpha=0.5, color='green', label='Si 2p (99.3 eV) 光電子')
    ax.fill_between(BE, C_KLL, alpha=0.5, color='purple', label=f'C KLL ({C_KLL_BE:.1f} eV) Auger')
    ax.fill_between(BE, O_KLL, alpha=0.5, color='orange', label=f'O KLL ({O_KLL_BE:.1f} eV) Auger')

    # 光電子とAuger電子の区別を強調
    ax.axvline(284.5, color='blue', linestyle='--', linewidth=1, alpha=0.5)
    ax.axvline(532.0, color='red', linestyle='--', linewidth=1, alpha=0.5)
    ax.axvline(99.3, color='green', linestyle='--', linewidth=1, alpha=0.5)
    ax.axvline(C_KLL_BE, color='purple', linestyle=':', linewidth=2)
    ax.axvline(O_KLL_BE, color='orange', linestyle=':', linewidth=2)

    ax.set_xlabel('結合エネルギー (eV)', fontsize=12)
    ax.set_ylabel('強度 (cps)', fontsize=12)
    ax.set_title(f'XPSワイドスキャン（X線エネルギー: {x_ray_energy:.1f} eV）', fontsize=14, fontweight='bold')
    ax.legend(loc='upper right', fontsize=10)
    ax.grid(alpha=0.3)
    ax.invert_xaxis()
    ax.set_xlim(1500, 0)

    plt.tight_layout()
    plt.show()

    print(f"X線エネルギー: {x_ray_energy:.1f} eV")
    print("\n光電子ピーク（結合エネルギー）:")
    print(f"  C 1s: 284.5 eV")
    print(f"  O 1s: 532.0 eV")
    print(f"  Si 2p: 99.3 eV")
    print("\nAuger電子ピーク（結合エネルギー表示）:")
    print(f"  C KLL: {C_KLL_BE:.1f} eV（運動エネルギー: {C_KLL_kinetic} eV）")
    print(f"  O KLL: {O_KLL_BE:.1f} eV（運動エネルギー: {O_KLL_kinetic} eV）")
    print("\n識別方法: X線エネルギーを変えると、Auger電子の結合エネルギー表示位置は変わるが、")
    print("         光電子の結合エネルギーは変わらない。")

# Al Kα線で測定
print("=== Al Kα線（1486.6 eV）での測定 ===")
simulate_xps_with_auger(x_ray_energy=1486.6)

# Mg Kα線で測定
print("\n=== Mg Kα線（1253.6 eV）での測定 ===")
simulate_xps_with_auger(x_ray_energy=1253.6)

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter

def xps_noise_reduction(BE, noisy_spectrum, window_length=11, polyorder=3):
    """
    Savitzky-Golayフィルタによるノイズ除去

    Parameters:
    -----------
    BE : array
        結合エネルギー
    noisy_spectrum : array
        ノイズを含むスペクトル
    window_length : int
        フィルタ窓の長さ（奇数）
    polyorder : int
        多項式の次数

    Returns:
    --------
    smoothed_spectrum : array
        平滑化されたスペクトル
    """
    smoothed = savgol_filter(noisy_spectrum, window_length=window_length, polyorder=polyorder)
    return smoothed

# シミュレーション：低カウントレートの測定（ノイズ大）
BE = np.linspace(280, 295, 1500)

# 真のスペクトル
true_spectrum = gaussian_peak(BE, amplitude=500, center=284.5, width=1.2) + \
                gaussian_peak(BE, amplitude=200, center=286.5, width=1.3)

# 高ノイズ
heavy_noise = np.random.normal(0, 30, len(BE))
noisy_spectrum = true_spectrum + heavy_noise

# 異なるパラメータでSGフィルタ適用
smoothed_sg5 = xps_noise_reduction(BE, noisy_spectrum, window_length=5, polyorder=2)
smoothed_sg11 = xps_noise_reduction(BE, noisy_spectrum, window_length=11, polyorder=3)
smoothed_sg21 = xps_noise_reduction(BE, noisy_spectrum, window_length=21, polyorder=3)

# プロット
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 元のノイズスペクトル
axes[0, 0].plot(BE, noisy_spectrum, 'gray', alpha=0.5, linewidth=0.5, label='ノイズあり')
axes[0, 0].plot(BE, true_spectrum, 'r-', linewidth=2, label='真のスペクトル')
axes[0, 0].set_xlabel('結合エネルギー (eV)')
axes[0, 0].set_ylabel('強度 (cps)')
axes[0, 0].set_title('元のスペクトル（高ノイズ）')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)
axes[0, 0].invert_xaxis()

# SGフィルタ（window_length=5）
axes[0, 1].plot(BE, noisy_spectrum, 'gray', alpha=0.3, linewidth=0.5)
axes[0, 1].plot(BE, smoothed_sg5, 'b-', linewidth=2, label='SG (window=5, order=2)')
axes[0, 1].plot(BE, true_spectrum, 'r--', linewidth=1.5, label='真のスペクトル')
axes[0, 1].set_xlabel('結合エネルギー (eV)')
axes[0, 1].set_ylabel('強度 (cps)')
axes[0, 1].set_title('Savitzky-Golay フィルタ（window=5）')
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.3)
axes[0, 1].invert_xaxis()

# SGフィルタ（window_length=11）
axes[1, 0].plot(BE, noisy_spectrum, 'gray', alpha=0.3, linewidth=0.5)
axes[1, 0].plot(BE, smoothed_sg11, 'g-', linewidth=2, label='SG (window=11, order=3)')
axes[1, 0].plot(BE, true_spectrum, 'r--', linewidth=1.5, label='真のスペクトル')
axes[1, 0].set_xlabel('結合エネルギー (eV)')
axes[1, 0].set_ylabel('強度 (cps)')
axes[1, 0].set_title('Savitzky-Golay フィルタ（window=11）')
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.3)
axes[1, 0].invert_xaxis()

# SGフィルタ（window_length=21）
axes[1, 1].plot(BE, noisy_spectrum, 'gray', alpha=0.3, linewidth=0.5)
axes[1, 1].plot(BE, smoothed_sg21, 'purple', linewidth=2, label='SG (window=21, order=3)')
axes[1, 1].plot(BE, true_spectrum, 'r--', linewidth=1.5, label='真のスペクトル')
axes[1, 1].set_xlabel('結合エネルギー (eV)')
axes[1, 1].set_ylabel('強度 (cps)')
axes[1, 1].set_title('Savitzky-Golay フィルタ（window=21）')
axes[1, 1].legend()
axes[1, 1].grid(alpha=0.3)
axes[1, 1].invert_xaxis()

plt.tight_layout()
plt.show()

# 誤差評価
mse_sg5 = np.mean((smoothed_sg5 - true_spectrum)**2)
mse_sg11 = np.mean((smoothed_sg11 - true_spectrum)**2)
mse_sg21 = np.mean((smoothed_sg21 - true_spectrum)**2)

print("平滑化の評価（平均二乗誤差）:")
print(f"  SG (window=5, order=2): MSE = {mse_sg5:.2f}")
print(f"  SG (window=11, order=3): MSE = {mse_sg11:.2f}（最良）")
print(f"  SG (window=21, order=3): MSE = {mse_sg21:.2f}")
print("\n推奨: window_length = 11-15, polyorder = 2-3")

import numpy as np
import matplotlib.pyplot as plt

def charge_shift_correction(BE, spectrum, reference_peak_position, true_reference_BE):
    """
    帯電シフトの補正

    Parameters:
    -----------
    BE : array
        測定された結合エネルギー
    spectrum : array
        測定スペクトル
    reference_peak_position : float
        参照ピークの観測位置 (eV)
    true_reference_BE : float
        参照ピークの真の結合エネルギー (eV)

    Returns:
    --------
    corrected_BE : array
        補正後の結合エネルギー
    shift : float
        帯電シフト量 (eV)
    """
    shift = reference_peak_position - true_reference_BE
    corrected_BE = BE - shift
    return corrected_BE, shift

# シミュレーション：絶縁体試料の帯電
BE_charged = np.linspace(280, 540, 2600)

# 帯電により +3.0 eV シフトしたスペクトル
charge_shift = 3.0
C_1s_charged = gaussian_peak(BE_charged, amplitude=800, center=284.5 + charge_shift, width=1.2)
O_1s_charged = gaussian_peak(BE_charged, amplitude=600, center=532.0 + charge_shift, width=1.5)
spectrum_charged = C_1s_charged + O_1s_charged + np.random.normal(0, 10, len(BE_charged))

# C 1s ピーク位置の検出（最大値）
C_1s_region = (BE_charged >= 284.5 + charge_shift - 5) & (BE_charged <= 284.5 + charge_shift + 5)
C_1s_observed_pos = BE_charged[C_1s_region][np.argmax(spectrum_charged[C_1s_region])]

# 帯電補正
corrected_BE, detected_shift = charge_shift_correction(BE_charged, spectrum_charged,
                                                        reference_peak_position=C_1s_observed_pos,
                                                        true_reference_BE=284.5)

# 補正後のスペクトル（同じ強度、横軸のみ補正）
spectrum_corrected = spectrum_charged

# プロット
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 補正前
ax1.plot(BE_charged, spectrum_charged, 'r-', linewidth=1.5, label='帯電スペクトル')
ax1.axvline(284.5 + charge_shift, color='blue', linestyle='--', linewidth=2, label=f'C 1s（観測: {284.5 + charge_shift:.1f} eV）')
ax1.axvline(532.0 + charge_shift, color='green', linestyle='--', linewidth=2, label=f'O 1s（観測: {532.0 + charge_shift:.1f} eV）')
ax1.set_xlabel('結合エネルギー（測定値, eV）', fontsize=12)
ax1.set_ylabel('強度 (cps)', fontsize=12)
ax1.set_title('帯電シフト前（補正前）', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)
ax1.invert_xaxis()

# 補正後
ax2.plot(corrected_BE, spectrum_corrected, 'b-', linewidth=1.5, label='補正後スペクトル')
ax2.axvline(284.5, color='blue', linestyle='--', linewidth=2, label='C 1s（補正後: 284.5 eV）')
ax2.axvline(532.0, color='green', linestyle='--', linewidth=2, label='O 1s（補正後: 532.0 eV）')
ax2.set_xlabel('結合エネルギー（補正値, eV）', fontsize=12)
ax2.set_ylabel('強度 (cps)', fontsize=12)
ax2.set_title('帯電シフト補正後', fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)
ax2.invert_xaxis()

plt.tight_layout()
plt.show()

print("帯電補正結果:")
print(f"  検出された帯電シフト: {detected_shift:.2f} eV")
print(f"  真の帯電シフト: {charge_shift:.2f} eV")
print(f"  C 1s ピーク位置: {C_1s_observed_pos:.2f} eV → 284.5 eV に補正")
print(f"  O 1s ピーク位置: {532.0 + charge_shift:.2f} eV → 532.0 eV に補正")

h_nu = 1486.6  # eV (Al Kα)
E_kinetic = 1202.1  # eV
phi = 4.5  # eV
E_binding = h_nu - E_kinetic - phi
print(f"C 1s 結合エネルギー: {E_binding:.1f} eV")

I_Si = 15000
I_O = 35000
RSF_Si = 0.27
RSF_O = 0.66

norm_Si = I_Si / RSF_Si
norm_O = I_O / RSF_O
total = norm_Si + norm_O

C_Si = (norm_Si / total) * 100
C_O = (norm_O / total) * 100

print(f"Si: {C_Si:.1f} at%")
print(f"O: {C_O:.1f} at%")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Fe 2p3/2 スペクトルのシミュレーション
BE_Fe = np.linspace(700, 720, 2000)

# 3つの化学状態
Fe0 = gaussian_peak(BE_Fe, amplitude=300, center=707.0, width=2.0)
Fe2 = gaussian_peak(BE_Fe, amplitude=500, center=709.5, width=2.0)
Fe3 = gaussian_peak(BE_Fe, amplitude=400, center=710.8, width=2.0)
observed_Fe = Fe0 + Fe2 + Fe3 + np.random.normal(0, 15, len(BE_Fe))

# フィッティング（ガウス関数を使用）
def three_gaussian(x, A1, c1, w1, A2, c2, w2, A3, c3, w3):
    return (gaussian_peak(x, A1, c1, w1) +
            gaussian_peak(x, A2, c2, w2) +
            gaussian_peak(x, A3, c3, w3))

p0 = [300, 707.0, 2.0, 500, 709.5, 2.0, 400, 710.8, 2.0]
popt, _ = curve_fit(three_gaussian, BE_Fe, observed_Fe, p0=p0, maxfev=10000)

# 各成分の再構成
Fe0_fit = gaussian_peak(BE_Fe, popt[0], popt[1], popt[2])
Fe2_fit = gaussian_peak(BE_Fe, popt[3], popt[4], popt[5])
Fe3_fit = gaussian_peak(BE_Fe, popt[6], popt[7], popt[8])

# ピーク面積
area_Fe0 = np.trapz(Fe0_fit, BE_Fe)
area_Fe2 = np.trapz(Fe2_fit, BE_Fe)
area_Fe3 = np.trapz(Fe3_fit, BE_Fe)
total_area = area_Fe0 + area_Fe2 + area_Fe3

# 各酸化状態の割合
ratio_Fe0 = (area_Fe0 / total_area) * 100
ratio_Fe2 = (area_Fe2 / total_area) * 100
ratio_Fe3 = (area_Fe3 / total_area) * 100

print("Fe酸化状態の割合:")
print(f"  Fe⁰（金属鉄）: {ratio_Fe0:.1f}%")
print(f"  Fe²⁺（FeO）: {ratio_Fe2:.1f}%")
print(f"  Fe³⁺（Fe₂O₃）: {ratio_Fe3:.1f}%")

# プロット
plt.figure(figsize=(12, 6))
plt.plot(BE_Fe, observed_Fe, 'ko', markersize=2, alpha=0.5, label='観測データ')
plt.fill_between(BE_Fe, Fe0_fit, alpha=0.5, color='blue', label=f'Fe⁰ ({ratio_Fe0:.1f}%)')
plt.fill_between(BE_Fe, Fe2_fit, alpha=0.5, color='green', label=f'Fe²⁺ ({ratio_Fe2:.1f}%)')
plt.fill_between(BE_Fe, Fe3_fit, alpha=0.5, color='red', label=f'Fe³⁺ ({ratio_Fe3:.1f}%)')
plt.xlabel('結合エネルギー (eV)', fontsize=12)
plt.ylabel('強度 (cps)', fontsize=12)
plt.title('Fe 2p₃/₂ 多成分フィッティング', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.gca().invert_xaxis()
plt.tight_layout()
plt.show()

C_1s_observed = 287.5  # eV
C_1s_reference = 284.5  # eV
Si_2p_observed = 102.3  # eV

charge_shift = C_1s_observed - C_1s_reference
Si_2p_corrected = Si_2p_observed - charge_shift

print(f"帯電シフト: {charge_shift:.1f} eV")
print(f"補正後 Si 2p: {Si_2p_corrected:.1f} eV（金属Si）")

import numpy as np

I_90 = 1000  # cps
I_30 = 200   # cps
lambda_imfp = 3.0  # nm
theta_90 = np.radians(90)
theta_30 = np.radians(30)

# 層厚さの計算
d = -lambda_imfp * np.log(I_30 / I_90) / (1/np.sin(theta_30) - 1/np.sin(theta_90))

print(f"SiO₂層厚さ: {d:.2f} nm")

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix
import seaborn as sns

# トレーニングデータ生成（特徴量: ピーク中心位置、FWHM、ピーク高さ）
np.random.seed(42)

# 金属（BE ≈ 99, FWHM ≈ 1.0, 高さ ≈ 800）
metal_features = np.random.normal([99.0, 1.0, 800], [0.3, 0.1, 50], (10, 3))

# 酸化物（BE ≈ 103, FWHM ≈ 1.5, 高さ ≈ 600）
oxide_features = np.random.normal([103.0, 1.5, 600], [0.3, 0.15, 40], (10, 3))

# 窒化物（BE ≈ 101, FWHM ≈ 1.2, 高さ ≈ 700）
nitride_features = np.random.normal([101.0, 1.2, 700], [0.3, 0.12, 45], (10, 3))

# データ統合
X = np.vstack([metal_features, oxide_features, nitride_features])
y = np.array([0]*10 + [1]*10 + [2]*10)  # 0: 金属, 1: 酸化物, 2: 窒化物

# トレーニング・テストデータ分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42, stratify=y)

# ランダムフォレスト分類器
rf_classifier = RandomForestClassifier(n_estimators=100, max_depth=10, random_state=42)
rf_classifier.fit(X_train, y_train)

# 予測
y_pred = rf_classifier.predict(X_test)

# 評価
class_names = ['金属', '酸化物', '窒化物']
print("分類レポート:")
print(classification_report(y_test, y_pred, target_names=class_names))

# 混同行列
cm = confusion_matrix(y_test, y_pred)
plt.figure(figsize=(8, 6))
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', xticklabels=class_names, yticklabels=class_names)
plt.xlabel('予測ラベル', fontsize=12)
plt.ylabel('真のラベル', fontsize=12)
plt.title('混同行列（XPS分類）', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# 特徴量重要度
importances = rf_classifier.feature_importances_
feature_names = ['ピーク中心 (eV)', 'FWHM (eV)', 'ピーク高さ (cps)']
plt.figure(figsize=(8, 5))
plt.barh(feature_names, importances, color='skyblue', edgecolor='black')
plt.xlabel('重要度', fontsize=12)
plt.title('特徴量の重要度', fontsize=14, fontweight='bold')
plt.grid(alpha=0.3, axis='x')
plt.tight_layout()
plt.show()

print("\n特徴量重要度:")
for name, importance in zip(feature_names, importances):
    print(f"  {name}: {importance:.3f}")

import numpy as np
import matplotlib.pyplot as plt

def identify_photoelectron_auger(BE_AlKa, intensity_AlKa, BE_MgKa, intensity_MgKa, threshold=1.0):
    """
    光電子ピークとAuger電子ピークを識別

    Parameters:
    -----------
    BE_AlKa, BE_MgKa : array
        各X線源での結合エネルギー
    intensity_AlKa, intensity_MgKa : array
        各スペクトルの強度
    threshold : float
        ピーク位置の変化閾値（eV）

    Returns:
    --------
    peak_types : dict
        各ピークのタイプ（'photoelectron' or 'Auger'）
    """
    # ピーク検出（簡易的に最大値検出）
    from scipy.signal import find_peaks

    peaks_AlKa, _ = find_peaks(intensity_AlKa, height=100, distance=50)
    peaks_MgKa, _ = find_peaks(intensity_MgKa, height=100, distance=50)

    # 各ピークの対応を探索
    peak_types = {}
    for i, peak_Al in enumerate(peaks_AlKa):
        BE_Al = BE_AlKa[peak_Al]

        # MgKaスペクトルで対応ピークを探す（±5 eV以内）
        matched = False
        for peak_Mg in peaks_MgKa:
            BE_Mg = BE_MgKa[peak_Mg]
            if abs(BE_Al - BE_Mg) < threshold:
                # 結合エネルギーが一致 → 光電子ピーク
                peak_types[f"Peak_{i+1} ({BE_Al:.1f} eV)"] = "Photoelectron"
                matched = True
                break

        if not matched:
            # 結合エネルギーが変化 → Auger電子ピーク
            # MgKaでの期待位置を計算
            E_shift = 1486.6 - 1253.6  # Al Kα - Mg Kα = 233 eV
            expected_BE_Mg = BE_Al - E_shift
            peak_types[f"Peak_{i+1} ({BE_Al:.1f} eV, Mg: {expected_BE_Mg:.1f} eV)"] = "Auger"

    return peak_types

# 実装例（コード例5のデータを使用）
# （省略: 実際にはコード例5のシミュレーションデータを使用）

print("識別アルゴリズム:")
print("  光電子: 結合エネルギーがX線エネルギーに依存しない")
print("  Auger: 結合エネルギー表示がX線エネルギーに依存する")