Skip to content

Physics and Material Constants

This section documents the physical foundations of the simulation, including target properties and spline-interpolated interaction cross-sections.

lib.physics_base

Fundamental Physical Constants and QED Formulations Module.

Provides standard CODATA / Particle Data Group (PDG) constants and exact quantum electrodynamics (QED) cross-section equations. All energy calculations are calibrated in keV.

adrv_bound(u, k)

Analytical antiderivative for dσ/du. Used to calculate bounds without slow numerical integration.

Parameters:

Name Type Description Default
u float

Ratio (E-E_prime)/E of energy transfer to initial photon energy.

 required
k float

Ratio of initial photon energy to the electron rest mass.

 required

Returns:

Type Description
float

Antiderivative of dσ/du in barns.
Source code in lib/physics_base.py 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def adrv_bound(u, k):
    """
    Analytical antiderivative for dσ/du.
    Used to calculate bounds without slow numerical integration.

    Parameters
    ----------
    u : float
        Ratio (E-E_prime)/E of energy transfer to initial photon energy.
    k : float
        Ratio of initial photon energy to the electron rest mass.

    Returns
    -------
    float
        Antiderivative of dσ/du in barns.

    """
    if u >= 1.0 or u <= 0.0:
        return 0.0
    return se / (k**3) * (1 / (1 - u) - ((u * k) ** 2) / 2 + u * ((1 + k) ** 2) + np.log(1 - u) * (3 - (1 - k) ** 2))

sample_bound_compton_energy(E, E_bind)

Samples the scattered photon energy E_prime using Kahn's method, recommitting the sampled kinematics into the energy-transfer fraction scale u = (E - E_prime) / E bounded by strict atomic subshell thresholds.

Parameters:

Name Type Description Default
E float

Incident photon energy in keV.

 required
E_bind float

Binding energy of the target atomic subshell in keV.

 required

Returns:

Type Description
float

Scattered photon energy (E_prime) in keV.
Source code in lib/physics_base.py 
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def sample_bound_compton_energy(E, E_bind):
    """
    Samples the scattered photon energy E_prime using Kahn's method,
    recommitting the sampled kinematics into the energy-transfer fraction scale
    u = (E - E_prime) / E bounded by strict atomic subshell thresholds.

    Parameters
    ----------
    E : float
        Incident photon energy in keV.
    E_bind : float
        Binding energy of the target atomic subshell in keV.

    Returns
    -------
    float
        Scattered photon energy (E_prime) in keV.
    """
    k = E / me
    α1 = 1.0 + 2.0 * k
    α2 = 1.0 + 1.0 / α1

    # Boundaries mapped directly from your Klein_Nishina threshold limits
    u_min = E_bind / E
    u_max = 2.0 * k / α1

    # Physical gate: if the photon energy cannot overcome binding, return unmodified
    if u_min >= u_max or E <= E_bind:
        return E

    while True:
        r1, r2, r3 = np.random.random(3)

        # Branch A: Sampling from the 1/u envelope
        if r1 < (1.0 / α2):
            u_kahn = 1.0 / (1.0 + 2.0 * k * r2)
            if r3 <= 0.5 * (u_kahn**2 + (1.0 - (1.0 - u_kahn) / k)**2):
                u_transfer = 1.0 - u_kahn
                if u_min <= u_transfer <= u_max:
                    return u_kahn * E
        # Branch B: Sampling from the flat envelope
        else:
            u_kahn = 1.0 - 2.0 * k * r2 / α1
            if r3 <= (1.0 + u_kahn**2) / (1.0 + (1.0 - u_kahn) / k):
                u_transfer = 1.0 - u_kahn
                if u_min <= u_transfer <= u_max:
                    return u_kahn * E

σ_KN_bound_shell(E, E_bind, n_elec=1)

Calculates the integrated cross-section (in barns) for a specific bound shell using YOUR fast analytical antiderivative with strict energy threshold check.

Parameters:

Name Type Description Default
E float

Incident photon energy in keV.

 required
E_bind float

Binding energy of the target atomic subshell in keV.

 required
n_elec

Number of electrons as multiplication factor.

 1

Returns:

Type Description
float

