add changes

README.md

@@ -12,7 +12,7 @@ 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)
+track structure theory](docs/RadMeas_paper.pdf)
 by Jeppe Brage Christensen.
 
 > **Note** - due to its dependency on

calibration/lib_optimiser.py

@@ -1207,20 +1207,27 @@ def estimate_parameter_uncertainties(
         aggregation=aggregation,
     )
 
+    # The Hessian must reflect data constraint only — smoothness_weight=0 here.
+    # The regulariser is a prior that helps the optimiser converge but its
+    # curvature is artificial: it makes unconstrained high-dose anchors look
+    # well-identified, producing falsely *small* uncertainties for the anchors
+    # least constrained by the data.  Dropping it from the Hessian lets those
+    # anchors correctly surface as large-uncertainty parameters.
     objective_fn = partial(
         _objective_direct_logv,
         etch_model=etch_model,
         models_dict=models_dict,
         log_v_bulk=log_v_bulk,
         log_v_max=log_v_max,
         n_free=n_free,
-        smoothness_weight=smoothness_weight,
+        smoothness_weight=0.0,
         fit_debris_damping=fit_debris_damping,
         **cost_kwargs,
     )
 
+    # Evaluate cost at optimum using the data-only objective for scaling.
     f0 = objective_fn(x0)
-    print(f"Cost at optimum: {f0:.6e}")
+    print(f"Cost at optimum (data-only, no regularisation): {f0:.6e}")
     print(f"Parameters: {n_params}")
 
     # --- Compute Hessian via central finite differences ---

calibration/run_Vd_calibration.py

