Skip to content

Commit deeaa5e

Browse files
rustdf: address Codex review of the calibration converter
- tof_to_mz: solve the ModelType-2 quadratic in the numerically stable ("citardauq") form instead of Newton / the naive quadratic formula, which lost precision to cancellation when |c2| is tiny (kept ModelType-1 bit-exact). - from_d_folder: reject unsupported MzCalibration/TimsCalibration ModelType values instead of silently falling through to the model-2 branch. - Document BrukerFormulaConverter as a fixed-calibration converter that ignores the per-call frame_id (consistent with the other SDK-free converters). - Fix stale read_mz_calibration comment that still claimed model 2 ignores c2. Deferred (defensive-only / matches existing style): fallible constructors vs .unwrap() on metadata reads, and guards for non-finite/zero coefficients that real instrument calibration never produces.
1 parent 3bd2389 commit deeaa5e

3 files changed

Lines changed: 42 additions & 9 deletions

File tree

rustdf/src/data/calibration.rs

Lines changed: 22 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -91,10 +91,11 @@ impl MzCalibrator {
9191
/// TOF index -> m/z. Inverts `t = C0 + b*s + c2*s^2 + c3*s^3` for `s`.
9292
pub fn tof_to_mz(&self, tof_index: u32) -> f64 {
9393
let t = self.tof_index_to_time(tof_index as f64);
94-
// Base (linear-in-sqrt) closed form, exact when c2 = c3 = 0.
95-
let mut s = (t - self.c0) / self.b;
96-
if self.c2 != 0.0 || self.c3 != 0.0 {
97-
// Newton refinement for the ModelType-1 cubic curve.
94+
// Linear-in-sqrt estimate; exact when c2 = c3 = 0.
95+
let s0 = (t - self.c0) / self.b;
96+
let s = if self.c3 != 0.0 {
97+
// ModelType-1 cubic: Newton refinement from the linear estimate.
98+
let mut s = s0;
9899
for _ in 0..8 {
99100
let f = self.c0 + self.b * s + self.c2 * s * s + self.c3 * s * s * s - t;
100101
let df = self.b + 2.0 * self.c2 * s + 3.0 * self.c3 * s * s;
@@ -107,7 +108,23 @@ impl MzCalibrator {
107108
break;
108109
}
109110
}
110-
}
111+
s
112+
} else if self.c2 != 0.0 {
113+
// ModelType-2 quadratic `c2*s^2 + b*s + (c0 - t) = 0`, solved in the
114+
// numerically stable ("citardauq") form so the physical root does not
115+
// lose precision to cancellation when |c2| is tiny. b > 0 always, so
116+
// q < 0 and is never zero. Falls back to the linear estimate if the
117+
// discriminant is negative (out-of-range tof).
118+
let disc = self.b * self.b - 4.0 * self.c2 * (self.c0 - t);
119+
if disc < 0.0 {
120+
s0
121+
} else {
122+
let q = -0.5 * (self.b + disc.sqrt());
123+
(self.c0 - t) / q
124+
}
125+
} else {
126+
s0
127+
};
111128
s * s - self.c4
112129
}
113130

rustdf/src/data/handle.rs

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -363,6 +363,13 @@ impl IndexConverter for BrukerLibTimsDataConverter {
363363
/// the same C0/C1/C2 quadratic-in-sqrt(m) curve and reproduces the SDK to a
364364
/// few ppm (the proprietary ModelType-2 C8..C14 fine correction is not
365365
/// modelled).
366+
///
367+
/// This is a **fixed-calibration** converter: it captures one frame's
368+
/// calibration rows + temperatures at build time and applies them to every
369+
/// frame, ignoring the per-call `frame_id` (like the other SDK-free converters
370+
/// `Simple`/`Calibrated`/`Lookup`). Valid because the coefficients are
371+
/// effectively constant across a run; use the SDK-backed `BrukerLib` converter
372+
/// if a run genuinely carries multiple calibration rows.
366373
pub struct BrukerFormulaConverter {
367374
pub mz: crate::data::calibration::MzCalibrator,
368375
pub im: crate::data::calibration::MobilityCalibrator,
@@ -396,6 +403,15 @@ impl BrukerFormulaConverter {
396403
.find(|c| c.id == tims_id)
397404
.ok_or("TimsCalibration row not found")?;
398405

406+
// Reject calibration models this converter does not implement, rather
407+
// than silently mis-computing (m/z model 2 is the default branch).
408+
if mzc.model_type != 1 && mzc.model_type != 2 {
409+
return Err(format!("unsupported MzCalibration ModelType {}", mzc.model_type).into());
410+
}
411+
if tc.model_type != 2 {
412+
return Err(format!("unsupported TimsCalibration ModelType {}", tc.model_type).into());
413+
}
414+
399415
let mz = MzCalibrator::new(
400416
mzc.model_type,
401417
mzc.digitizer_timebase,

rustdf/src/data/meta.rs

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -440,10 +440,10 @@ pub fn read_dia_ms_ms_windows(
440440
/// tof_time = c0 + b*sqrt(m) + c2*m + c3*m^1.5, b = sqrt(1e12 / c1)
441441
/// inverted to give m/z from a TOF index (see `data::calibration::MzCalibrator`).
442442
///
443-
/// For model_type 2 (modern instruments) the base curve reduces to
444-
/// sqrt(mz) = (tof_time - c0) / sqrt(1e12 / c1)
445-
/// (c2/c3 are NOT curve terms in model 2). This base is accurate to ~10 ppm; the
446-
/// remaining error is a proprietary degree-6 correction polynomial in c8..c14.
443+
/// Model 2 (modern instruments) uses the same `c0 + b*sqrt(m) + c2*m` curve
444+
/// (c2 IS a curve term; only the c3 cubic term is model-1 only). That base is
445+
/// accurate to a few ppm; the remaining error is a proprietary correction in
446+
/// c8..c14 (windowed to [c5, c6]) that is not modelled.
447447
pub fn read_mz_calibration(
448448
bruker_d_folder_name: &str,
449449
) -> Result<Vec<MzCalibration>, Box<dyn std::error::Error>> {

0 commit comments

Comments
 (0)