rustdf/data/calibration.rs
1//! SDK-free Bruker timsTOF axis calibration.
2//!
3//! Pure-Rust ports of the calibration formulas Bruker publishes for the TDF
4//! format, so that TOF-index -> m/z and scan -> 1/K0 can be computed without
5//! loading the proprietary `libtimsdata` SDK. The algorithms mirror the
6//! implementation in PAPPSO's GPL library `libpappsomspp`
7//! (`mzcalibrationmodel1.cpp`, `timsframebase.cpp`); the coefficient meanings
8//! are cross-checked against Bruker's own `tims_calibration.py` reference.
9//!
10//! Two independent models are involved, each carrying its own `ModelType`:
11//! * m/z : `MzCalibration` table (this data set: ModelType 2)
12//! * 1/K0 : `TimsCalibration` table (this data set: ModelType 2)
13//!
14//! IMPORTANT (m/z): PAPPSO only implements m/z *ModelType 1*. Modern instruments
15//! write *ModelType 2*, whose first coefficients (C0,C1) still describe the same
16//! `t = C0 + sqrt(1e12/C1)*sqrt(m)` base curve, but which adds a degree-6
17//! correction polynomial (C8..C14) that neither PAPPSO nor this module models.
18//! We therefore reproduce the *base* curve exactly (few-ppm agreement with the
19//! SDK) and, for genuine ModelType-1 data, the full cubic-in-sqrt(m) curve.
20
21/// m/z axis calibration (Bruker "model type 1" base curve + optional cubic).
22///
23/// Flight time from a TOF index: `t = index * timebase + delay`
24/// Calibration curve (time as fn of mass): `t = C0 + b*s + c2*s^2 + c3*s^3`
25/// with `s = sqrt(m + c4)` and `b = sqrt(1e12 / C1_tempcomp)`.
26///
27/// Coefficient meaning (columns of the `MzCalibration` table):
28/// * `timebase` = `DigitizerTimebase` — ns per digitizer sample.
29/// * `delay` = `DigitizerDelay` — fixed time offset (samples) before t0.
30/// * `C0` — constant term of the time/mass curve (~ the t-intercept).
31/// * `C1` — governs the dominant sqrt term; `b = sqrt(1e12 / C1)`.
32/// * `c2` — quadratic term `C2*s^2` of the curve; used by BOTH models.
33/// * `c3` — cubic term `C3*s^3`; ModelType 1 only (in ModelType 2 the C3
34/// column is a duplicate of C0 and is dropped).
35/// * `c4` — "reduced mass" shift m0 (`x = m - m0`); patent US7,851,746.
36/// ModelType 1 only (in ModelType 2 the C4 column is a duplicate of
37/// C2 and is dropped, exactly like C3).
38/// Temperature compensation (`T1/T2` = reference temps in `MzCalibration`,
39/// `dC1/dC2` its sensitivities, `T1f/T2f` = per-frame `Frames.T1/Frames.T2`):
40/// `tc = 1 + (dC1*(T1-T1f) + dC2*(T2-T2f)) / 1e6`, applied as `C1 *= tc`.
41#[derive(Debug, Clone)]
42pub struct MzCalibrator {
43 pub timebase: f64,
44 pub delay: f64,
45 pub c0: f64,
46 pub b: f64, // sqrt(1e12 / (C1 * tc))
47 pub c2: f64,
48 pub c3: f64,
49 pub c4: f64,
50}
51
52impl MzCalibrator {
53 /// Build a calibrator from raw `MzCalibration` columns + per-frame temps.
54 ///
55 /// `model_type` selects whether the C2/C3 curve terms are honoured (type 1)
56 /// or zeroed (type 2, base curve only).
57 #[allow(clippy::too_many_arguments)]
58 pub fn new(
59 model_type: i64,
60 timebase: f64,
61 delay: f64,
62 t1_ref: f64,
63 t2_ref: f64,
64 dc1: f64,
65 dc2: f64,
66 c0: f64,
67 c1: f64,
68 c2: f64,
69 c3: f64,
70 c4: f64,
71 t1_frame: f64,
72 t2_frame: f64,
73 ) -> Self {
74 let tc = 1.0 + (dc1 * (t1_ref - t1_frame) + dc2 * (t2_ref - t2_frame)) / 1.0e6;
75 let b = (1.0e12 / (c1 * tc)).sqrt();
76 // Both models share the quadratic-in-sqrt(m) curve `t = C0 + b*s + C2*s^2`
77 // (empirically: a0->C0, a1->sqrt(1e12/C1), a2->C2). The cubic `C3*s^3`
78 // term is real only for ModelType 1; in ModelType 2 the C3 column is a
79 // duplicate of C0 and must be dropped. ModelType 2 additionally carries a
80 // C8..C14 fine correction (~few ppm, worst at low m/z) that is NOT an
81 // additive polynomial in m/z and is left unmodelled here.
82 let c2 = c2 / tc;
83 // ModelType 2 reuses the C3 *and* C4 columns as duplicates of C0/C2
84 // (verified bit-for-bit on real files: C3 == C0, C4 == C2), so neither is
85 // the cubic term nor the reduced-mass shift there and both must be
86 // dropped. Subtracting the duplicated C4 as if it were m0 costs ~1.4 ppm
87 // mean / 6 ppm max against the SDK on ModelType-2 data.
88 let (c3, c4) = if model_type == 1 { (c3, c4) } else { (0.0, 0.0) };
89 Self { timebase, delay, c0, b, c2, c3, c4 }
90 }
91
92 /// Build a calibrator straight from an `MzCalibration` row plus the
93 /// per-frame digitizer temperatures (`Frames.T1`, `Frames.T2`).
94 pub fn from_calibration(
95 cal: &crate::data::meta::MzCalibration,
96 t1_frame: f64,
97 t2_frame: f64,
98 ) -> Self {
99 Self::new(
100 cal.model_type,
101 cal.digitizer_timebase,
102 cal.digitizer_delay,
103 cal.t1,
104 cal.t2,
105 cal.dc1,
106 cal.dc2,
107 cal.c0,
108 cal.c1,
109 cal.c2,
110 cal.c3,
111 cal.c4,
112 t1_frame,
113 t2_frame,
114 )
115 }
116
117 /// Flight time (digitizer units) for a TOF index.
118 #[inline]
119 fn tof_index_to_time(&self, tof_index: f64) -> f64 {
120 tof_index * self.timebase + self.delay
121 }
122
123 /// TOF index -> m/z. Inverts `t = C0 + b*s + c2*s^2 + c3*s^3` for `s`.
124 pub fn tof_to_mz(&self, tof_index: u32) -> f64 {
125 let t = self.tof_index_to_time(tof_index as f64);
126 // Linear-in-sqrt estimate; exact when c2 = c3 = 0.
127 let s0 = (t - self.c0) / self.b;
128 let s = if self.c3 != 0.0 {
129 // ModelType-1 cubic: Newton refinement from the linear estimate.
130 let mut s = s0;
131 for _ in 0..8 {
132 let f = self.c0 + self.b * s + self.c2 * s * s + self.c3 * s * s * s - t;
133 let df = self.b + 2.0 * self.c2 * s + 3.0 * self.c3 * s * s;
134 if df == 0.0 {
135 break;
136 }
137 let step = f / df;
138 s -= step;
139 if step.abs() < 1e-12 {
140 break;
141 }
142 }
143 s
144 } else if self.c2 != 0.0 {
145 // ModelType-2 quadratic `c2*s^2 + b*s + (c0 - t) = 0`, solved in the
146 // numerically stable ("citardauq") form so the physical root does not
147 // lose precision to cancellation when |c2| is tiny. b > 0 always, so
148 // q < 0 and is never zero. Falls back to the linear estimate if the
149 // discriminant is negative (out-of-range tof).
150 let disc = self.b * self.b - 4.0 * self.c2 * (self.c0 - t);
151 if disc < 0.0 {
152 s0
153 } else {
154 let q = -0.5 * (self.b + disc.sqrt());
155 (self.c0 - t) / q
156 }
157 } else {
158 s0
159 };
160 s * s - self.c4
161 }
162
163 /// m/z -> TOF index (forward direction, always closed form).
164 pub fn mz_to_tof(&self, mz: f64) -> u32 {
165 let s = (mz + self.c4).max(0.0).sqrt();
166 let t = self.c0 + self.b * s + self.c2 * s * s + self.c3 * s * s * s;
167 (((t - self.delay) / self.timebase).round()).max(0.0) as u32
168 }
169}
170
171/// Ion-mobility axis calibration (Bruker "model type 2", the only TIMS model).
172///
173/// Two steps, both exact ports of PAPPSO `timsframebase.cpp`:
174/// 1. scan -> trapping voltage: `V = dv_start + slope*(scan - ttrans - ndelay)`
175/// with `slope = (dv_end - dv_start) / ncycles`. V must lie in [vmin,vmax].
176/// 2. voltage -> inverse mobility: `1/K0 = 1 / (C0m + C1m / V)`.
177///
178/// Coefficient meaning (columns of the `TimsCalibration` table, ModelType 2):
179/// * `C0` = `ndelay` — scan offset (delay), subtracted before scaling.
180/// * `C1` = `ncycles` — number of TIMS cycles; sets the voltage-vs-scan slope.
181/// * `C2` = `dv_start` — trapping voltage at the start of the ramp.
182/// * `C3` = `dv_end` — trapping voltage at the end of the ramp.
183/// * `C4` = `ttrans` — transit time in cycles, subtracted before scaling.
184/// * `C5` — unused by the mobility formula (polynomial grade flag).
185/// * `C6` = `C0m` — additive constant of the mobility reciprocal.
186/// * `C7` = `C1m` — voltage-scaled term of the mobility reciprocal.
187/// * `C8` = `vmin` — lower voltage validity bound.
188/// * `C9` = `vmax` — upper voltage validity bound.
189#[derive(Debug, Clone)]
190pub struct MobilityCalibrator {
191 pub ndelay: f64,
192 pub dv_start: f64,
193 pub ttrans: f64,
194 pub c0m: f64,
195 pub c1m: f64,
196 pub vmin: f64,
197 pub vmax: f64,
198 pub slope: f64,
199}
200
201impl MobilityCalibrator {
202 /// Build from raw `TimsCalibration` C0..C9 (ModelType must be 2).
203 #[allow(clippy::too_many_arguments)]
204 pub fn new(
205 c0: f64,
206 c1: f64,
207 c2: f64,
208 c3: f64,
209 c4: f64,
210 _c5: f64,
211 c6: f64,
212 c7: f64,
213 c8: f64,
214 c9: f64,
215 ) -> Self {
216 Self {
217 ndelay: c0,
218 dv_start: c2,
219 ttrans: c4,
220 c0m: c6,
221 c1m: c7,
222 vmin: c8,
223 vmax: c9,
224 slope: (c3 - c2) / c1,
225 }
226 }
227
228 /// Build a mobility calibrator straight from a `TimsCalibration` row.
229 pub fn from_calibration(cal: &crate::data::meta::TimsCalibration) -> Self {
230 Self::new(
231 cal.c0, cal.c1, cal.c2, cal.c3, cal.c4, cal.c5, cal.c6, cal.c7, cal.c8, cal.c9,
232 )
233 }
234
235 /// scan index -> trapping voltage (clamped to the valid window).
236 #[inline]
237 fn voltage(&self, scan: f64) -> f64 {
238 let v = self.dv_start + self.slope * (scan - self.ttrans - self.ndelay);
239 v.clamp(self.vmin, self.vmax)
240 }
241
242 /// scan index -> 1/K0 (inverse reduced ion mobility).
243 pub fn scan_to_one_over_k0(&self, scan: u32) -> f64 {
244 1.0 / (self.c0m + self.c1m / self.voltage(scan as f64))
245 }
246
247 /// 1/K0 -> nearest scan index (exact algebraic inverse, then round).
248 pub fn one_over_k0_to_scan(&self, one_over_k0: f64) -> u32 {
249 // invert 1/K0 = 1/(C0m + C1m/V) -> V, then V -> scan
250 let inv = 1.0 / one_over_k0;
251 let v = self.c1m / (inv - self.c0m);
252 let scan = (v - self.dv_start) / self.slope + self.ttrans + self.ndelay;
253 scan.round().max(0.0) as u32
254 }
255}
256
257#[cfg(test)]
258mod tests {
259 use super::*;
260
261 /// `synchro-hela.d`: `MzCalibration` ModelType 1, frame 1 temperatures.
262 fn model1() -> MzCalibrator {
263 MzCalibrator::new(
264 1,
265 0.2,
266 25585.2,
267 25.432518231649194,
268 22.91046717419208,
269 21.0,
270 0.0,
271 319.4836862507276,
272 156650.78463959479,
273 -4.7797946742077594e-05,
274 0.0,
275 2.6639979911791643e-05,
276 25.454778878245186,
277 23.065621750666164,
278 )
279 }
280
281 /// `G8602.d`: `MzCalibration` ModelType 2, frame 1 temperatures. Note that
282 /// C3 duplicates C0 and C4 duplicates C2 in this model — the calibrator has
283 /// to ignore both.
284 fn model2() -> MzCalibrator {
285 MzCalibrator::new(
286 2,
287 0.2,
288 18290.2,
289 25.3488367593942,
290 25.618275892137202,
291 20.0,
292 0.0,
293 319.95445850882476,
294 154831.91077331622,
295 -0.0005550791433026118,
296 319.95445850882476,
297 -0.0005550791433026118,
298 25.37849141449188,
299 26.171382195221447,
300 )
301 }
302
303 fn ppm(got: f64, want: f64) -> f64 {
304 (got - want).abs() / want * 1e6
305 }
306
307 /// ModelType 1 is fully modelled, so it must reproduce
308 /// `tims_index_to_mz` bit-for-bit (reference values read from the SDK).
309 #[test]
310 fn model1_is_bit_exact_against_the_sdk() {
311 let cal = model1();
312 for (tof, sdk) in [
313 (1u32, 100.00058181811438),
314 (50000, 194.8219835433269),
315 (150000, 478.458543155091),
316 (250000, 887.4159830461258),
317 (350000, 1421.6944158188035),
318 ] {
319 let got = cal.tof_to_mz(tof);
320 assert!(
321 (got - sdk).abs() / sdk < 1e-12,
322 "tof {tof}: got {got}, SDK {sdk}"
323 );
324 }
325 }
326
327 /// In ModelType 2 the C3/C4 columns are duplicates of C0/C2, not the cubic
328 /// term and the reduced-mass shift, so the calibrator must zero both.
329 #[test]
330 fn model2_drops_duplicated_c3_and_c4() {
331 let cal = model2();
332 assert_eq!(cal.c3, 0.0);
333 assert_eq!(cal.c4, 0.0);
334 }
335
336 /// The ModelType-2 fine correction is windowed to `[C5, C6]` in m/z, so
337 /// *below* that window (here C5 = 225.95) the C0/C1/C2 base curve is the
338 /// whole model and must match the SDK exactly. This is what regresses --
339 /// by 4.6 to 11.1 ppm -- if the duplicated C4 is subtracted as if it were m0.
340 #[test]
341 fn model2_is_exact_below_the_correction_window() {
342 let cal = model2();
343 for (tof, sdk) in [(1u32, 50.00106408611359), (50000, 121.13087702783326)] {
344 let got = cal.tof_to_mz(tof);
345 assert!(ppm(got, sdk) < 0.01, "tof {tof}: got {got}, SDK {sdk}");
346 }
347 }
348
349 /// Inside the correction window the unmodelled C8..C14 polynomial leaves a
350 /// small residual; it stays within a few ppm of the SDK.
351 #[test]
352 fn model2_is_within_a_few_ppm_inside_the_correction_window() {
353 let cal = model2();
354 for (tof, sdk) in [
355 (150000u32, 356.29356381674006),
356 (250000, 715.3229449956262),
357 (350000, 1198.2257163072832),
358 ] {
359 let got = cal.tof_to_mz(tof);
360 assert!(ppm(got, sdk) < 3.0, "tof {tof}: got {got}, SDK {sdk}");
361 }
362 }
363
364 #[test]
365 fn mz_tof_round_trips() {
366 for cal in [model1(), model2()] {
367 for tof in [1u32, 50_000, 150_000, 250_000, 350_000] {
368 let back = cal.mz_to_tof(cal.tof_to_mz(tof));
369 assert!(
370 back.abs_diff(tof) <= 1,
371 "round trip {tof} -> {} -> {back}",
372 cal.tof_to_mz(tof)
373 );
374 }
375 }
376 }
377
378 /// `synchro-hela.d` `TimsCalibration` (ModelType 2) against
379 /// `tims_scannum_to_oneoverk0`: the mobility model is the SDK computation,
380 /// so agreement is at the floating-point rounding limit.
381 #[test]
382 fn mobility_matches_the_sdk() {
383 let cal = MobilityCalibrator::new(
384 1.0,
385 926.0,
386 174.99089000590644,
387 89.34134412781896,
388 33.333333333333336,
389 1.0,
390 0.031163511286580903,
391 129.15504635187241,
392 12.75853947905552,
393 3135.217073581985,
394 );
395 for (scan, sdk) in [
396 (0u32, 1.3226192435624091),
397 (1, 1.321960900558156),
398 (100, 1.2566451818439914),
399 (400, 1.0570143319893166),
400 (900, 0.7184888610899344),
401 ] {
402 let got = cal.scan_to_one_over_k0(scan);
403 assert!(
404 (got - sdk).abs() / sdk < 1e-12,
405 "scan {scan}: got {got}, SDK {sdk}"
406 );
407 assert_eq!(cal.one_over_k0_to_scan(sdk), scan);
408 }
409 }
410}