Skip to main content

rustdf/sim/
dda.rs

1use mscore::algorithm::isotope::calculate_precursor_transmission_factor;
2use crate::sim::containers::IsotopeTransmissionMode;
3use mscore::data::peptide::{PeptideIon, PeptideProductIonSeriesCollection};
4use mscore::data::spectrum::{IndexedMzSpectrum, MsType, MzSpectrum};
5use mscore::simulation::annotation::{
6    MzSpectrumAnnotated, PeakAnnotation, TimsFrameAnnotated, TimsSpectrumAnnotated,
7};
8use mscore::timstof::frame::TimsFrame;
9use mscore::timstof::quadrupole::{IonTransmission, TimsTransmissionDDA};
10use mscore::timstof::spectrum::TimsSpectrum;
11use std::collections::{BTreeMap, HashSet};
12use std::path::Path;
13use std::sync::Arc;
14
15use rand::Rng;
16use rayon::prelude::*;
17use crate::sim::containers::IsotopeTransmissionConfig;
18use crate::sim::scheme::InstrumentCapabilities;
19use crate::sim::handle::{FragmentIonsWithComplementary, TimsTofSyntheticsDataHandle};
20use crate::sim::precursor::TimsTofSyntheticsPrecursorFrameBuilder;
21
22pub struct TimsTofSyntheticsFrameBuilderDDA {
23    pub path: String,
24    pub precursor_frame_builder: TimsTofSyntheticsPrecursorFrameBuilder,
25    pub transmission_settings: TimsTransmissionDDA,
26    pub fragment_ions:
27        Option<BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrum>)>>,
28    pub fragment_ions_annotated: Option<
29        BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrumAnnotated>)>,
30    >,
31    /// Configuration for quad-selection dependent isotope transmission (already
32    /// gated by the instrument capabilities at construction — P5e).
33    pub isotope_transmission_config: IsotopeTransmissionConfig,
34    /// Fragment ions with complementary distribution data (only populated when isotope_transmission_config.enabled)
35    pub fragment_ions_with_complementary: Option<BTreeMap<(u32, i8, i32), FragmentIonsWithComplementary>>,
36    /// Physical instrument capabilities (P5e). Default = Bruker timsTOF.
37    pub capabilities: InstrumentCapabilities,
38}
39
40impl TimsTofSyntheticsFrameBuilderDDA {
41    /// Create a new DDA frame builder.
42    ///
43    /// # Arguments
44    ///
45    /// * `path` - Path to the simulation database
46    /// * `with_annotations` - Whether to include annotations
47    /// * `num_threads` - Number of threads for parallel processing
48    /// * `isotope_config` - Optional configuration for quad-dependent isotope transmission
49    pub fn new(
50        path: &Path,
51        with_annotations: bool,
52        num_threads: usize,
53        isotope_config: Option<IsotopeTransmissionConfig>,
54    ) -> Self {
55        Self::new_with_capabilities(
56            path,
57            with_annotations,
58            num_threads,
59            isotope_config,
60            InstrumentCapabilities::default(),
61        )
62    }
63
64    /// Like [`Self::new`], but with explicit instrument capabilities (P5e). The
65    /// isotope-transmission config is gated by them (an instrument without
66    /// mobility-dependent quad isotope transmission forces the mode to None).
67    /// Default capabilities = Bruker timsTOF, so [`Self::new`] is unchanged.
68    pub fn new_with_capabilities(
69        path: &Path,
70        with_annotations: bool,
71        num_threads: usize,
72        isotope_config: Option<IsotopeTransmissionConfig>,
73        capabilities: InstrumentCapabilities,
74    ) -> Self {
75        let handle = TimsTofSyntheticsDataHandle::new(path).unwrap();
76        // P5b: refuse to render fragments stored under an incompatible prediction
77        // set (e.g. a future Thermo set whose CE encoding the render keying can't
78        // resolve). Bruker/legacy sets pass; this never fires for current DBs.
79        handle
80            .read_prediction_set()
81            .expect("read prediction set")
82            .assert_render_compatible()
83            .expect("incompatible fragment prediction set for this renderer");
84        let fragment_ions_raw = handle.read_fragment_ions().unwrap();
85        let transmission_settings = handle.get_transmission_dda();
86
87        let synthetics = TimsTofSyntheticsPrecursorFrameBuilder::new(path).unwrap();
88        let config = isotope_config.unwrap_or_default().gated_by(capabilities);
89
90        // Build fragment ions with complementary data if any transmission mode is enabled
91        let fragment_ions_with_complementary = if config.is_enabled() {
92            Some(TimsTofSyntheticsDataHandle::build_fragment_ions_with_transmission_data(
93                &synthetics.peptides,
94                &fragment_ions_raw,
95                num_threads,
96            ))
97        } else {
98            None
99        };
100
101        match with_annotations {
102            true => {
103                let fragment_ions =
104                    Some(TimsTofSyntheticsDataHandle::build_fragment_ions_annotated(
105                        &synthetics.peptides,
106                        &fragment_ions_raw,
107                        num_threads,
108                    ));
109                Self {
110                    path: path.to_str().unwrap().to_string(),
111                    precursor_frame_builder: synthetics,
112                    transmission_settings,
113                    fragment_ions: None,
114                    fragment_ions_annotated: fragment_ions,
115                    isotope_transmission_config: config,
116                    fragment_ions_with_complementary,
117                    capabilities,
118                }
119            }
120            false => {
121                let fragment_ions = Some(TimsTofSyntheticsDataHandle::build_fragment_ions(
122                    &synthetics.peptides,
123                    &fragment_ions_raw,
124                    num_threads,
125                ));
126                Self {
127                    path: path.to_str().unwrap().to_string(),
128                    precursor_frame_builder: synthetics,
129                    transmission_settings,
130                    fragment_ions,
131                    fragment_ions_annotated: None,
132                    isotope_transmission_config: config,
133                    fragment_ions_with_complementary,
134                    capabilities,
135                }
136            }
137        }
138    }
139    /// Construct an eager DDA frame builder from already-loaded, in-memory
140    /// entities instead of reading the whole database.
141    ///
142    /// Used by the lazy DDA builder so it can delegate per-batch frame
143    /// construction to the *same* (correct) eager algorithm instead of keeping a
144    /// second, divergent copy. Only the non-annotated path is built (lazy loading
145    /// does not support annotations). `fragment_ions_raw` is the per-batch slice
146    /// of fragment ions for the supplied peptides, so memory stays bounded to the
147    /// batch.
148    pub fn from_entities(
149        precursor_frame_builder: TimsTofSyntheticsPrecursorFrameBuilder,
150        transmission_settings: TimsTransmissionDDA,
151        fragment_ions_raw: Vec<crate::sim::containers::FragmentIonSim>,
152        isotope_config: Option<IsotopeTransmissionConfig>,
153        fragmentation: bool,
154        num_threads: usize,
155    ) -> Self {
156        // Lazy DDA (delegates here) is Bruker timsTOF; gate is a no-op. P6 will
157        // thread real capabilities through the lazy path.
158        let config = isotope_config
159            .unwrap_or_default()
160            .gated_by(InstrumentCapabilities::default());
161
162        // The fragment isotope map is only consumed when fragmentation is on
163        // (build_ms2_frame's `true` branch). Expanding all predicted fragment
164        // intensities into isotope spectra is expensive, so skip it entirely in
165        // no-fragmentation mode (the `false` branch never touches these maps).
166        let (fragment_ions, fragment_ions_with_complementary) = if fragmentation {
167            let with_complementary = if config.is_enabled() {
168                Some(TimsTofSyntheticsDataHandle::build_fragment_ions_with_transmission_data(
169                    &precursor_frame_builder.peptides,
170                    &fragment_ions_raw,
171                    num_threads,
172                ))
173            } else {
174                None
175            };
176            let fragment_ions = Some(TimsTofSyntheticsDataHandle::build_fragment_ions(
177                &precursor_frame_builder.peptides,
178                &fragment_ions_raw,
179                num_threads,
180            ));
181            (fragment_ions, with_complementary)
182        } else {
183            (None, None)
184        };
185
186        Self {
187            path: String::new(),
188            precursor_frame_builder,
189            transmission_settings,
190            fragment_ions,
191            fragment_ions_annotated: None,
192            isotope_transmission_config: config,
193            fragment_ions_with_complementary,
194            capabilities: InstrumentCapabilities::default(),
195        }
196    }
197
198    /// Build a frame for DDA synthetic experiment
199    ///
200    /// # Arguments
201    ///
202    /// * `frame_id` - The frame id
203    /// * `fragmentation` - A boolean indicating if fragmentation is enabled, if false, the frame has same mz distribution as the precursor frame but will be quadrupole filtered
204    ///
205    /// # Returns
206    ///
207    /// A TimsFrame
208    ///
209    pub fn build_frame(
210        &self,
211        frame_id: u32,
212        fragmentation: bool,
213        mz_noise_precursor: bool,
214        uniform: bool,
215        precursor_noise_ppm: f64,
216        mz_noise_fragment: bool,
217        fragment_noise_ppm: f64,
218        right_drag: bool,
219    ) -> TimsFrame {
220        // determine if the frame is a precursor frame
221        match self
222            .precursor_frame_builder
223            .precursor_frame_id_set
224            .contains(&frame_id)
225        {
226            true => self.build_ms1_frame(
227                frame_id,
228                mz_noise_precursor,
229                uniform,
230                precursor_noise_ppm,
231                right_drag,
232            ),
233            false => self.build_ms2_frame(
234                frame_id,
235                fragmentation,
236                mz_noise_fragment,
237                uniform,
238                fragment_noise_ppm,
239                right_drag,
240            ),
241        }
242    }
243
244    pub fn build_frame_annotated(
245        &self,
246        frame_id: u32,
247        fragmentation: bool,
248        mz_noise_precursor: bool,
249        uniform: bool,
250        precursor_noise_ppm: f64,
251        mz_noise_fragment: bool,
252        fragment_noise_ppm: f64,
253        right_drag: bool,
254    ) -> TimsFrameAnnotated {
255        match self
256            .precursor_frame_builder
257            .precursor_frame_id_set
258            .contains(&frame_id)
259        {
260            true => self.build_ms1_frame_annotated(
261                frame_id,
262                mz_noise_precursor,
263                uniform,
264                precursor_noise_ppm,
265                right_drag,
266            ),
267            false => self.build_ms2_frame_annotated(
268                frame_id,
269                fragmentation,
270                mz_noise_fragment,
271                uniform,
272                fragment_noise_ppm,
273                right_drag,
274            ),
275        }
276    }
277
278    pub fn get_fragment_ion_ids(&self, precursor_frame_ids: Vec<u32>) -> Vec<u32> {
279        let mut peptide_ids: HashSet<u32> = HashSet::new();
280        // get all peptide ids for the precursor frame ids
281        for frame_id in precursor_frame_ids {
282            for (peptide_id, peptide) in self.precursor_frame_builder.peptides.iter() {
283                if peptide.frame_start <= frame_id && peptide.frame_end >= frame_id {
284                    peptide_ids.insert(*peptide_id);
285                }
286            }
287        }
288        // get all ion ids for the peptide ids
289        let mut result: Vec<u32> = Vec::new();
290        for item in peptide_ids {
291            let ions = self.precursor_frame_builder.ions.get(&item).unwrap();
292            for ion in ions.iter() {
293                result.push(ion.ion_id);
294            }
295        }
296        result
297    }
298
299    pub fn build_frames(
300        &self,
301        frame_ids: Vec<u32>,
302        fragmentation: bool,
303        mz_noise_precursor: bool,
304        uniform: bool,
305        precursor_noise_ppm: f64,
306        mz_noise_fragment: bool,
307        fragment_noise_ppm: f64,
308        right_drag: bool,
309        num_threads: usize,
310    ) -> Vec<TimsFrame> {
311        // Use thread pool with custom parallelism
312        let pool = rayon::ThreadPoolBuilder::new()
313            .num_threads(num_threads)
314            .build()
315            .unwrap();
316
317        pool.install(|| {
318            // Use indexed parallel iteration to maintain order, avoiding post-sort
319            let mut tims_frames: Vec<TimsFrame> = Vec::with_capacity(frame_ids.len());
320            unsafe { tims_frames.set_len(frame_ids.len()); }
321
322            frame_ids.par_iter().enumerate().for_each(|(idx, frame_id)| {
323                let frame = self.build_frame(
324                    *frame_id,
325                    fragmentation,
326                    mz_noise_precursor,
327                    uniform,
328                    precursor_noise_ppm,
329                    mz_noise_fragment,
330                    fragment_noise_ppm,
331                    right_drag,
332                );
333                unsafe {
334                    let ptr = tims_frames.as_ptr() as *mut TimsFrame;
335                    std::ptr::write(ptr.add(idx), frame);
336                }
337            });
338
339            tims_frames
340        })
341    }
342    pub fn build_frames_annotated(
343        &self,
344        frame_ids: Vec<u32>,
345        fragmentation: bool,
346        mz_noise_precursor: bool,
347        uniform: bool,
348        precursor_noise_ppm: f64,
349        mz_noise_fragment: bool,
350        fragment_noise_ppm: f64,
351        right_drag: bool,
352        num_threads: usize,
353    ) -> Vec<TimsFrameAnnotated> {
354        // Use thread pool with custom parallelism
355        let pool = rayon::ThreadPoolBuilder::new()
356            .num_threads(num_threads)
357            .build()
358            .unwrap();
359
360        pool.install(|| {
361            // Use indexed parallel iteration to maintain order, avoiding post-sort
362            let mut tims_frames: Vec<TimsFrameAnnotated> = Vec::with_capacity(frame_ids.len());
363            unsafe { tims_frames.set_len(frame_ids.len()); }
364
365            frame_ids.par_iter().enumerate().for_each(|(idx, frame_id)| {
366                let frame = self.build_frame_annotated(
367                    *frame_id,
368                    fragmentation,
369                    mz_noise_precursor,
370                    uniform,
371                    precursor_noise_ppm,
372                    mz_noise_fragment,
373                    fragment_noise_ppm,
374                    right_drag,
375                );
376                unsafe {
377                    let ptr = tims_frames.as_ptr() as *mut TimsFrameAnnotated;
378                    std::ptr::write(ptr.add(idx), frame);
379                }
380            });
381
382            tims_frames
383        })
384    }
385
386    fn build_ms1_frame(
387        &self,
388        frame_id: u32,
389        mz_noise_precursor: bool,
390        uniform: bool,
391        precursor_ppm: f64,
392        right_drag: bool,
393    ) -> TimsFrame {
394        let mut tims_frame = self.precursor_frame_builder.build_precursor_frame(
395            frame_id,
396            mz_noise_precursor,
397            uniform,
398            precursor_ppm,
399            right_drag,
400        );
401        let intensities_rounded = tims_frame
402            .ims_frame
403            .intensity
404            .iter()
405            .map(|x| x.round())
406            .collect::<Vec<_>>();
407        tims_frame.ims_frame.intensity = Arc::new(intensities_rounded);
408        tims_frame
409    }
410
411    fn build_ms1_frame_annotated(
412        &self,
413        frame_id: u32,
414        mz_noise_precursor: bool,
415        uniform: bool,
416        precursor_ppm: f64,
417        right_drag: bool,
418    ) -> TimsFrameAnnotated {
419        let mut tims_frame = self
420            .precursor_frame_builder
421            .build_precursor_frame_annotated(
422                frame_id,
423                mz_noise_precursor,
424                uniform,
425                precursor_ppm,
426                right_drag,
427            );
428        let intensities_rounded = tims_frame
429            .intensity
430            .iter()
431            .map(|x| x.round())
432            .collect::<Vec<_>>();
433        tims_frame.intensity = intensities_rounded;
434        tims_frame
435    }
436
437    fn build_ms2_frame(
438        &self,
439        frame_id: u32,
440        fragmentation: bool,
441        mz_noise_fragment: bool,
442        uniform: bool,
443        fragment_ppm: f64,
444        right_drag: bool,
445    ) -> TimsFrame {
446        match fragmentation {
447            false => {
448                let mut frame = self.transmission_settings.transmit_tims_frame(
449                    &self.build_ms1_frame(
450                        frame_id,
451                        mz_noise_fragment,
452                        uniform,
453                        fragment_ppm,
454                        right_drag,
455                    ),
456                    None,
457                );
458                let intensities_rounded = frame
459                    .ims_frame
460                    .intensity
461                    .iter()
462                    .map(|x| x.round())
463                    .collect::<Vec<_>>();
464                frame.ims_frame.intensity = Arc::new(intensities_rounded);
465                // DDA MS2 frames are PASEF fragment frames (MsMsType 8). The
466                // no-fragmentation mode still produces a DDA fragment frame
467                // (quad-filtered precursor), so it must keep the DDA type, not
468                // FragmentDia (9). The fragmentation=true branch already uses
469                // FragmentDda; this matches it.
470                frame.ms_type = MsType::FragmentDda;
471                frame
472            }
473            true => {
474                let mut frame = self.build_fragment_frame(
475                    frame_id,
476                    &self.fragment_ions.as_ref().unwrap(),
477                    mz_noise_fragment,
478                    uniform,
479                    fragment_ppm,
480                    None,
481                    None,
482                    None,
483                    Some(right_drag),
484                );
485                let intensities_rounded = frame
486                    .ims_frame
487                    .intensity
488                    .iter()
489                    .map(|x| x.round())
490                    .collect::<Vec<_>>();
491                frame.ims_frame.intensity = Arc::new(intensities_rounded);
492                frame
493            }
494        }
495    }
496
497    fn build_ms2_frame_annotated(
498        &self,
499        frame_id: u32,
500        fragmentation: bool,
501        mz_noise_fragment: bool,
502        uniform: bool,
503        fragment_ppm: f64,
504        right_drag: bool,
505    ) -> TimsFrameAnnotated {
506        match fragmentation {
507            false => {
508                let mut frame = self.transmission_settings.transmit_tims_frame_annotated(
509                    &self.build_ms1_frame_annotated(
510                        frame_id,
511                        mz_noise_fragment,
512                        uniform,
513                        fragment_ppm,
514                        right_drag,
515                    ),
516                    None,
517                );
518                let intensities_rounded = frame
519                    .intensity
520                    .iter()
521                    .map(|x| x.round())
522                    .collect::<Vec<_>>();
523                frame.intensity = intensities_rounded;
524                // See build_ms2_frame: DDA no-frag MS2 stays FragmentDda (8).
525                frame.ms_type = MsType::FragmentDda;
526                frame
527            }
528            true => {
529                let mut frame = self.build_fragment_frame_annotated(
530                    frame_id,
531                    &self.fragment_ions_annotated.as_ref().unwrap(),
532                    mz_noise_fragment,
533                    uniform,
534                    fragment_ppm,
535                    None,
536                    None,
537                    None,
538                    Some(right_drag),
539                );
540                let intensities_rounded = frame
541                    .intensity
542                    .iter()
543                    .map(|x| x.round())
544                    .collect::<Vec<_>>();
545                frame.intensity = intensities_rounded;
546                frame
547            }
548        }
549    }
550
551    /// Build a fragment frame
552    ///
553    /// # Arguments
554    ///
555    /// * `frame_id` - The frame id
556    /// * `mz_min` - The minimum m/z value in fragment spectrum
557    /// * `mz_max` - The maximum m/z value in fragment spectrum
558    /// * `intensity_min` - The minimum intensity value in fragment spectrum
559    ///
560    /// # Returns
561    ///
562    /// A TimsFrame
563    ///
564    fn build_fragment_frame(
565        &self,
566        frame_id: u32,
567        fragment_ions: &BTreeMap<
568            (u32, i8, i32),
569            (PeptideProductIonSeriesCollection, Vec<MzSpectrum>),
570        >,
571        mz_noise_fragment: bool,
572        uniform: bool,
573        fragment_ppm: f64,
574        mz_min: Option<f64>,
575        mz_max: Option<f64>,
576        intensity_min: Option<f64>,
577        right_drag: Option<bool>,
578    ) -> TimsFrame {
579        // Cache frame-level lookups
580        let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
581            MsType::Unknown
582        } else {
583            MsType::FragmentDda
584        };
585
586        let rt = *self.precursor_frame_builder.frame_to_rt.get(&frame_id)
587            .expect("frame_to_rt should always have this frame") as f64;
588        let right_drag_val = right_drag.unwrap_or(false);
589        let mz_min_val = mz_min.unwrap_or(100.0);
590        let mz_max_val = mz_max.unwrap_or(1700.0);
591        let intensity_min_val = intensity_min.unwrap_or(1.0);
592
593        // Get PASEF meta for this frame - these define which precursors were SELECTED
594        let Some(pasef_meta) = self.transmission_settings.pasef_meta.get(&(frame_id as i32)) else {
595            return TimsFrame::new(frame_id as i32, ms_type, rt, vec![], vec![], vec![], vec![], vec![]);
596        };
597
598        // Preallocate with estimated capacity
599        let estimated_capacity = pasef_meta.len() * 4;
600        let mut tims_spectra: Vec<TimsSpectrum> = Vec::with_capacity(estimated_capacity);
601
602        // Iterate over PASEF meta entries - each represents a SELECTED precursor
603        for meta in pasef_meta.iter() {
604            // Get the selected ion_id from the PASEF precursor field
605            let ion_id = meta.precursor as u32;
606
607            // Look up peptide_id and charge from the ion_id
608            let Some(&(peptide_id, charge_state)) = self.precursor_frame_builder.ion_id_to_peptide_charge.get(&ion_id) else {
609                continue;
610            };
611
612            // Get peptide info
613            let Some(peptide) = self.precursor_frame_builder.peptides.get(&peptide_id) else {
614                continue;
615            };
616
617            // Check if this peptide is active in this frame
618            if frame_id < peptide.frame_start || frame_id > peptide.frame_end {
619                continue;
620            }
621
622            // Get ion info from peptide_to_ions
623            let Some((ion_abundances, scan_occurrences, scan_abundances, charges, spectra)) = self
624                .precursor_frame_builder
625                .peptide_to_ions
626                .get(&peptide_id)
627            else {
628                continue;
629            };
630
631            // Find the index of this specific charge state
632            let Some(ion_index) = charges.iter().position(|&c| c == charge_state) else {
633                continue;
634            };
635
636            let ion_abundance = ion_abundances[ion_index];
637            let all_scan_occurrence = &scan_occurrences[ion_index];
638            let all_scan_abundance = &scan_abundances[ion_index];
639            let spectrum = &spectra[ion_index];
640            let total_events = *self.precursor_frame_builder.peptide_to_events.get(&peptide_id).unwrap();
641
642            // Get frame abundance for this peptide
643            let frame_abundance = self.precursor_frame_builder.frame_to_abundances
644                .get(&frame_id)
645                .and_then(|(pep_ids, abundances)| {
646                    pep_ids.iter().position(|&p| p == peptide_id)
647                        .map(|idx| abundances[idx])
648                })
649                .unwrap_or(0.0);
650
651            if frame_abundance == 0.0 {
652                continue;
653            }
654
655            // Collision energy from the PASEF meta
656            let collision_energy = meta.collision_energy;
657            // Resolve the fragment CE key, tolerant to ~0.1 eV quantization noise
658            // (DDA pasef CE is full-precision; the stored key is f32-quantized).
659            let Some(collision_energy_quantized) = crate::sim::handle::resolve_fragment_ce_key(
660                fragment_ions, peptide_id, charge_state, collision_energy,
661            ) else {
662                // Fail loud if fragments exist for this precursor but none near the
663                // applied CE: the prediction set does not cover this instrument's CE
664                // (P5b). A precursor with NO predicted fragments at all is a
665                // legitimate skip. (For Bruker after the ±0.1 eV probe this never
666                // fires — verified 0 misses.)
667                if crate::sim::handle::fragment_prefix_exists(fragment_ions, peptide_id, charge_state) {
668                    panic!(
669                        "DDA fragment lookup miss: peptide {} charge {} applied CE {:.4} eV has \
670                         predicted fragments, but none within 0.1 eV — the prediction set does \
671                         not cover this instrument's collision energy",
672                        peptide_id, charge_state, collision_energy,
673                    );
674                }
675                continue;
676            };
677            let (_, fragment_series_vec) = fragment_ions
678                .get(&(peptide_id, charge_state, collision_energy_quantized))
679                .expect("resolve_fragment_ce_key returned a present key");
680
681            // Process scans within the PASEF selection window
682            for (scan, scan_abundance) in all_scan_occurrence.iter().zip(all_scan_abundance.iter()) {
683                // Check if scan is within the PASEF selection window
684                let scan_i32 = *scan as i32;
685                if scan_i32 < meta.scan_start || scan_i32 > meta.scan_end {
686                    continue;
687                }
688
689                // Get transmitted isotope indices based on config mode
690                let transmitted_indices = match self.isotope_transmission_config.mode {
691                    IsotopeTransmissionMode::None => HashSet::new(),
692                    IsotopeTransmissionMode::PrecursorScaling | IsotopeTransmissionMode::PerFragment => {
693                        self.transmission_settings.get_transmission_set(
694                            frame_id as i32,
695                            scan_i32,
696                            &spectrum.mz,
697                            Some(self.isotope_transmission_config.min_probability),
698                        )
699                    },
700                };
701
702                // Calculate abundance factor
703                let fraction_events = frame_abundance * scan_abundance * ion_abundance * total_events;
704
705                // Cache scan mobility
706                let scan_mobility = *self.precursor_frame_builder.scan_to_mobility.get(scan).unwrap() as f64;
707
708                // Calculate transmission factor for PrecursorScaling mode
709                let transmission_factor = if self.isotope_transmission_config.mode == IsotopeTransmissionMode::PrecursorScaling {
710                    if let Some(comp_data) = self.fragment_ions_with_complementary.as_ref() {
711                        if let Some(frag_data) = comp_data.get(&(peptide_id, charge_state, collision_energy_quantized)) {
712                            calculate_precursor_transmission_factor(
713                                &frag_data.precursor_isotope_distribution,
714                                &transmitted_indices,
715                            )
716                        } else {
717                            1.0
718                        }
719                    } else {
720                        1.0
721                    }
722                } else {
723                    1.0
724                };
725
726                // Complementary per-(peptide,charge,CE) data for PerFragment;
727                // looked up once (same for all series) and passed to the shared
728                // kernel — same lookup the transmission_factor block above uses.
729                let frag_data = self
730                    .fragment_ions_with_complementary
731                    .as_ref()
732                    .and_then(|c| c.get(&(peptide_id, charge_state, collision_energy_quantized)));
733
734                for (series_idx, fragment_ion_series) in fragment_series_vec.iter().enumerate() {
735                    let final_spectrum = crate::sim::dia::fragment_series_spectrum(
736                        self.isotope_transmission_config.mode,
737                        fragment_ion_series,
738                        series_idx,
739                        fraction_events,
740                        transmission_factor,
741                        frag_data,
742                        &transmitted_indices,
743                        self.isotope_transmission_config.max_isotopes,
744                    );
745
746                    let mz_spectrum = if mz_noise_fragment {
747                        if uniform {
748                            final_spectrum.add_mz_noise_uniform(fragment_ppm, right_drag_val)
749                        } else {
750                            final_spectrum.add_mz_noise_normal(fragment_ppm)
751                        }
752                    } else {
753                        final_spectrum
754                    };
755
756                    let spectrum_len = mz_spectrum.mz.len();
757                    tims_spectra.push(TimsSpectrum::new(
758                        frame_id as i32,
759                        *scan as i32,
760                        rt,
761                        scan_mobility,
762                        ms_type.clone(),
763                        IndexedMzSpectrum::from_mz_spectrum(
764                            vec![0; spectrum_len],
765                            mz_spectrum,
766                        ).filter_ranged(100.0, 1700.0, 1.0, 1e9),
767                    ));
768                }
769
770                // Add unfragmented precursor ions (survival) if configured
771                if self.isotope_transmission_config.has_precursor_survival() {
772                    let mut rng = rand::thread_rng();
773                    let survival_fraction = rng.gen_range(
774                        self.isotope_transmission_config.precursor_survival_min
775                        ..=self.isotope_transmission_config.precursor_survival_max
776                    );
777
778                    if survival_fraction > 0.0 {
779                        // Transmit the precursor spectrum through the quadrupole
780                        let precursor_transmitted = self.transmission_settings.transmit_spectrum(
781                            frame_id as i32,
782                            *scan as i32,
783                            spectrum.clone(),
784                            Some(self.isotope_transmission_config.min_probability),
785                        );
786
787                        if !precursor_transmitted.mz.is_empty() {
788                            // Scale by survival fraction and event count
789                            let precursor_scaled = precursor_transmitted * (fraction_events as f64 * survival_fraction);
790
791                            let precursor_mz_spectrum = if mz_noise_fragment {
792                                if uniform {
793                                    precursor_scaled.add_mz_noise_uniform(fragment_ppm, right_drag_val)
794                                } else {
795                                    precursor_scaled.add_mz_noise_normal(fragment_ppm)
796                                }
797                            } else {
798                                precursor_scaled
799                            };
800
801                            let precursor_len = precursor_mz_spectrum.mz.len();
802                            tims_spectra.push(TimsSpectrum::new(
803                                frame_id as i32,
804                                *scan as i32,
805                                rt,
806                                scan_mobility,
807                                ms_type.clone(),
808                                IndexedMzSpectrum::from_mz_spectrum(
809                                    vec![0; precursor_len],
810                                    precursor_mz_spectrum,
811                                ).filter_ranged(100.0, 1700.0, 1.0, 1e9),
812                            ));
813                        }
814                    }
815                }
816            }
817        }
818
819        if tims_spectra.is_empty() {
820            return TimsFrame::new(frame_id as i32, ms_type, rt, vec![], vec![], vec![], vec![], vec![]);
821        }
822
823        let tims_frame = TimsFrame::from_tims_spectra(tims_spectra);
824        tims_frame.filter_ranged(
825            mz_min_val,
826            mz_max_val,
827            0,
828            1000,
829            0.0,
830            10.0,
831            intensity_min_val,
832            1e9,
833            0,
834            i32::MAX,
835        )
836    }
837
838    pub fn build_fragment_frame_annotated(
839        &self,
840        frame_id: u32,
841        fragment_ions: &BTreeMap<
842            (u32, i8, i32),
843            (PeptideProductIonSeriesCollection, Vec<MzSpectrumAnnotated>),
844        >,
845        mz_noise_fragment: bool,
846        uniform: bool,
847        fragment_ppm: f64,
848        mz_min: Option<f64>,
849        mz_max: Option<f64>,
850        intensity_min: Option<f64>,
851        right_drag: Option<bool>,
852    ) -> TimsFrameAnnotated {
853        // Cache frame-level lookups
854        let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
855            MsType::Unknown
856        } else {
857            MsType::FragmentDda
858        };
859
860        let rt = *self.precursor_frame_builder.frame_to_rt.get(&frame_id).unwrap() as f64;
861        let right_drag_val = right_drag.unwrap_or(false);
862        let mz_min_val = mz_min.unwrap_or(100.0);
863        let mz_max_val = mz_max.unwrap_or(1700.0);
864        let intensity_min_val = intensity_min.unwrap_or(1.0);
865
866        // Get PASEF meta for this frame - these define which precursors were SELECTED
867        let Some(pasef_meta) = self.transmission_settings.pasef_meta.get(&(frame_id as i32)) else {
868            return TimsFrameAnnotated::new(frame_id as i32, rt, ms_type, vec![], vec![], vec![], vec![], vec![], vec![]);
869        };
870
871        // Preallocate with estimated capacity
872        let estimated_capacity = pasef_meta.len() * 4;
873        let mut tims_spectra: Vec<TimsSpectrumAnnotated> = Vec::with_capacity(estimated_capacity);
874
875        // Iterate over PASEF meta entries - each represents a SELECTED precursor
876        for meta in pasef_meta.iter() {
877            // Get the selected ion_id from the PASEF precursor field
878            let ion_id = meta.precursor as u32;
879
880            // Look up peptide_id and charge from the ion_id
881            let Some(&(peptide_id, charge_state)) = self.precursor_frame_builder.ion_id_to_peptide_charge.get(&ion_id) else {
882                continue;
883            };
884
885            // Get peptide info
886            let Some(peptide) = self.precursor_frame_builder.peptides.get(&peptide_id) else {
887                continue;
888            };
889
890            // Check if this peptide is active in this frame
891            if frame_id < peptide.frame_start || frame_id > peptide.frame_end {
892                continue;
893            }
894
895            // Get ion info from peptide_to_ions
896            let Some((ion_abundances, scan_occurrences, scan_abundances, charges, _)) = self
897                .precursor_frame_builder
898                .peptide_to_ions
899                .get(&peptide_id)
900            else {
901                continue;
902            };
903
904            // Find the index of this specific charge state
905            let Some(ion_index) = charges.iter().position(|&c| c == charge_state) else {
906                continue;
907            };
908
909            let ion_abundance = ion_abundances[ion_index];
910            let all_scan_occurrence = &scan_occurrences[ion_index];
911            let all_scan_abundance = &scan_abundances[ion_index];
912            let total_events = *self.precursor_frame_builder.peptide_to_events.get(&peptide_id).unwrap();
913
914            // Get frame abundance for this peptide
915            let frame_abundance = self.precursor_frame_builder.frame_to_abundances
916                .get(&frame_id)
917                .and_then(|(pep_ids, abundances)| {
918                    pep_ids.iter().position(|&p| p == peptide_id)
919                        .map(|idx| abundances[idx])
920                })
921                .unwrap_or(0.0);
922
923            if frame_abundance == 0.0 {
924                continue;
925            }
926
927            // Collision energy from the PASEF meta
928            let collision_energy = meta.collision_energy;
929            // Resolve the fragment CE key, tolerant to ~0.1 eV quantization noise
930            // (DDA pasef CE is full-precision; the stored key is f32-quantized).
931            let Some(collision_energy_quantized) = crate::sim::handle::resolve_fragment_ce_key(
932                fragment_ions, peptide_id, charge_state, collision_energy,
933            ) else {
934                // Fail loud if fragments exist for this precursor but none near the
935                // applied CE: the prediction set does not cover this instrument's CE
936                // (P5b). A precursor with NO predicted fragments at all is a
937                // legitimate skip. (For Bruker after the ±0.1 eV probe this never
938                // fires — verified 0 misses.)
939                if crate::sim::handle::fragment_prefix_exists(fragment_ions, peptide_id, charge_state) {
940                    panic!(
941                        "DDA fragment lookup miss: peptide {} charge {} applied CE {:.4} eV has \
942                         predicted fragments, but none within 0.1 eV — the prediction set does \
943                         not cover this instrument's collision energy",
944                        peptide_id, charge_state, collision_energy,
945                    );
946                }
947                continue;
948            };
949            let (_, fragment_series_vec) = fragment_ions
950                .get(&(peptide_id, charge_state, collision_energy_quantized))
951                .expect("resolve_fragment_ce_key returned a present key");
952
953            // Create ion for annotation and precursor spectrum calculation
954            let ion = PeptideIon::new(
955                peptide.sequence.sequence.clone(),
956                charge_state as i32,
957                ion_abundance as f64,
958                Some(peptide_id as i32),
959            );
960            // Calculate isotopic spectrum for precursor survival
961            let precursor_spectrum_annotated = ion.calculate_isotopic_spectrum_annotated(1e-3, 1e-8, 200, 1e-4);
962
963            // Process scans within the PASEF selection window
964            for (scan, scan_abundance) in all_scan_occurrence.iter().zip(all_scan_abundance.iter()) {
965                // Check if scan is within the PASEF selection window
966                let scan_i32 = *scan as i32;
967                if scan_i32 < meta.scan_start || scan_i32 > meta.scan_end {
968                    continue;
969                }
970
971                // Calculate abundance factor
972                let fraction_events = frame_abundance * scan_abundance * ion_abundance * total_events;
973
974                // Cache scan mobility
975                let scan_mobility = *self.precursor_frame_builder.scan_to_mobility.get(scan).unwrap() as f64;
976
977                for fragment_ion_series in fragment_series_vec.iter() {
978                    let scaled_spec = fragment_ion_series.clone() * fraction_events as f64;
979
980                    let mz_spectrum = if mz_noise_fragment {
981                        if uniform {
982                            scaled_spec.add_mz_noise_uniform(fragment_ppm, right_drag_val)
983                        } else {
984                            scaled_spec.add_mz_noise_normal(fragment_ppm)
985                        }
986                    } else {
987                        scaled_spec
988                    };
989
990                    let spectrum_len = mz_spectrum.mz.len();
991                    tims_spectra.push(TimsSpectrumAnnotated::new(
992                        frame_id as i32,
993                        *scan,
994                        rt,
995                        scan_mobility,
996                        ms_type.clone(),
997                        vec![0; spectrum_len],
998                        mz_spectrum,
999                    ));
1000                }
1001
1002                // Add unfragmented precursor ions (survival) if configured
1003                if self.isotope_transmission_config.has_precursor_survival() {
1004                    let mut rng = rand::thread_rng();
1005                    let survival_fraction = rng.gen_range(
1006                        self.isotope_transmission_config.precursor_survival_min
1007                        ..=self.isotope_transmission_config.precursor_survival_max
1008                    );
1009
1010                    if survival_fraction > 0.0 {
1011                        // Create a non-annotated spectrum for transmission
1012                        let precursor_mz_spectrum = MzSpectrum::new(
1013                            precursor_spectrum_annotated.mz.clone(),
1014                            precursor_spectrum_annotated.intensity.clone(),
1015                        );
1016
1017                        // Transmit through the quadrupole
1018                        let precursor_transmitted = self.transmission_settings.transmit_spectrum(
1019                            frame_id as i32,
1020                            scan_i32,
1021                            precursor_mz_spectrum,
1022                            Some(self.isotope_transmission_config.min_probability),
1023                        );
1024
1025                        if !precursor_transmitted.mz.is_empty() {
1026                            // Scale by survival fraction and event count
1027                            let precursor_scaled = precursor_transmitted * (fraction_events as f64 * survival_fraction);
1028
1029                            let precursor_final = if mz_noise_fragment {
1030                                if uniform {
1031                                    precursor_scaled.add_mz_noise_uniform(fragment_ppm, right_drag_val)
1032                                } else {
1033                                    precursor_scaled.add_mz_noise_normal(fragment_ppm)
1034                                }
1035                            } else {
1036                                precursor_scaled
1037                            };
1038
1039                            // Convert to annotated spectrum (with precursor annotations)
1040                            let annotations: Vec<PeakAnnotation> = precursor_final.mz.iter()
1041                                .map(|_| PeakAnnotation { contributions: vec![] })
1042                                .collect();
1043                            let precursor_annotated = MzSpectrumAnnotated::new(
1044                                precursor_final.mz.to_vec(),
1045                                precursor_final.intensity.to_vec(),
1046                                annotations,
1047                            );
1048
1049                            let precursor_len = precursor_annotated.mz.len();
1050                            tims_spectra.push(TimsSpectrumAnnotated::new(
1051                                frame_id as i32,
1052                                *scan,
1053                                rt,
1054                                scan_mobility,
1055                                ms_type.clone(),
1056                                vec![0; precursor_len],
1057                                precursor_annotated,
1058                            ));
1059                        }
1060                    }
1061                }
1062            }
1063        }
1064
1065        if tims_spectra.is_empty() {
1066            return TimsFrameAnnotated::new(frame_id as i32, rt, ms_type, vec![], vec![], vec![], vec![], vec![], vec![]);
1067        }
1068
1069        TimsFrameAnnotated::from_tims_spectra_annotated(tims_spectra).filter_ranged(
1070            mz_min_val, mz_max_val, 0.0, 10.0, 0, 1000, intensity_min_val, 1e9,
1071        )
1072    }
1073
1074    pub fn get_collision_energy(&self, frame_id: i32, scan_id: i32) -> f64 {
1075        self.transmission_settings.get_collision_energy(frame_id, scan_id).unwrap_or(0.0)
1076    }
1077
1078    pub fn get_collision_energies(&self, frame_ids: Vec<i32>, scan_ids: Vec<i32>) -> Vec<f64> {
1079        let mut collision_energies: Vec<f64> = Vec::new();
1080        for frame_id in frame_ids {
1081            for scan_id in &scan_ids {
1082                collision_energies.push(self.get_collision_energy(frame_id, *scan_id));
1083            }
1084        }
1085        collision_energies
1086    }
1087}