Scattering x-section in barns.
Source code in lib/physics_base.py 
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
def σ_KN_bound_shell(E, E_bind, n_elec=1):
    """
    Calculates the integrated cross-section (in barns) for a specific bound shell
    using YOUR fast analytical antiderivative with strict energy threshold check.

    Parameters
    ----------
    E : float
        Incident photon energy in keV.
    E_bind : float
        Binding energy of the target atomic subshell in keV.
    n_elec: float
        Number of electrons as multiplication factor.

    Returns
    -------
    float
        Scattering x-section in barns.

    """
    k = E / me
    u_min = E_bind / E
    u_max = 2 * k / (1 + 2 * k)

    # If photon energy is below binding threshold or physically impossible
    if u_min >= u_max or E <= E_bind or k < 1e-4:
        return 0.0

    # Lightning fast analytical boundary subtraction!
    res = adrv_bound(u_max, k) - adrv_bound(u_min, k)
    return n_elec * res

σ_KN_free(E)

Calculates the exact analytical total Klein-Nishina cross-section for a single FREE electron (in barns) using the general QED expression. Safe and stable near zero energy limit (converges to Thomson cross-section).

Parameters:

Name Type Description Default
E float

Incident photon energy in keV.

 required

Returns:

Type Description
float

Scattering x-section in barns.
Source code in lib/physics_base.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def σ_KN_free(E):
    """
    Calculates the exact analytical total Klein-Nishina cross-section 
    for a single FREE electron (in barns) using the general QED expression.
    Safe and stable near zero energy limit (converges to Thomson cross-section).

    Parameters
    ----------
    E : float
        Incident photon energy in keV.

    Returns
    -------
    float
        Scattering x-section in barns.

    """
    k = E / me
    if k < 1e-4:
        return se * (8.0 / 3.0)  # Pure classical Thomson limit: ~0.6652 barn

    # Standard integrated Klein-Nishina formulation
    part1 = ((1 + k) / (k**2)) * ((2 * (1 + k)) / (1 + 2*k) - np.log(1 + 2*k) / k)
    part2 = np.log(1 + 2*k) / (2 * k)
    part3 = (1 + 3*k) / ((1 + 2*k)**2)

    return 2 * se * (part1 + part2 - part3)

Target Atom Management

The target parameters and cubic spline interpolators are managed by the Atom configuration engine.

atom_prepare

Database Builder and Physics Compiler Module for GEIM Simulation.

Parses processed EPICS2023 raw data structures, builds high-precision CubicSpline interpolators, calculates atomic relaxation tables, computes calibrated multi-shell Compton cross-sections using physics_base analytics, and dumps everything into a single performance-optimized pickle database.

build_database(source_pckl='EPICS2023.pckl', output_pckl='geim_data.pckl')

Loads raw EPICS parser output, compiles CubicSplines and relaxation tables, integrates physics_base formulations, and exports the final simulation database.

