Docs update

README.md

@@ -1,76 +1,77 @@
-[![Documentation](https://img.shields.io/badge/docs-latest-blue.svg)](https://dosimetry.pages.psi.ch/tracketch/)
+[![Documentation](https://img.shields.io/badge/docs-latest-blue.svg)](https://jbrage.github.io/tracketch/)
 [![License](https://img.shields.io/badge/license-MIT-green.svg)](LICENSE)
 [![Python](https://img.shields.io/badge/python-3.10+-blue.svg)](https://www.python.org/)
 [![Tests](https://img.shields.io/badge/tests-pytest-brightgreen.svg)](tests/)
-[![HitCount](https://hits.dwyl.com/jbrage/tracketch.svg?style=flat-square&show=unique)](http://hits.dwyl.com/jbrage/tracketch)
-
-
 
 # tracketch
 
 A Python toolkit for simulating ion tracks in CR-39 plastic nuclear track
-detectors.  Given an ion species, energy, and etching conditions, **tracketch**
-predicts the track countours, diameter, and depth that would be observed under an
-optical microscope after chemical etching.
-
-<p align="center">
-   <img src="docs/icons/logo.png" alt="tracketch logo" width="500">
-</p>
+detectors. Given an ion species, energy, and etching conditions, `tracketch`
+predicts the track countours, diameter and length for a given etching duration.
 
+The work is published in the article
+[A unified model for etched-track formation in CR-39 detectors using amorphous
+track structure theory](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=6469212)
+by Jeppe Brage Christensen.
 
 > **Note** - due to its dependency on
-> [libamtrack](https://libamtrack.github.io/), **tracketch** currently runs on
+> [libamtrack](https://libamtrack.github.io/), `tracketch` currently runs on
 > **Linux only**. Most Windows installations support WSL installations.
 
+<table align="center">
+   <tr>
+      <td bgcolor="white">
+         <img src="docs/icons/logo.png" alt="tracketch logo" width="500">
+      </td>
+   </tr>
+</table>
+
+
+Cite as 
+```bibtex
+@article{christensen2026unified,
+  title={A unified model for etched-track formation in CR-39 detectors using amorphous track structure theory},
+  author={Christensen, Jeppe Brage},
+  journal={SSRN preprint http://dx.doi.org/10.2139/ssrn.6469212},
+  year={2026}
+}
+```
 
-> **Note** - This version only supports simulations of ion impinging perpendicular to the CR39 detector surface.  Support for 3D track geometries and non-axisymmetric effects is planned for future releases.
-> Furthermore, this etch-rate model is calibrated to the detector type and etching conditions used in [Dörschel et al (2003)](https://www.sciencedirect.com/science/article/abs/pii/S1350448702000471). For other detector types and etching conditions, the calibration workflow described in the documentation can be used to re-calibrate the model against experimental data.
+> **Note** - This version only supports simulations of ion impinging
+> perpendicular to the CR39 detector surface. Support for 3D track geometries
+> and non-axisymmetric effects is planned for future releases.
+> Furthermore, this etch-rate model is calibrated to the detector type and
+> etching conditions used in
+> [Dörschel et al (2003)](https://www.sciencedirect.com/science/article/abs/pii/S1350448702000471).
+> For other detector types and etching conditions, the calibration workflow
+> described in the documentation can be used to re-calibrate the model against
+> experimental data.
 
 
 ## Documentation
 
-The documentation is hosted online and includes user guides, API references, and examples:  
-[Online documentation](https://dosimetry.pages.psi.ch/tracketch/)  
+The documentation is hosted online and includes user guides, API references,
+and examples: [Documentation](https://jbrage.github.io/tracketch/)
+A PDF version of the documentation is also available in the `docs/` directory: [tracketch_docs.pdf](docs/tracketch.pdf)
 
 ## How it works
 
 The simulation follows the same physics chain as the real experiment:
 
-1. **Radial dose distribution (RDD)**: the ion deposits energy radially via
-   delta-electrons.  tracketch computes the 2-D dose map along the ion path using
-   models from libamtrack (e.g. Cucinotta) and SRIM stopping-power tables.
+1. **Dose maps**: the ion deposits energy radially via delta-electrons.
+   `tracketch` computes the dose map along the ion path using amorphous track
+   structure models from `libamtrack` (e.g. Cucinotta) and slows down the ion
+   in the CR-39 detector using `SRIM` stopping-power tables.
 2. **Etch-rate model**: a calibrated function converts local dose to local
-   etch velocity *V(D)*.  The bulk (undamaged) material etches at a slower
-   constant rate *V_bulk*.
-3. **Wavefront propagation**: the etchant front is propagated from the
-   detector surface into the bulk using shortest-path algorithms (Dijkstra or
-   Fast Marching).  The result is an *arrival-time map*.
+   etch velocity $v_d$. The bulk (undamaged) material etches at a slower
+   constant rate $v_\text{bulk}$.
+3. **Wavefront propagation**: the etchant front is propagated from the detector
+   surface into the bulk using shortest-path algorithms (Dijkstra or Fast
+   Marching), analogous to the propagation of a wavefront. The result is an
+   *arrival-time map*.
 4. **Track observables**: iso-time contours of the arrival-time map give the
    etched track shape at any desired etching duration.
 
-## Etch-rate model selection
-
-tracketch uses a single default calibration model for all particles:
-`"Doerschel_etching"`.
-
-- The etch-rate model is not selected by `particle_name`.
-- Particle selection affects stopping-power and dose physics, not which
-   etch-rate model file is loaded.
-- Users can still load a different etch-rate model explicitly (for example,
-   a custom re-calibrated model).
-
-```python
-from tracketch import TrackSimulator, load_etchrate_model
-
-# Explicitly load a custom calibration model if needed.
-etch_model = load_etchrate_model("my_custom_model")
-
-sim = TrackSimulator(
-      particle_name="12C",
-      start_energy_MeV_u=270.0,
-      etch_model=etch_model,
-)
-```
 
 ## Installation
 
@@ -79,13 +80,14 @@ sim = TrackSimulator(
 git clone <repo-url>/tracketch
 cd tracketch
 
-# create a virtual environment and install
+# create a virtual environment and install; change "all" to "numba" or "cpp"
+# to only install one backend or to skip testing/docs dependencies
 python3 -m venv .venv
 source .venv/bin/activate
-pip install -e ".[numba]"
+pip install -e ".[all]"
 ```
 
-This installs tracketch with the Numba-accelerated Dijkstra backend, which
+This installs `tracketch` with the Numba-accelerated Dijkstra backend, which
 works out of the box on any Linux machine.
 
 ### Optional: C++ Dijkstra backend
@@ -94,11 +96,10 @@ To also install the faster CPP backend, use
 pip install -e ".[numba,cpp]"
 cd tracketch/wavefront/dijkstra/cpp
 python setup_dijkstra.py build_ext --inplace
-cd -
 ```
 
 ## Minimal working example
-Check out the tracketch examples  in the [`examples/`](examples/) directory.
+Check out the `tracketch` examples  in the [`examples/`](examples/) directory.
 ```python
 import matplotlib.pyplot as plt
 from tracketch import TrackSimulator
@@ -142,7 +143,7 @@ and `symbol` is the element symbol:
 
 The available ions depend on the **stopping-power source**:
 
-- **`stopping_power_source="SRIM"` (default)** - uses tabulated SRIM data.
+- **`stopping_power_source="SRIM"` (default)** - uses tabulated `SRIM` data.
   Only the following particles are supported:
   ```python
   import tracketch
@@ -152,8 +153,8 @@ The available ions depend on the **stopping-power source**:
   Attempting any other particle raises a `ValueError` with a suggestion to
   switch source.
 
-- **`stopping_power_source="libamtrack"`** - uses the libamtrack/PSTAR
-  parametrisation. Accepts any nuclide name recognised by libamtrack, e.g.
+- **`stopping_power_source="libamtrack"`** - uses the `libamtrack`
+  parametrisation. Accepts any nuclide name recognised by `libamtrack`, e.g.
   `"56Fe"`, `"238U"`, `"28Si"`.
 
 ```python
@@ -180,7 +181,7 @@ Pass them as `material_name="CR39"` (default) or `material_name="water"`.
 
 ## Simulation grid
 
-The simulation operates on a 2D cylindrical grid (*r*, *z*) with these defaults:
+The simulation operates on a cylindrical grid (*r*, *z*) with these defaults:
 
 | Parameter | Default | Description |
 |-----------|---------|-------------|
@@ -214,11 +215,16 @@ publication-quality results consider doubling `n_points_r` and `n_points_z`.
 Full API reference and usage guides can be re-built with Sphinx:
 
 ```bash
+pip install -e ".[docs]"
 cd docs
 make html
 ```
-
 Then open `docs/_build/html/index.html` in your browser.
+Or 
+```bash
+make pdlatex
+```
+to build the PDF version of the documentation.
 
 ## Project layout
 
@@ -229,7 +235,8 @@ Then open `docs/_build/html/index.html` in your browser.
 │   │   └── stopping_power/    SRIM data tables and parsing
 │   ├── etching/               EtchRateModel -- dose -> etch velocity
 │   ├── simulation/            TrackSimulator -- full model
-│   └── wavefront/             Track contour calculation through arrival-time solvers (Dijkstra, FMM)
+│   └── wavefront/             Track contour calculation through arrival-time
+│                              solvers (Dijkstra, FMM)
 ├── calibration/               Etch-model calibration against experimental data
 ├── examples/                  Runnable example scripts
 ├── tests/                     pytest test suite

calibration/benchmark_bragg_corrections.py

@@ -44,7 +44,7 @@
     PARTICLES = ["1H", "4He", "7Li", "12C"]
     SEEDS = [11, 22]
     POPSIZE = 50
-    MAXITER = 300
+    MAXITER = 50
     TOL = 1e-6
     OUTPUT_DIR_NAME = "bragg_benchmark"
 else:
@@ -60,18 +60,14 @@
 
 SMOOTHNESS_WEIGHT = 0.001
 
-# More general seeded initialization for V(d):
-# - build a broader log-dose anchor grid
-# - perturb log-velocity increments per seed
-ANCHOR_MIN_DOSE_GY = 1e2
-ANCHOR_MAX_DOSE_GY = 1e9
-# Piecewise anchor placement with highest density in 1e5..1e7 Gy.
-ANCHOR_MID_MIN_DOSE_GY = 1e4
-ANCHOR_DENSE_MIN_DOSE_GY = 1e5
-ANCHOR_DENSE_MAX_DOSE_GY = 1e8
-ANCHOR_DENSITY_WEIGHTS = (2.0, 3.0, 30.0, 5.0)
-# Run one or multiple anchor-count configurations.
-ANCHOR_COUNT_USER_OPTIONS = [40]
+# Match run_Vd_calibration.py anchor buckets exactly.
+# d_Gy = concat([
+#   logspace(2, 4, endpoint=False, num=1),
+#   logspace(4, 5, endpoint=False, num=2),
+#   logspace(5, 8, endpoint=False, num=8),
+#   logspace(8, 9, endpoint=True,  num=2),
+# ])
+ANCHOR_COUNT_USER_OPTIONS = [13]
 INIT_LOG_INCREMENT_NOISE_STD = 0.08
 
 # Run debris damping in separate configurations.
@@ -123,52 +119,24 @@ def _save_track_shape_pngs(models_dict: dict, run_dir: Path) -> list[str]:
 
 
 def _build_user_anchor_dose_grid(anchor_count_user: int) -> np.ndarray:
-    """Build piecewise log-dose anchors with dense spacing in 1e5..1e8 Gy."""
-    if anchor_count_user < 4:
-        return np.logspace(
-            np.log10(ANCHOR_MIN_DOSE_GY),
-            np.log10(ANCHOR_MAX_DOSE_GY),
-            anchor_count_user,
-        )
-
-    dose_edges = np.asarray(
+    """Build the same anchor-dose grid as run_Vd_calibration.py."""
+    user_doses = np.concatenate(
         [
-            ANCHOR_MIN_DOSE_GY,
-            ANCHOR_MID_MIN_DOSE_GY,
-            ANCHOR_DENSE_MIN_DOSE_GY,
-            ANCHOR_DENSE_MAX_DOSE_GY,
-            ANCHOR_MAX_DOSE_GY,
-        ],
-        dtype=float,
+            np.logspace(2, 4, endpoint=False, num=1),
+            np.logspace(4, 5, endpoint=False, num=2),
+            np.logspace(5, 8, endpoint=False, num=8),
+            np.logspace(8, 9, endpoint=True, num=2),
+        ]
     )
-    weights = np.asarray(ANCHOR_DENSITY_WEIGHTS, dtype=float)
+    user_doses = np.unique(user_doses)
+    user_doses = np.sort(user_doses)
 
-    seg_counts = np.maximum(
-        1, np.floor(anchor_count_user * weights / weights.sum()).astype(int)
-    )
-    while int(seg_counts.sum()) < anchor_count_user:
-        seg_counts[int(np.argmax(weights))] += 1
-    while int(seg_counts.sum()) > anchor_count_user:
-        reducible = np.where(seg_counts > 1)[0]
-        if reducible.size == 0:
-            break
-        seg_counts[int(reducible[np.argmin(weights[reducible])])] -= 1
-
-    segments: list[np.ndarray] = []
-    for i, n_pts in enumerate(seg_counts):
-        seg = np.logspace(
-            np.log10(dose_edges[i]),
-            np.log10(dose_edges[i + 1]),
-            int(n_pts),
-            endpoint=(i == len(seg_counts) - 1),
+    if anchor_count_user != len(user_doses):
+        raise ValueError(
+            "anchor_count_user must match run_Vd_calibration anchor count "
+            f"({len(user_doses)}), got {anchor_count_user}"
         )
-        segments.append(seg)
-
-    user_doses = np.sort(np.concatenate(segments))
 
-    # Enforce exact endpoints for reproducibility and requested full range.
-    user_doses[0] = ANCHOR_MIN_DOSE_GY
-    user_doses[-1] = ANCHOR_MAX_DOSE_GY
     return user_doses
 
 

calibration/run_Vd_calibration.py

@@ -71,8 +71,8 @@ def starting_guess_v_um_h(
 
 if final_fit:
     # Final-fit profile aligned with benchmark settings
-    de_popsize = 50
-    de_maxiter = 50
+    de_popsize = 60
+    de_maxiter = 200
     de_tol = 1e-6
     de_mutation = (0.5, 1.4)
     de_recombination = 0.9
@@ -100,10 +100,10 @@ def starting_guess_v_um_h(
 d_Gy = np.concatenate(
     [
         # Match benchmark anchor buckets exactly.
-        np.logspace(2, 4, endpoint=False, num=1),
-        np.logspace(4, 5, endpoint=False, num=2),
-        np.logspace(5, 8, endpoint=False, num=10),
-        np.logspace(8, 9, endpoint=True, num=3),
+        np.logspace(2, 4, endpoint=False, num=2),
+        np.logspace(4, 5, endpoint=False, num=3),
+        np.logspace(5, 8, endpoint=False, num=9),
+        np.logspace(8, 9, endpoint=True, num=2),
     ]
 )
 d_Gy = np.unique(d_Gy)  # Remove duplicates
@@ -208,9 +208,7 @@ def starting_guess_v_um_h(
 # uncertainty: when anchors saturate at V_max, this is the most relevant
 # constraint on the saturation velocity.
 if np.isfinite(etch_model.V_max_um_h):
-    etch_model.V_max_uncertainty_um_h = float(
-        unc["velocity_uncertainties_um_h"][-1]
-    )
+    etch_model.V_max_uncertainty_um_h = float(unc["velocity_uncertainties_um_h"][-1])
 
 # Save results
 etch_model.name = CALIBRATED_MODEL_NAME
@@ -222,27 +220,28 @@ def starting_guess_v_um_h(
     for entry in sub_dict.values():
         entry["simulator"].recalculate_from_current_etch_model()
 
+folder_name = "figures_calibration" if final_fit else "figures_calibration_dev"
+os.makedirs(folder_name, exist_ok=True)
+
 if objective_mode in ["both", "shape"]:
     figs = plot_track_shapes(models_dict["track_shape"], xscale="linear")
     shape_keys = list(models_dict["track_shape"].keys())
     for i, particle in enumerate(shape_keys):
-        figs[i].savefig(f"figures/track_shape_{particle}.png", dpi=300)
+        figs[i].savefig(f"{folder_name}/track_shape_{particle}.png", dpi=300)
 
     # Debris comparison: contours with vs without damping
     debris_figs = plot_track_shapes_debris_comparison(
         models_dict["track_shape"], xscale="linear"
     )
     for i, particle in enumerate(shape_keys):
         debris_figs[i].savefig(
-            f"figures/debris_comparison_{particle}.png",
+            f"{folder_name}/debris_comparison_{particle}.png",
             dpi=300,
         )
 
 if objective_mode in ["both", "length"]:
     figs_by_ion = plot_track_length(models_dict["track_length"])
     for ion, fig in figs_by_ion.items():
-        folder_name = "figures_calibration" if final_fit else "figures_calibration_dev"
-        os.makedirs(folder_name, exist_ok=True)
         fig.savefig(f"{folder_name}/track_length_{ion}.png", dpi=300)
 
     # %%