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/// Temperature compensation (`T1/T2` = reference temps in `MzCalibration`,
37/// `dC1/dC2` its sensitivities, `T1f/T2f` = per-frame `Frames.T1/Frames.T2`):
38/// `tc = 1 + (dC1*(T1-T1f) + dC2*(T2-T2f)) / 1e6`, applied as `C1 *= tc`.
39#[derive(Debug, Clone)]
40pub struct MzCalibrator {
41 pub timebase: f64,
42 pub delay: f64,
43 pub c0: f64,
44 pub b: f64, // sqrt(1e12 / (C1 * tc))
45 pub c2: f64,
46 pub c3: f64,
47 pub c4: f64,
48}
49
50impl MzCalibrator {
51 /// Build a calibrator from raw `MzCalibration` columns + per-frame temps.
52 ///
53 /// `model_type` selects whether the C2/C3 curve terms are honoured (type 1)
54 /// or zeroed (type 2, base curve only).
55 #[allow(clippy::too_many_arguments)]
56 pub fn new(
57 model_type: i64,
58 timebase: f64,
59 delay: f64,
60 t1_ref: f64,
61 t2_ref: f64,
62 dc1: f64,
63 dc2: f64,
64 c0: f64,
65 c1: f64,
66 c2: f64,
67 c3: f64,
68 c4: f64,
69 t1_frame: f64,
70 t2_frame: f64,
71 ) -> Self {
72 let tc = 1.0 + (dc1 * (t1_ref - t1_frame) + dc2 * (t2_ref - t2_frame)) / 1.0e6;
73 let b = (1.0e12 / (c1 * tc)).sqrt();
74 // Both models share the quadratic-in-sqrt(m) curve `t = C0 + b*s + C2*s^2`
75 // (empirically: a0->C0, a1->sqrt(1e12/C1), a2->C2). The cubic `C3*s^3`
76 // term is real only for ModelType 1; in ModelType 2 the C3 column is a
77 // duplicate of C0 and must be dropped. ModelType 2 additionally carries a
78 // C8..C14 fine correction (~few ppm, worst at low m/z) that is NOT an
79 // additive polynomial in m/z and is left unmodelled here.
80 let c2 = c2 / tc;
81 let c3 = if model_type == 1 { c3 } else { 0.0 };
82 Self { timebase, delay, c0, b, c2, c3, c4 }
83 }
84
85 /// Flight time (digitizer units) for a TOF index.
86 #[inline]
87 fn tof_index_to_time(&self, tof_index: f64) -> f64 {
88 tof_index * self.timebase + self.delay
89 }
90
91 /// TOF index -> m/z. Inverts `t = C0 + b*s + c2*s^2 + c3*s^3` for `s`.
92 pub fn tof_to_mz(&self, tof_index: u32) -> f64 {
93 let t = self.tof_index_to_time(tof_index as f64);
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;
99 for _ in 0..8 {
100 let f = self.c0 + self.b * s + self.c2 * s * s + self.c3 * s * s * s - t;
101 let df = self.b + 2.0 * self.c2 * s + 3.0 * self.c3 * s * s;
102 if df == 0.0 {
103 break;
104 }
105 let step = f / df;
106 s -= step;
107 if step.abs() < 1e-12 {
108 break;
109 }
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 };
128 s * s - self.c4
129 }
130
131 /// m/z -> TOF index (forward direction, always closed form).
132 pub fn mz_to_tof(&self, mz: f64) -> u32 {
133 let s = (mz + self.c4).max(0.0).sqrt();
134 let t = self.c0 + self.b * s + self.c2 * s * s + self.c3 * s * s * s;
135 (((t - self.delay) / self.timebase).round()).max(0.0) as u32
136 }
137}
138
139/// Ion-mobility axis calibration (Bruker "model type 2", the only TIMS model).
140///
141/// Two steps, both exact ports of PAPPSO `timsframebase.cpp`:
142/// 1. scan -> trapping voltage: `V = dv_start + slope*(scan - ttrans - ndelay)`
143/// with `slope = (dv_end - dv_start) / ncycles`. V must lie in [vmin,vmax].
144/// 2. voltage -> inverse mobility: `1/K0 = 1 / (C0m + C1m / V)`.
145///
146/// Coefficient meaning (columns of the `TimsCalibration` table, ModelType 2):
147/// * `C0` = `ndelay` — scan offset (delay), subtracted before scaling.
148/// * `C1` = `ncycles` — number of TIMS cycles; sets the voltage-vs-scan slope.
149/// * `C2` = `dv_start` — trapping voltage at the start of the ramp.
150/// * `C3` = `dv_end` — trapping voltage at the end of the ramp.
151/// * `C4` = `ttrans` — transit time in cycles, subtracted before scaling.
152/// * `C5` — unused by the mobility formula (polynomial grade flag).
153/// * `C6` = `C0m` — additive constant of the mobility reciprocal.
154/// * `C7` = `C1m` — voltage-scaled term of the mobility reciprocal.
155/// * `C8` = `vmin` — lower voltage validity bound.
156/// * `C9` = `vmax` — upper voltage validity bound.
157#[derive(Debug, Clone)]
158pub struct MobilityCalibrator {
159 pub ndelay: f64,
160 pub dv_start: f64,
161 pub ttrans: f64,
162 pub c0m: f64,
163 pub c1m: f64,
164 pub vmin: f64,
165 pub vmax: f64,
166 pub slope: f64,
167}
168
169impl MobilityCalibrator {
170 /// Build from raw `TimsCalibration` C0..C9 (ModelType must be 2).
171 #[allow(clippy::too_many_arguments)]
172 pub fn new(
173 c0: f64,
174 c1: f64,
175 c2: f64,
176 c3: f64,
177 c4: f64,
178 _c5: f64,
179 c6: f64,
180 c7: f64,
181 c8: f64,
182 c9: f64,
183 ) -> Self {
184 Self {
185 ndelay: c0,
186 dv_start: c2,
187 ttrans: c4,
188 c0m: c6,
189 c1m: c7,
190 vmin: c8,
191 vmax: c9,
192 slope: (c3 - c2) / c1,
193 }
194 }
195
196 /// scan index -> trapping voltage (clamped to the valid window).
197 #[inline]
198 fn voltage(&self, scan: f64) -> f64 {
199 let v = self.dv_start + self.slope * (scan - self.ttrans - self.ndelay);
200 v.clamp(self.vmin, self.vmax)
201 }
202
203 /// scan index -> 1/K0 (inverse reduced ion mobility).
204 pub fn scan_to_one_over_k0(&self, scan: u32) -> f64 {
205 1.0 / (self.c0m + self.c1m / self.voltage(scan as f64))
206 }
207
208 /// 1/K0 -> nearest scan index (exact algebraic inverse, then round).
209 pub fn one_over_k0_to_scan(&self, one_over_k0: f64) -> u32 {
210 // invert 1/K0 = 1/(C0m + C1m/V) -> V, then V -> scan
211 let inv = 1.0 / one_over_k0;
212 let v = self.c1m / (inv - self.c0m);
213 let scan = (v - self.dv_start) / self.slope + self.ttrans + self.ndelay;
214 scan.round().max(0.0) as u32
215 }
216}