Source code in atom_prepare.py 
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
def build_database(source_pckl='EPICS2023.pckl', output_pckl='geim_data.pckl'):
    """
    Loads raw EPICS parser output, compiles CubicSplines and relaxation tables,
    integrates physics_base formulations, and exports the final simulation database.
    """
    print(f"Loading raw parser data from {source_pckl}...")
    with open(source_pckl, 'rb') as fi:
        Z, M, DATA = pickle.load(fi)

    print("Compiling atomic subshell configurations...")
    subshells = DATA[0][0][91][912][0]['shell']
    shell_labels = ['1s¹⁄₂','2s¹⁄₂','2p¹⁄₂','2p³⁄₂','3s¹⁄₂','3p¹⁄₂','3p³⁄₂','3d³⁄₂','3d⁵⁄₂','4s¹⁄₂','4p¹⁄₂','4p³⁄₂']

    # Pack basic configurations
    shell_names = dict(zip(subshells, shell_labels))
    n_elec = dict(zip(subshells, DATA[0][0][91][912][0][' num_electrons']))
    e_bind = dict(zip(subshells, DATA[0][0][91][913][0]['binding_energy']))
    e_loc = dict(zip(subshells, DATA[0][0][92][935][0]['energy to atom']))

    print("Compiling photon integrated interaction splines...")
    t_coh = DATA[7][0][71][0][0]
    t_inc = DATA[7][0][72][0][0]
    t_pe  = DATA[7][0][73][0][0]

    splines = {
        'COHERENT': CubicSpline(t_coh['Eγ'], t_coh['σ_coherent'], extrapolate=False),
        'INCOHERENT': CubicSpline(t_inc['Eγ'], t_inc['σ_incoherent'], extrapolate=False),
        'PHOTOELECT': CubicSpline(t_pe['Eγ'], t_pe['σ_photoelect'], extrapolate=False),

        'FormFactor': CubicSpline(DATA[7][0][93][941][0]['x'], DATA[7][0][93][941][0]['form_factor '], extrapolate=True),
        'Scatt_Func': CubicSpline(DATA[7][0][93][942][0]['x'], DATA[7][0][93][942][0]['scatter_func'], extrapolate=True),
        'ImAnScFact': CubicSpline(DATA[7][0][93][943][0]['Eγ'], DATA[7][0][93][943][0]['Im_an_sc_fac'], extrapolate=True),
        'ReAnScFact': CubicSpline(DATA[7][0][93][944][0]['Eγ'], DATA[7][0][93][944][0]['Re_an_sc_fac'], extrapolate=True)
    }

    print("Compiling electron integrated interaction splines...")
    splines.update({
        'BREMSSTRAH': CubicSpline(DATA[9][0][82][0][0]['Ee'], DATA[9][0][82][0][0]['σ_bremsstrah'], extrapolate=True),
        'AvEeBremss': CubicSpline(DATA[9][9][82][10][0]['Ee'], DATA[9][9][82][10][0]['<Ee_bremsst>'], extrapolate=True),
        'AvEγBremss': CubicSpline(DATA[9][7][82][10][0]['Ee'], DATA[9][7][82][10][0]['<Eγ_bremsst>'], extrapolate=True),

        'ELASTTRANS': CubicSpline(DATA[9][0][7][0][0]['Ee'], DATA[9][0][7][0][0]['σ_elastic_trans'], extrapolate=True),
        'ELALASCATT': CubicSpline(DATA[9][0][8][0][0]['Ee'], DATA[9][0][8][0][0]['σ_LA_elastic_sc'], extrapolate=True),
        'ELASTSCATT': CubicSpline(DATA[9][0][10][0][0]['Ee'], DATA[9][0][10][0][0]['σ_elastic_scatt'], extrapolate=True)
    })

    print("Compiling shell-specific interaction splines...")
    phel_shell, ionz_shell, reco_ebind = {}, {}, {}
    for ss in subshells:
        t_phel = DATA[7][0][73][0][ss]
        phel_shell[ss] = CubicSpline(t_phel['Eγ'][1:], t_phel['σ_photoelect'][1:], extrapolate=False)

        t_ionz = DATA[9][0][81][0][ss]
        ionz_shell[ss] = CubicSpline(t_ionz['Ee'][1:], t_ionz['σ_ionization'][1:], extrapolate=False)

        t_reco = DATA[9][19][81][10][ss]
        reco_ebind[ss] = CubicSpline(t_reco['Ee'], t_reco['<Ee_recoil>'] / e_bind[ss], extrapolate=False)

    splines['PHEL_SHELL'] = phel_shell
    splines['IONZ_SHELL'] = ionz_shell
    splines['RECO_EBIND'] = reco_ebind

    # --- GEIM MULTI-SHELL COMPTON HYBRID CALIBRATION ---
    print("Compiling calibrated shell-specific Compton cross-section splines...")
    compton_grid_Eγ = t_inc['Eγ']  # Keep the exact, unmodified original EPICS2023 grid

    comp_shell = {}
    for ss in subshells:
        Eb = e_bind[ss]
        Ne = n_elec[ss]

        # Allocate the calibrated array over the exact master grid size
        shell_σ_calibrated = np.zeros(len(compton_grid_Eγ))

        for idx, Eγ in enumerate(compton_grid_Eγ):
            # If the photon energy is below or equal to the binding threshold, 
            # the cross-section is strictly physical zero.
            if Eγ <= Eb:
                shell_σ_calibrated[idx] = 0.0
                continue

            σ_epics_total = t_inc['σ_incoherent'][idx]

            # Sum up theoretical bound-shell cross-sections for the denominator
            Σ_step = 0.0
            for current_ss in subshells:
                Σ_step += physics_base.σ_KN_bound_shell(Eγ, e_bind[current_ss], n_elec[current_ss])

            # Compute analytical weight for the current target subshell
            σ_target_theory = physics_base.σ_KN_bound_shell(Eγ, Eb, Ne)

            if Σ_step > 1e-6:
                shell_σ_calibrated[idx] = σ_epics_total * (σ_target_theory / Σ_step)
            else:
                shell_σ_calibrated[idx] = 0.0

        # Construct the CubicSpline using the unbroken master energetic grid.
        # This guarantees alignment and synchronization across all atomic processes.
        comp_shell[ss] = CubicSpline(compton_grid_Eγ, shell_σ_calibrated, extrapolate=False)

    splines['COMP_SHELL'] = comp_shell


    print("Computing atomic relaxation cascade tables...")
    relaxation_table = calculate_relaxation_table(subshells, DATA)

    # Pack everything into the final compiled dict structure
    compiled_db = {
        'Z': Z,
        'A_raw': M,
        'subshells': subshells,
        'shell_names': shell_names,
        'n_elec': n_elec,
        'e_bind': e_bind,
        'e_loc': e_loc,
        'splines': splines,
        'relaxation': relaxation_table
    }

    print(f"Saving compiled performance database to {output_pckl}...")
    with open(output_pckl, 'wb') as fo:
        pickle.dump(compiled_db, fo, protocol=pickle.HIGHEST_PROTOCOL)
    print("Database built successfully!\n")

