121 lines
4.1 KiB
Python
121 lines
4.1 KiB
Python
import pycatima as catima
|
|
import numpy as np, sys
|
|
import matplotlib.pyplot as plt
|
|
|
|
# --- 1. Constants ---
|
|
P_TORR = 250
|
|
P_CO2 = 3
|
|
if(len(sys.argv)==3):
|
|
P_TORR=int(sys.argv[1])
|
|
P_CO2 = int(sys.argv[2])
|
|
|
|
TEMP_K = 293.15
|
|
R = 8.3144
|
|
MEV2U = 1.0 / 931.494
|
|
|
|
# Gas Density Calculations
|
|
p_pa = P_TORR * 133.322
|
|
molar_density = p_pa / (R * TEMP_K)
|
|
m_he, m_c, m_o= 4.0026, 12.0000, 15.9949
|
|
m_mix_avg = ((1 - P_CO2 / 100) * m_he) + (P_CO2 / 100 * (m_c + 2*m_o))
|
|
rho_g_cm3 = (molar_density * m_mix_avg) / 1e6
|
|
print(f"Gas density at {P_TORR} Torr: {rho_g_cm3:.6e} g/cm^3")
|
|
|
|
# --- 2. Material & Step Setup ---
|
|
material_def = [(m_he, 2, (1 - P_CO2 / 100)), (m_c, 6, P_CO2 / 100), (m_o, 8, 2*P_CO2 / 100)]
|
|
gas_mix = catima.Material(material_def)
|
|
gas_mix.density(rho_g_cm3)
|
|
|
|
# Thickness step settings
|
|
step_mg_cm2 = 0.001 # 1 ug/cm2 steps as per your example -- kept fine for
|
|
# numerical accuracy of the dedx integration itself.
|
|
step_g_cm2 = step_mg_cm2 / 1000.0
|
|
max_steps = 1000000000 # Adjust based on how far you want to track
|
|
|
|
coarse_step_cm = 0.2 # row spacing over most of the track
|
|
fine_step_cm = 0.03 # row spacing near the Bragg peak
|
|
fine_zone_frac = 0.085 # fraction of the *total* range treated as "near the peak"
|
|
|
|
# Set relative integration tolerance (lower = higher precision, slower calculation)
|
|
catima.Config.epsrel = 1e-6
|
|
# Set absolute integration tolerance
|
|
catima.Config.epsabs = 1e-9
|
|
|
|
def generate_lookup(z, mass_u, e_start_mev, label):
|
|
filename = f"Eloss/E_vs_x_{label}.dat"
|
|
header = f"Energy(MeV) \tmg/cm2 \tcm\tEin\tsigE\tsigA\tsigR\tsigX\tcov\ttof\tsp\nStarting Energy: {e_start_mev} MeV"
|
|
|
|
projectile = catima.Projectile(mass_u, z)
|
|
e_u_init = e_start_mev / mass_u
|
|
|
|
# 1. Get exact analytical range directly from CATIMA
|
|
projectile.T(e_u_init)
|
|
total_range_g_cm2 = catima.range(projectile, gas_mix)
|
|
total_range_cm = total_range_g_cm2 / rho_g_cm3
|
|
fine_zone_start_cm = total_range_cm * (1.0 - fine_zone_frac)
|
|
|
|
output = []
|
|
current_dist_cm = 0.0
|
|
|
|
# 2. Step directly at checkpoint resolution
|
|
while current_dist_cm <= total_range_cm * 1.05: # Track slightly past range
|
|
thickness_g_cm2 = current_dist_cm * rho_g_cm3
|
|
|
|
# Configure target thickness for this checkpoint
|
|
gas_mix.density(rho_g_cm3).thickness(thickness_g_cm2)
|
|
projectile.T(e_u_init)
|
|
result = catima.calculate(projectile, gas_mix)
|
|
|
|
# Save checkpoint row
|
|
e_total_out = result.Eout * mass_u
|
|
output.append([
|
|
e_total_out, thickness_g_cm2 * 1000.0, current_dist_cm,
|
|
result.Ein, result.sigma_E, result.sigma_a, result.sigma_r,
|
|
result.sigma_x, result.cov, result.tof, result.sp
|
|
])
|
|
|
|
if result.Eout == 0:
|
|
break
|
|
|
|
# Adaptive spatial step
|
|
step = fine_step_cm if current_dist_cm >= fine_zone_start_cm else coarse_step_cm
|
|
current_dist_cm += step
|
|
|
|
np.savetxt(filename, output, fmt='%.6f', delimiter='\t', header=header)
|
|
print(f"Lookup table created: {filename} ({len(output)} rows)")
|
|
|
|
data = np.array(output)
|
|
energy_mev, dist_cm, sigma_e, tof = data[:, 0], data[:, 2], data[:, 4], data[:, 9]
|
|
|
|
fig, (ax_e, ax_sig, ax_tof) = plt.subplots(3, 1, figsize=(8, 12), sharex=True)
|
|
|
|
ax_e.plot(dist_cm, energy_mev)
|
|
ax_e.fill_between(dist_cm, energy_mev - sigma_e, energy_mev + sigma_e, alpha=0.3)
|
|
ax_e.set_ylabel("Energy (MeV)")
|
|
ax_e.set_title(f"Energy Loss Curve {label.capitalize()}")
|
|
ax_e.grid(True)
|
|
|
|
ax_sig.plot(dist_cm, sigma_e)
|
|
ax_sig.set_ylabel("Energy straggle $\\sigma_E$ (MeV)")
|
|
ax_sig.grid(True)
|
|
|
|
ax_tof.plot(dist_cm, tof)
|
|
ax_tof.set_xlabel("Distance (cm)")
|
|
ax_tof.set_ylabel("Time of flight (ns)")
|
|
ax_tof.grid(True)
|
|
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
|
|
|
|
# --- 3. Run ---
|
|
# Format: generate_lookup(Z, mass_u, E_start_MeV, label)
|
|
generate_lookup(1, 1.0078, 30, "proton")
|
|
generate_lookup(1, 2.01355, 30, "deuteron")
|
|
generate_lookup(2, 4.0026, 50, "alpha")
|
|
generate_lookup(13,26.9815, 80, "aluminum")
|
|
#generate_lookup(9,17.0021, 70, "fluorine")
|
|
#generate_lookup(8,15.9949, 70, "oxygen")
|
|
|