Examples
The examples/ directory contains runnable scripts that demonstrate the main
features of tracketch. Each section below shows the code and explains what it
does.
Simulating track shapes
The most common use-case: create a TrackSimulator, then
extract the etched track contour at several etching durations.
"""
Track shapes for different ions and etching times.
Simulates the full pipeline: dose map -> etch rate -> arrival time -> track
contour. Shows how track shape evolves with etching time for a 270 MeV/u
carbon ion in CR39.
If the loaded etch model carries anchor velocity uncertainties, a 1-sigma
shaded band is added around each contour.
"""
import matplotlib.pyplot as plt
from tracketch import TrackSimulator
sim = TrackSimulator(
particle_name="12C",
start_energy_MeV_u=270.0,
r_max_um=5,
)
# calculate the track contours for different etching times
etch_times_h = [0.5, 1.0, 2.0, 3.0]
fig, ax = plt.subplots()
for idx, t_h in enumerate(etch_times_h):
color = plt.cm.viridis(idx / len(etch_times_h))
unc_result = sim.get_iso_time_contour_with_uncertainty(t_h)
if unc_result is not None:
(r, z), z_band, r_lo, r_hi = unc_result
ax.fill_betweenx(z_band, r_lo, r_hi, alpha=0.25, color=color)
ax.fill_betweenx(z_band, -r_hi, -r_lo, alpha=0.25, color=color)
else:
r, z = sim.get_iso_time_contour(etching_time_h=t_h)
ax.plot(r, z, label=f"{t_h}", color=color)
ax.plot(-r, z, color=color) # mirror for full track profile
ax.set_xlabel("r / um")
ax.set_ylabel("z / um")
ax.set_title("12C @ 270 MeV/u: track shape evolution")
ax.invert_yaxis()
ax.legend(title="Etch time / h")
ax.set_aspect("equal")
The contour is returned as (r, z) arrays in um. Mirroring -r gives
the full axially-symmetric profile.
Hydrogen isotopes – proton, deuteron, triton
Added in version 0.1.
tracketch fully supports all three hydrogen isotopes ("1H", "2H",
"3H") as particle names. The physics is correct for each isotope:
LET is identical at the same kinetic energy per nucleon (MeV/u) – it depends only on charge and velocity, not on mass number.
CSDA range scales linearly with mass number A, because a heavier isotope carries more momentum at the same MeV/u.
Total kinetic energy (MeV) = energy per nucleon x A.
The example below computes track radii for proton, deuteron, and triton across a range of energies, illustrating the range difference at the same MeV/u:
"""
Plotting track radii as a function of ion energy.
This example computes the track radius in CR-39 for a grid of ion energies
and two etching times, using :class:`~tracketch.TrackSimulator`. Results
are collected in a :class:`pandas.DataFrame` and visualised with
``seaborn.lineplot``.
"""
import numpy as np
import itertools
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from tracketch import TrackSimulator
import tqdm
etch_times_h = [5, 10]
particles = ["1H", "2H", "3H"]
energies_MeV_u = np.logspace(-2, 2, num=300)
# Create all unique combinations
combinations = list(itertools.product(particles, energies_MeV_u))
combinations_df = pd.DataFrame(combinations, columns=["particle_name", "energy_MeV_u"])
result_df = pd.DataFrame(
columns=[
"particle_name",
"energy_MeV_u",
"energy_MeV",
"CSDA_um",
"etch_time_h",
"radius_um",
]
)
# iterate through combinations and calculate track radius for each
for idx, row in tqdm.tqdm(combinations_df.iterrows(), total=len(combinations_df)):
particle = row["particle_name"]
energy_MeV_u = row["energy_MeV_u"]
# create simulator
sim = TrackSimulator(
particle_name=particle,
start_energy_MeV_u=energy_MeV_u,
)
for etching_time_h in etch_times_h:
# get results for this etching time
new_row = pd.DataFrame(
{
"Particle": [particle],
"Energy (MeV/u)": [energy_MeV_u],
"Etch time / h": [etching_time_h],
"Energy / MeV": [sim.start_energy_MeV],
"CSDA / um": [sim.CSDA_range_um],
"Track radius / um": [
sim.get_track_radius_um(etch_time_h=etching_time_h)
],
"Track length / um": [
sim.get_track_length_um(etch_time_h=etching_time_h)
],
}
)
result_df = pd.concat([result_df, new_row], ignore_index=True)
# %%
fig, axes = plt.subplots(ncols=2, nrows=2, sharey=True, figsize=(10, 10))
fig.suptitle("Track radius in CR-39 for hydrogen isotopes", y=0.92)
fig.subplots_adjust(wspace=0.05)
(ax1, ax2, ax3, ax4) = axes.flatten()
sns.lineplot(
data=result_df,
x="CSDA / um",
y="Track radius / um",
hue="Particle",
style="Etch time / h",
ax=ax1,
)
sns.lineplot(
data=result_df,
x="Energy (MeV/u)",
y="Track radius / um",
hue="Particle",
style="Etch time / h",
ax=ax2,
)
sns.lineplot(
data=result_df,
x="CSDA / um",
y="Track length / um",
hue="Particle",
style="Etch time / h",
ax=ax3,
)
sns.lineplot(
data=result_df,
x="Energy (MeV/u)",
y="Track length / um",
hue="Particle",
style="Etch time / h",
ax=ax4,
)
for ax in axes.flatten():
ax.grid(alpha=0.3)
ax.set_xscale("log")
plt.show()
# %%
Radial dose distributions
The radial dose distribution (RDD) describes how dose falls off with distance from the ion path. tracketch wraps the Cucinotta model from libamtrack and scales it to the target material.
"""
Radial dose distributions (RDDs) for different ions.
Plots the Cucinotta RDD for protons, helium, and carbon ions at selected
energies. The RDD describes how dose falls off radially from the ion track.
"""
import numpy as np
import matplotlib.pyplot as plt
from tracketch import get_RDD_Gy
r_m = np.logspace(-9, -3, 200)
ions = [
("1H", 100.0),
("4He", 150.0),
("12C", 270.0),
]
fig, ax = plt.subplots()
for particle, E_MeV_u in ions:
dose = get_RDD_Gy(r_m, E_MeV_u, particle, RDD_name="Cucinotta")
ax.loglog(r_m * 1e6, dose, label=f"{particle} @ {E_MeV_u} MeV/u")
ax.set_xlabel("Radius / um")
ax.set_ylabel("Dose / Gy")
ax.legend()
ax.set_title("Radial dose distribution (Cucinotta)")
plt.tight_layout()
plt.show()
LET and CSDA range
Linear energy transfer (LET) and CSDA range can be queried from SRIM tables
or from libamtrack. The plot includes hydrogen isotopes to illustrate that
LET curves for 1H, 2H, and 3H coincide while CSDA ranges differ
by a factor of A.
"""
LET and CSDA range for different ions in CR39.
Plots linear energy transfer (LET) and continuous slowing-down approximation
(CSDA) range as a function of kinetic energy for hydrogen isotopes, helium,
and carbon.
Hydrogen isotopes (``1H``, ``2H``, ``3H``) share the same LET curve because
LET depends only on charge and velocity (MeV/u), not on mass number. Their
CSDA ranges differ by a factor equal to the mass number *A* -- a deuteron
travels twice as far as a proton at the same MeV/u.
"""
import numpy as np
import matplotlib.pyplot as plt
from tracketch import get_LET_keV_um, get_CSDA_um
energies = np.logspace(-1, 2.5, 50)
particles = ["1H", "2H", "3H", "4He", "12C"]
fig, (ax_let, ax_csda) = plt.subplots(1, 2, figsize=(10, 4))
for particle in particles:
LET = [get_LET_keV_um(E, particle, "CR39", source="SRIM") for E in energies]
CSDA = [get_CSDA_um(E, particle, "CR39", source="SRIM") for E in energies]
ax_let.loglog(energies, LET, label=particle)
ax_csda.loglog(energies, CSDA, label=particle)
ax_let.set_xlabel("Energy / (MeV/u)")
ax_let.set_ylabel("LET / (keV/um)")
ax_let.set_title("Linear energy transfer in CR39")
ax_let.legend()
ax_csda.set_xlabel("Energy / (MeV/u)")
ax_csda.set_ylabel("CSDA range / um")
ax_csda.set_title("CSDA range in CR39")
ax_csda.legend()
plt.tight_layout()
plt.show()
Customising the simulation grid
By default TrackSimulator uses a fixed radial and depth
grid. For ions with very long tracks, very short tracks, or when a different
resolution is needed, individual grid parameters can be passed as plain keyword
arguments:
sim = TrackSimulator(
particle_name="12C",
start_energy_MeV_u=10.0,
r_max_um=40, # wider radial grid
z_max_um=80, # deeper depth grid
n_points_z=200,
)
The full example below runs the same ion three times (default, wide, and coarse grids) and overlays the resulting iso-time contours:
"""
Customising the simulation grid.
By default, :class:`~tracketch.TrackSimulator` uses a grid that covers
r in [1e-4, 20] um and z in [0, 40] um. For ions with very long or very
short tracks, or when you need higher/lower resolution, you can override
individual grid parameters as plain keyword arguments:
* ``r_max_um`` -- maximum radial extent of the grid (um)
* ``z_max_um`` -- maximum depth (um); should exceed the CSDA range
* ``r_min_um`` -- inner radial boundary (um)
* ``n_points_r`` -- radial resolution (more points -> slower but finer RDD)
* ``n_points_z`` -- depth resolution
This example runs the same 10 MeV/u carbon ion three times:
1. Default grid
2. Wider, deeper grid (suitable for high-energy heavy ions)
3. Compact low-resolution grid (fast scan / debugging)
and overlays the three iso-time contours to show the effect.
"""
import matplotlib.pyplot as plt
from tracketch import TrackSimulator
particle = "12C"
energy_MeV_u = 10.0
etch_time_h = 5.0
configs = [
dict(label="default", color="C0"),
dict(label="wide grid", color="C1", r_max_um=40, z_max_um=80),
dict(label="coarse grid", color="C2", n_points_r=100, n_points_z=50),
]
fig, ax = plt.subplots(figsize=(5, 6))
for cfg in configs:
label = cfg.pop("label")
color = cfg.pop("color")
sim = TrackSimulator(
particle_name=particle,
start_energy_MeV_u=energy_MeV_u,
**cfg,
)
r, z = sim.get_iso_time_contour(etch_time_h)
ax.plot(r, z, label=label, color=color)
# Mark the bulk-etch depth for the default run only
if label == "default":
z_bulk = sim.etch_model.V_bulk_um_h * etch_time_h
ax.axhline(
z_bulk,
color="gray",
linestyle="--",
linewidth=1,
label=f"bulk surface ({z_bulk:.1f} um)",
)
ax.invert_yaxis()
ax.set_xlabel("r / um")
ax.set_ylabel("z / um")
ax.set_title(f"{particle} @ {energy_MeV_u} MeV/u, t = {etch_time_h} h")
ax.legend()
plt.tight_layout()
plt.show()
# %%
print(f"CSDA range: {sim.CSDA_range_um:.1f} um")
print(
f"Grid r: [{sim._r_lims_um[0]:.1e}, {sim._r_lims_um[1]:.0f}] um "
f"({len(sim._r_grid_um)} points)"
)
print(
f"Grid z: [{sim._z_lims_um[0]:.0f}, {sim._z_lims_um[1]:.0f}] um "
f"({len(sim._z_grid_um)} points)"
)
Comparing arrival-time solvers
tracketch ships two wavefront propagation backends: Fast Marching (FMM) and Dijkstra. This example compares them on a uniform-speed grid.
"""
Arrival time maps using different marching backends.
Compares FMM and Dijkstra (Numba) on a simple uniform-speed grid.
Both should produce similar results; Dijkstra supports non-uniform grids
natively while FMM requires interpolation to a uniform grid.
"""
import numpy as np
import matplotlib.pyplot as plt
from tracketch.wavefront import arrival_time_fmm, arrival_time_dijkstra_fast
# Simple grid with uniform etch rate
r_um = np.linspace(0.01, 5.0, 80)
z_um = np.linspace(0.0, 10.0, 150)
speed = np.ones((len(z_um), len(r_um)))
t_fmm = arrival_time_fmm(r_um, z_um, speed, r_is_logscaled=False)
t_dij = arrival_time_dijkstra_fast(r_um, z_um, speed, r_is_logscaled=False)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
for ax, t_map, title in [(ax1, t_fmm, "FMM"), (ax2, t_dij, "Dijkstra (Numba)")]:
im = ax.pcolormesh(r_um, z_um, t_map, shading="auto")
ax.set_xlabel("r / um")
ax.set_title(title)
fig.colorbar(im, ax=ax, label="Arrival time / hr")
ax1.set_ylabel("z / um")
ax1.invert_yaxis()
plt.tight_layout()
plt.show()