calculate_relaxation_table(subshells, raw_data)

Computes cumulative atomic relaxation probability tables for cascade logic.

Splits transition data into 2-step Monte Carlo lookups aligned with EADL: 1. Select the intermediate vacancy subshell. 2. Decide outcome (Radiative photon or Non-radiative Auger electron).

Source code in atom_prepare.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def calculate_relaxation_table(subshells, raw_data):
    """
    Computes cumulative atomic relaxation probability tables for cascade logic.

    Splits transition data into 2-step Monte Carlo lookups aligned with EADL:
    1. Select the intermediate vacancy subshell.
    2. Decide outcome (Radiative photon or Non-radiative Auger electron).
    """
    S, L = list(subshells), len(subshells)
    T = {}

    # Process only inner shells capable of relaxation (skip valence/outermost)
    for V0 in S[:-3]:
        # Access descriptors exactly matching the data structure
        Pp = raw_data[0][7][92][931][V0]  # Radiative transitions data
        Pe = raw_data[0][9][92][932][V0]  # Non-radiative transitions data

        P2D = np.zeros((L - 1, L - 2))

        # Fill radiative transitions
        for j, k in zip(Pp['efrom'], Pp['chance']):
            ind = min(S.index(j), L - 3)
            P2D[0][ind] += k  # Integrated probability row
            P2D[1][ind] += k  # Photon emission row

        # Fill non-radiative transitions
        for i, j, k in zip(Pe['efrom'], Pe['pe_from'], Pe['chance']):
            indr = min(1 + S.index(i), L - 2)
            indc = min(S.index(j), L - 3)
            P2D[0][indc] += k  # Integrated probability row
            P2D[indr][indc] += k

        # Normalize transition channels inside columns
        for j in range(L - 2):
            col_sum = np.sum(P2D[1:, j])
            if col_sum > 1e-6:
                for i in range(1, L - 2):
                    P2D[i][j] /= col_sum

        # Convert standardized matrix to clean dict lookup tables for Monte Carlo
        sub_table = {'V1': [], 'P1': []}
        for j in range(len(P2D[0])):
            if P2D[0][j] > 1e-5:
                V1 = subshells[j]
                sub_table['V1'].append(V1)
                sub_table['P1'].append(P2D[0][j])
                sub_table[V1] = {'V2': [], 'P2': []}
                for i in range(1, len(P2D[:, j])):
                    if P2D[i][j] > 1e-5:
                        V2 = 0 if i == 1 else subshells[i - 1]
                        sub_table[V1]['V2'].append(V2)
                        sub_table[V1]['P2'].append(P2D[i][j])

        # Normalize outer primary vacancy probabilities
        p1_sum = sum(sub_table['P1'])
        if p1_sum > 1e-10:
            sub_table['P1'] = [e / p1_sum for e in sub_table['P1']]

        # Normalize internal transition outcomes
        for k in sub_table['V1']:
            p2_sum = sum(sub_table[k]['P2'])
            if p2_sum > 1e-10:
                sub_table[k]['P2'] = [e / p2_sum for e in sub_table[k]['P2']]

        T[V0] = sub_table
    return T