@@ -37,7 +37,7 @@ def starting_guess_v_um_h(
     Interpolate anchor velocities from a previously saved optimized model.
 
     Loads the saved model JSON, then uses PCHIP interpolation to map its
-    V(D) curve onto the requested anchor doses.  Values are clipped to
+    v(d) curve onto the requested anchor doses.  Values are clipped to
     [V_bulk, max_v_um].
     """
     saved_model = load_etchrate_model(CALIBRATED_MODEL_NAME)
@@ -70,30 +70,33 @@ def starting_guess_v_um_h(
 final_fit = True
 
 if final_fit:
-    # Final-fit profile aligned with benchmark settings
-    de_popsize = 60
-    de_maxiter = 200
+    # Most thorough fit profile: broad DE search with extended convergence.
+    de_popsize = 70
+    de_maxiter = 300
     de_tol = 1e-6
     de_mutation = (0.5, 1.4)
     de_recombination = 0.9
     de_polish_method = "Powell"
     de_smoothness_weight = 0.001
-    fit_debris_damping = True
+    fit_debris_damping = True  # 5%-capped damping is now identifiable
 else:
     # Dev profile: faster iterations for quick checks
-    de_popsize = 15
-    de_maxiter = 30
+    de_popsize = 30
+    de_maxiter = 50
     de_tol = 1e-5
     de_mutation = (0.5, 1.0)
     de_recombination = 0.7
     de_polish_method = "Powell"
     de_smoothness_weight = 0.0
-    fit_debris_damping = False
+    fit_debris_damping = True  # 5%-capped damping is now identifiable
 
-# ---- Debris damping (co-optimized with V(D)) ----
+# ---- Debris damping ----
+# Not co-optimized: shape data alone cannot identify alpha/beta (optimizer
+# always escapes to alpha=max, beta=max, i.e. damping disabled).  Set fixed
+# values here if debris correction is desired; leave as None to disable.
 # eta = 1 / (1 + (aspect_ratio / alpha)**beta)
-debris_alpha_bounds_log10 = (0.1, 4.0)  # alpha in [0.1, 1e4]
-debris_beta_bounds = (0.3, 4.0)
+debris_alpha_bounds_log10 = (0.1, 3)
+debris_beta_bounds = (1.0, 6.0)
 
 # ---- Anchor setup ----
 # Only dose positions matter; DE finds the velocities automatically.
@@ -165,7 +168,7 @@ def starting_guess_v_um_h(
 )
 print(f"Initial cost: {J_initial:.6e}")
 
-# ---- Single-pass: co-optimize V(D) + debris damping ----
+# ---- Single-pass: co-optimize v(d) + debris damping ----
 result = optimize_etch_model(
     etch_model,
     models_dict,

examples/plot_debris_comparison.py

@@ -0,0 +1,49 @@
+# %%
+"""
+Track shapes with and without debris damping.
+
+Simulates a 7Li ion at low kinetic energy; a regime where
+the narrow, deep pit geometry makes debris-transport effects visible.
+Contours are shown for three etching times; each time gets a pair of
+curves: with debris (solid) and without (dashed).
+"""
+
+import matplotlib.pyplot as plt
+from tracketch import TrackSimulator, load_etchrate_model
+
+kinetic_energy_MeV_u = 1.0  # low energy to accentuate debris effect
+particle_name = "12C"
+
+# Load the calibrated model (includes debris parameters from calibration).
+sim = TrackSimulator(
+    particle_name=particle_name,
+    start_energy_MeV_u=kinetic_energy_MeV_u,
+    etch_model=load_etchrate_model("Doerschel_etching"),
+    r_max_um=1,
+)
+
+etch_times_h = [0.5]
+
+fig, ax = plt.subplots()
+for idx, t_h in enumerate(etch_times_h):
+    color = plt.cm.viridis(idx / len(etch_times_h))
+
+    # --- with debris (solid) ---
+    r, z = sim.get_iso_time_contour(t_h)
+    ax.plot(r, z, color=color, linestyle="-", label=f"{t_h} h  (debris)")
+    ax.plot(-r, z, color=color, linestyle="-")
+
+    # --- without debris (dashed) ---
+    r0, z0 = sim.get_iso_time_contour_nodebris(t_h)
+    ax.plot(r0, z0, color=color, linestyle="--", label=f"{t_h} h  (no debris)")
+    ax.plot(-r0, z0, color=color, linestyle="--")
+
+ax.grid(True, which="both", ls=":", color="gray", alpha=0.5)
+ax.set_xlabel("r / µm")
+ax.set_ylabel("z / µm")
+ax.set_title(f"{particle_name} @ {kinetic_energy_MeV_u} MeV/u\ndebris damping effect")
+ax.invert_yaxis()
+ax.legend(title="Etch time", bbox_to_anchor=(1.05, 1), loc="upper left")
+# ax.set_aspect("equal")
+plt.tight_layout()
+plt.show()

tracketch/etching/models/Doerschel_etching.json

@@ -2,56 +2,65 @@
   "name": "Doerschel_etching",
   "V_bulk_um_h": 1.73,
   "V_max_um_h": 100.0,
-  "V_max_uncertainty_um_h": 8.322281692861123,
+  "V_max_uncertainty_um_h": 5.115263265889307,
   "anchor_doses_Gy": [
     0.001,
     100.0,
+    1000.0,
     10000.0,
-    31622.776601683792,
+    21544.346900318822,
+    46415.88833612782,
     100000.0,
-    237137.37056616554,
-    562341.3251903491,
-    1333521.432163324,
-    3162277.6601683795,
-    7498942.093324558,
-    17782794.100389227,
-    42169650.342858225,
+    215443.46900318822,
+    464158.8833612782,
+    1000000.0,
+    2154434.6900318824,
+    4641588.833612773,
+    10000000.0,
+    21544346.900318824,
+    46415888.33612773,
     100000000.0,
     1000000000.0
   ],
   "anchor_velocities_um_h": [
     1.73,
-    1.7722745958779738,
-    1.8055135419252157,
-    1.934337911217451,
-    1.9429304969587473,
-    2.18881199713672,
-    2.199057329524819,
-    2.3626080097098825,
-    2.439817422350564,
-    2.5871586093465084,
-    2.954025314321479,
-    5.75319577937238,
-    10.089309684704848,
-    88.30593825964475
+    1.73016044103984,
+    1.7301604430845687,
+    1.7301604611795156,
+    1.7301604923124645,
+    1.7302141139654423,
+    1.7302654376187288,
+    1.7302654387840601,
+    1.8531846492207886,
+    1.8614834215674412,
+    2.063053727909972,
+    2.3429665835898352,
+    2.512707884627541,
+    3.7844470548250797,
+    5.796747249913689,
+    10.359467365079638,
+    99.99149095394164
   ],
   "anchor_velocity_uncertainties_um_h": [
     0.0346,
-    0.4670112122185652,
-    5.088253928348727,
-    6.265837626654851,
-    6.199807196196782,
-    3.378839880422098,
-    3.1014721779715413,
-    0.04725216019419765,
-    0.23512886353977802,
-    0.05174317218693017,
-    0.47407857540789694,
-    0.1150639155874476,
-    0.40138504004371267,
-    8.322281692861123
+    0.0975321268343697,
+    0.0975307026261057,
+    0.0975671376452783,
+    0.0977062356721335,
+    0.0825104217474753,
+    0.1085008616940364,
+    0.1085010421146749,
+    0.13706369298441577,
+    0.137229668431348825,
+    0.14126107455819944,
+    0.15554759796050469,
+    0.18154399948085776,
+    0.2756889410965016,
+    0.11593494499827378,
+    0.20718934730159277,
+    3.115263265889307
   ],
   "extrapolation_mode": "pchip",
-  "debris_alpha": 699.833485055439,
-  "debris_beta": 1.5103813534152948
+  "debris_alpha": 50.117704545098704,
+  "debris_beta": 2.7656543206968776
 }