rustdf/sim/

dda.rs

use mscore::simulation::noise_rng::noise_rng;
use mscore::algorithm::isotope::calculate_precursor_transmission_factor;
use crate::sim::containers::IsotopeTransmissionMode;
use mscore::data::peptide::{PeptideIon, PeptideProductIonSeriesCollection};
use mscore::data::spectrum::{IndexedMzSpectrum, MsType, MzSpectrum};
use mscore::simulation::annotation::{
    MzSpectrumAnnotated, PeakAnnotation, TimsFrameAnnotated, TimsSpectrumAnnotated,
};
use mscore::timstof::frame::TimsFrame;
use mscore::timstof::quadrupole::{IonTransmission, TimsTransmissionDDA};
use mscore::timstof::spectrum::TimsSpectrum;
use std::collections::{BTreeMap, HashSet};
use std::path::Path;
use std::sync::Arc;

use rand::Rng;
use rayon::prelude::*;
use crate::sim::containers::IsotopeTransmissionConfig;
use crate::sim::scheme::InstrumentCapabilities;
use crate::sim::handle::{FragmentIonsWithComplementary, TimsTofSyntheticsDataHandle};
use crate::sim::precursor::TimsTofSyntheticsPrecursorFrameBuilder;

pub struct TimsTofSyntheticsFrameBuilderDDA {
    pub path: String,
    pub precursor_frame_builder: TimsTofSyntheticsPrecursorFrameBuilder,
    pub transmission_settings: TimsTransmissionDDA,
    pub fragment_ions:
        Option<BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrum>)>>,
    pub fragment_ions_annotated: Option<
        BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrumAnnotated>)>,
    >,
    /// Configuration for quad-selection dependent isotope transmission (already
    /// gated by the instrument capabilities at construction — P5e).
    pub isotope_transmission_config: IsotopeTransmissionConfig,
    /// Fragment ions with complementary distribution data (only populated when isotope_transmission_config.enabled)
    pub fragment_ions_with_complementary: Option<BTreeMap<(u32, i8, i32), FragmentIonsWithComplementary>>,
    /// Physical instrument capabilities (P5e). Default = Bruker timsTOF.
    pub capabilities: InstrumentCapabilities,
    /// Master seed for the simulation's noise. 0 keeps the legacy `thread_rng` behaviour
    /// for API users that never set it; TimSim sets it from the run's `sample_seed`.
    pub noise_seed: u64,
}

impl TimsTofSyntheticsFrameBuilderDDA {

    /// Set the master seed for the m/z jitter, making a run's noise reproducible.
    pub fn set_noise_seed(&mut self, seed: u64) {
        self.noise_seed = seed;
    }
    /// Create a new DDA frame builder.
    ///
    /// # Arguments
    ///
    /// * `path` - Path to the simulation database
    /// * `with_annotations` - Whether to include annotations
    /// * `num_threads` - Number of threads for parallel processing
    /// * `isotope_config` - Optional configuration for quad-dependent isotope transmission
    pub fn new(
        path: &Path,
        with_annotations: bool,
        num_threads: usize,
        isotope_config: Option<IsotopeTransmissionConfig>,
    ) -> Self {
        Self::new_with_capabilities(
            path,
            with_annotations,
            num_threads,
            isotope_config,
            InstrumentCapabilities::default(),
        )
    }

    /// Like [`Self::new`], but with explicit instrument capabilities (P5e). The
    /// isotope-transmission config is gated by them (an instrument without
    /// mobility-dependent quad isotope transmission forces the mode to None).
    /// Default capabilities = Bruker timsTOF, so [`Self::new`] is unchanged.
    pub fn new_with_capabilities(
        path: &Path,
        with_annotations: bool,
        num_threads: usize,
        isotope_config: Option<IsotopeTransmissionConfig>,
        capabilities: InstrumentCapabilities,
    ) -> Self {
        let handle = TimsTofSyntheticsDataHandle::new(path).unwrap();
        // P5b: refuse to render fragments stored under an incompatible prediction
        // set (e.g. a future Thermo set whose CE encoding the render keying can't
        // resolve). Bruker/legacy sets pass; this never fires for current DBs.
        handle
            .read_prediction_set()
            .expect("read prediction set")
            .assert_render_compatible()
            .expect("incompatible fragment prediction set for this renderer");
        let fragment_ions_raw = handle.read_fragment_ions().unwrap();
        let transmission_settings = handle.get_transmission_dda();

        let synthetics = TimsTofSyntheticsPrecursorFrameBuilder::new(path).unwrap();
        let config = isotope_config.unwrap_or_default().gated_by(capabilities);

        // Build fragment ions with complementary data if any transmission mode is enabled
        let fragment_ions_with_complementary = if config.is_enabled() {
            Some(TimsTofSyntheticsDataHandle::build_fragment_ions_with_transmission_data(
                &synthetics.peptides,
                &fragment_ions_raw,
                num_threads,
            ))
        } else {
            None
        };

        match with_annotations {
            true => {
                let fragment_ions =
                    Some(TimsTofSyntheticsDataHandle::build_fragment_ions_annotated(
                        &synthetics.peptides,
                        &fragment_ions_raw,
                        num_threads,
                    ));
                Self {
                    noise_seed: 0,
                    path: path.to_str().unwrap().to_string(),
                    precursor_frame_builder: synthetics,
                    transmission_settings,
                    fragment_ions: None,
                    fragment_ions_annotated: fragment_ions,
                    isotope_transmission_config: config,
                    fragment_ions_with_complementary,
                    capabilities,
                }
            }
            false => {
                let fragment_ions = Some(TimsTofSyntheticsDataHandle::build_fragment_ions(
                    &synthetics.peptides,
                    &fragment_ions_raw,
                    num_threads,
                ));
                Self {
                    noise_seed: 0,
                    path: path.to_str().unwrap().to_string(),
                    precursor_frame_builder: synthetics,
                    transmission_settings,
                    fragment_ions,
                    fragment_ions_annotated: None,
                    isotope_transmission_config: config,
                    fragment_ions_with_complementary,
                    capabilities,
                }
            }
        }
    }
    /// Construct an eager DDA frame builder from already-loaded, in-memory
    /// entities instead of reading the whole database.
    ///
    /// Used by the lazy DDA builder so it can delegate per-batch frame
    /// construction to the *same* (correct) eager algorithm instead of keeping a
    /// second, divergent copy. Only the non-annotated path is built (lazy loading
    /// does not support annotations). `fragment_ions_raw` is the per-batch slice
    /// of fragment ions for the supplied peptides, so memory stays bounded to the
    /// batch.
    pub fn from_entities(
        precursor_frame_builder: TimsTofSyntheticsPrecursorFrameBuilder,
        transmission_settings: TimsTransmissionDDA,
        fragment_ions_raw: Vec<crate::sim::containers::FragmentIonSim>,
        isotope_config: Option<IsotopeTransmissionConfig>,
        fragmentation: bool,
        num_threads: usize,
    ) -> Self {
        // Lazy DDA (delegates here) is Bruker timsTOF; gate is a no-op. P6 will
        // thread real capabilities through the lazy path.
        let config = isotope_config
            .unwrap_or_default()
            .gated_by(InstrumentCapabilities::default());

        // The fragment isotope map is only consumed when fragmentation is on
        // (build_ms2_frame's `true` branch). Expanding all predicted fragment
        // intensities into isotope spectra is expensive, so skip it entirely in
        // no-fragmentation mode (the `false` branch never touches these maps).
        let (fragment_ions, fragment_ions_with_complementary) = if fragmentation {
            let with_complementary = if config.is_enabled() {
                Some(TimsTofSyntheticsDataHandle::build_fragment_ions_with_transmission_data(
                    &precursor_frame_builder.peptides,
                    &fragment_ions_raw,
                    num_threads,
                ))
            } else {
                None
            };
            let fragment_ions = Some(TimsTofSyntheticsDataHandle::build_fragment_ions(
                &precursor_frame_builder.peptides,
                &fragment_ions_raw,
                num_threads,
            ));
            (fragment_ions, with_complementary)
        } else {
            (None, None)
        };

        Self {
            noise_seed: 0,
            path: String::new(),
            precursor_frame_builder,
            transmission_settings,
            fragment_ions,
            fragment_ions_annotated: None,
            isotope_transmission_config: config,
            fragment_ions_with_complementary,
            capabilities: InstrumentCapabilities::default(),
        }
    }

    /// Build a frame for DDA synthetic experiment
    ///
    /// # Arguments
    ///
    /// * `frame_id` - The frame id
    /// * `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
    ///
    /// # Returns
    ///
    /// A TimsFrame
    ///
    pub fn build_frame(
        &self,
        frame_id: u32,
        fragmentation: bool,
        mz_noise_precursor: bool,
        uniform: bool,
        precursor_noise_ppm: f64,
        mz_noise_fragment: bool,
        fragment_noise_ppm: f64,
        right_drag: bool,
    ) -> TimsFrame {
        // determine if the frame is a precursor frame
        match self
            .precursor_frame_builder
            .precursor_frame_id_set
            .contains(&frame_id)
        {
            true => self.build_ms1_frame(
                frame_id,
                mz_noise_precursor,
                uniform,
                precursor_noise_ppm,
                right_drag,
            ),
            false => self.build_ms2_frame(
                frame_id,
                fragmentation,
                mz_noise_fragment,
                uniform,
                fragment_noise_ppm,
                right_drag,
            ),
        }
    }

    pub fn build_frame_annotated(
        &self,
        frame_id: u32,
        fragmentation: bool,
        mz_noise_precursor: bool,
        uniform: bool,
        precursor_noise_ppm: f64,
        mz_noise_fragment: bool,
        fragment_noise_ppm: f64,
        right_drag: bool,
    ) -> TimsFrameAnnotated {
        match self
            .precursor_frame_builder
            .precursor_frame_id_set
            .contains(&frame_id)
        {
            true => self.build_ms1_frame_annotated(
                frame_id,
                mz_noise_precursor,
                uniform,
                precursor_noise_ppm,
                right_drag,
            ),
            false => self.build_ms2_frame_annotated(
                frame_id,
                fragmentation,
                mz_noise_fragment,
                uniform,
                fragment_noise_ppm,
                right_drag,
            ),
        }
    }

    pub fn get_fragment_ion_ids(&self, precursor_frame_ids: Vec<u32>) -> Vec<u32> {
        let mut peptide_ids: HashSet<u32> = HashSet::new();
        // get all peptide ids for the precursor frame ids
        for frame_id in precursor_frame_ids {
            for (peptide_id, peptide) in self.precursor_frame_builder.peptides.iter() {
                if peptide.frame_start <= frame_id && peptide.frame_end >= frame_id {
                    peptide_ids.insert(*peptide_id);
                }
            }
        }
        // get all ion ids for the peptide ids
        let mut result: Vec<u32> = Vec::new();
        for item in peptide_ids {
            let ions = self.precursor_frame_builder.ions.get(&item).unwrap();
            for ion in ions.iter() {
                result.push(ion.ion_id);
            }
        }
        result
    }

    pub fn build_frames(
        &self,
        frame_ids: Vec<u32>,
        fragmentation: bool,
        mz_noise_precursor: bool,
        uniform: bool,
        precursor_noise_ppm: f64,
        mz_noise_fragment: bool,
        fragment_noise_ppm: f64,
        right_drag: bool,
        num_threads: usize,
    ) -> Vec<TimsFrame> {
        // Use thread pool with custom parallelism
        let pool = rayon::ThreadPoolBuilder::new()
            .num_threads(num_threads)
            .build()
            .unwrap();

        pool.install(|| {
            // Use indexed parallel iteration to maintain order, avoiding post-sort
            let mut tims_frames: Vec<TimsFrame> = Vec::with_capacity(frame_ids.len());
            unsafe { tims_frames.set_len(frame_ids.len()); }

            frame_ids.par_iter().enumerate().for_each(|(idx, frame_id)| {
                let frame = self.build_frame(
                    *frame_id,
                    fragmentation,
                    mz_noise_precursor,
                    uniform,
                    precursor_noise_ppm,
                    mz_noise_fragment,
                    fragment_noise_ppm,
                    right_drag,
                );
                unsafe {
                    let ptr = tims_frames.as_ptr() as *mut TimsFrame;
                    std::ptr::write(ptr.add(idx), frame);
                }
            });

            tims_frames
        })
    }
    pub fn build_frames_annotated(
        &self,
        frame_ids: Vec<u32>,
        fragmentation: bool,
        mz_noise_precursor: bool,
        uniform: bool,
        precursor_noise_ppm: f64,
        mz_noise_fragment: bool,
        fragment_noise_ppm: f64,
        right_drag: bool,
        num_threads: usize,
    ) -> Vec<TimsFrameAnnotated> {
        // Use thread pool with custom parallelism
        let pool = rayon::ThreadPoolBuilder::new()
            .num_threads(num_threads)
            .build()
            .unwrap();

        pool.install(|| {
            // Use indexed parallel iteration to maintain order, avoiding post-sort
            let mut tims_frames: Vec<TimsFrameAnnotated> = Vec::with_capacity(frame_ids.len());
            unsafe { tims_frames.set_len(frame_ids.len()); }

            frame_ids.par_iter().enumerate().for_each(|(idx, frame_id)| {
                let frame = self.build_frame_annotated(
                    *frame_id,
                    fragmentation,
                    mz_noise_precursor,
                    uniform,
                    precursor_noise_ppm,
                    mz_noise_fragment,
                    fragment_noise_ppm,
                    right_drag,
                );
                unsafe {
                    let ptr = tims_frames.as_ptr() as *mut TimsFrameAnnotated;
                    std::ptr::write(ptr.add(idx), frame);
                }
            });

            tims_frames
        })
    }

    fn build_ms1_frame(
        &self,
        frame_id: u32,
        mz_noise_precursor: bool,
        uniform: bool,
        precursor_ppm: f64,
        right_drag: bool,
    ) -> TimsFrame {
        let mut tims_frame = self.precursor_frame_builder.build_precursor_frame(
            frame_id,
            mz_noise_precursor,
            uniform,
            precursor_ppm,
            right_drag,
        );
        let intensities_rounded = tims_frame
            .ims_frame
            .intensity
            .iter()
            .map(|x| x.round())
            .collect::<Vec<_>>();
        tims_frame.ims_frame.intensity = Arc::new(intensities_rounded);
        tims_frame
    }

    fn build_ms1_frame_annotated(
        &self,
        frame_id: u32,
        mz_noise_precursor: bool,
        uniform: bool,
        precursor_ppm: f64,
        right_drag: bool,
    ) -> TimsFrameAnnotated {
        let mut tims_frame = self
            .precursor_frame_builder
            .build_precursor_frame_annotated(
                frame_id,
                mz_noise_precursor,
                uniform,
                precursor_ppm,
                right_drag,
            );
        let intensities_rounded = tims_frame
            .intensity
            .iter()
            .map(|x| x.round())
            .collect::<Vec<_>>();
        tims_frame.intensity = intensities_rounded;
        tims_frame
    }

    fn build_ms2_frame(
        &self,
        frame_id: u32,
        fragmentation: bool,
        mz_noise_fragment: bool,
        uniform: bool,
        fragment_ppm: f64,
        right_drag: bool,
    ) -> TimsFrame {
        match fragmentation {
            false => {
                let mut frame = self.transmission_settings.transmit_tims_frame(
                    &self.build_ms1_frame(
                        frame_id,
                        mz_noise_fragment,
                        uniform,
                        fragment_ppm,
                        right_drag,
                    ),
                    None,
                );
                let intensities_rounded = frame
                    .ims_frame
                    .intensity
                    .iter()
                    .map(|x| x.round())
                    .collect::<Vec<_>>();
                frame.ims_frame.intensity = Arc::new(intensities_rounded);
                // DDA MS2 frames are PASEF fragment frames (MsMsType 8). The
                // no-fragmentation mode still produces a DDA fragment frame
                // (quad-filtered precursor), so it must keep the DDA type, not
                // FragmentDia (9). The fragmentation=true branch already uses
                // FragmentDda; this matches it.
                frame.ms_type = MsType::FragmentDda;
                frame
            }
            true => {
                let mut frame = self.build_fragment_frame(
                    frame_id,
                    &self.fragment_ions.as_ref().unwrap(),
                    mz_noise_fragment,
                    uniform,
                    fragment_ppm,
                    None,
                    None,
                    None,
                    Some(right_drag),
                );
                let intensities_rounded = frame
                    .ims_frame
                    .intensity
                    .iter()
                    .map(|x| x.round())
                    .collect::<Vec<_>>();
                frame.ims_frame.intensity = Arc::new(intensities_rounded);
                frame
            }
        }
    }

    fn build_ms2_frame_annotated(
        &self,
        frame_id: u32,
        fragmentation: bool,
        mz_noise_fragment: bool,
        uniform: bool,
        fragment_ppm: f64,
        right_drag: bool,
    ) -> TimsFrameAnnotated {
        match fragmentation {
            false => {
                let mut frame = self.transmission_settings.transmit_tims_frame_annotated(
                    &self.build_ms1_frame_annotated(
                        frame_id,
                        mz_noise_fragment,
                        uniform,
                        fragment_ppm,
                        right_drag,
                    ),
                    None,
                );
                let intensities_rounded = frame
                    .intensity
                    .iter()
                    .map(|x| x.round())
                    .collect::<Vec<_>>();
                frame.intensity = intensities_rounded;
                // See build_ms2_frame: DDA no-frag MS2 stays FragmentDda (8).
                frame.ms_type = MsType::FragmentDda;
                frame
            }
            true => {
                let mut frame = self.build_fragment_frame_annotated(
                    frame_id,
                    &self.fragment_ions_annotated.as_ref().unwrap(),
                    mz_noise_fragment,
                    uniform,
                    fragment_ppm,
                    None,
                    None,
                    None,
                    Some(right_drag),
                );
                let intensities_rounded = frame
                    .intensity
                    .iter()
                    .map(|x| x.round())
                    .collect::<Vec<_>>();
                frame.intensity = intensities_rounded;
                frame
            }
        }
    }

    /// Build a fragment frame
    ///
    /// # Arguments
    ///
    /// * `frame_id` - The frame id
    /// * `mz_min` - The minimum m/z value in fragment spectrum
    /// * `mz_max` - The maximum m/z value in fragment spectrum
    /// * `intensity_min` - The minimum intensity value in fragment spectrum
    ///
    /// # Returns
    ///
    /// A TimsFrame
    ///
    fn build_fragment_frame(
        &self,
        frame_id: u32,
        fragment_ions: &BTreeMap<
            (u32, i8, i32),
            (PeptideProductIonSeriesCollection, Vec<MzSpectrum>),
        >,
        mz_noise_fragment: bool,
        uniform: bool,
        fragment_ppm: f64,
        mz_min: Option<f64>,
        mz_max: Option<f64>,
        intensity_min: Option<f64>,
        right_drag: Option<bool>,
    ) -> TimsFrame {
        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
        // Cache frame-level lookups
        let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
            MsType::Unknown
        } else {
            MsType::FragmentDda
        };

        let rt = *self.precursor_frame_builder.frame_to_rt.get(&frame_id)
            .expect("frame_to_rt should always have this frame") as f64;
        let right_drag_val = right_drag.unwrap_or(false);
        let mz_min_val = mz_min.unwrap_or(100.0);
        let mz_max_val = mz_max.unwrap_or(1700.0);
        let intensity_min_val = intensity_min.unwrap_or(1.0);

        // Get PASEF meta for this frame - these define which precursors were SELECTED
        let Some(pasef_meta) = self.transmission_settings.pasef_meta.get(&(frame_id as i32)) else {
            return TimsFrame::new(frame_id as i32, ms_type, rt, vec![], vec![], vec![], vec![], vec![]);
        };

        // Preallocate with estimated capacity
        let estimated_capacity = pasef_meta.len() * 4;
        let mut tims_spectra: Vec<TimsSpectrum> = Vec::with_capacity(estimated_capacity);

        // Iterate over PASEF meta entries - each represents a SELECTED precursor
        for meta in pasef_meta.iter() {
            // Get the selected ion_id from the PASEF precursor field
            let ion_id = meta.precursor as u32;

            // Look up peptide_id and charge from the ion_id
            let Some(&(peptide_id, charge_state)) = self.precursor_frame_builder.ion_id_to_peptide_charge.get(&ion_id) else {
                continue;
            };

            // Get peptide info
            let Some(peptide) = self.precursor_frame_builder.peptides.get(&peptide_id) else {
                continue;
            };

            // Check if this peptide is active in this frame
            if frame_id < peptide.frame_start || frame_id > peptide.frame_end {
                continue;
            }

            // Get ion info from peptide_to_ions
            let Some((ion_abundances, scan_occurrences, scan_abundances, charges, spectra)) = self
                .precursor_frame_builder
                .peptide_to_ions
                .get(&peptide_id)
            else {
                continue;
            };

            // Find the index of this specific charge state
            let Some(ion_index) = charges.iter().position(|&c| c == charge_state) else {
                continue;
            };

            let ion_abundance = ion_abundances[ion_index];
            let all_scan_occurrence = &scan_occurrences[ion_index];
            let all_scan_abundance = &scan_abundances[ion_index];
            let spectrum = &spectra[ion_index];
            let total_events = *self.precursor_frame_builder.peptide_to_events.get(&peptide_id).unwrap();

            // Get frame abundance for this peptide
            let frame_abundance = self.precursor_frame_builder.frame_to_abundances
                .get(&frame_id)
                .and_then(|(pep_ids, abundances)| {
                    pep_ids.iter().position(|&p| p == peptide_id)
                        .map(|idx| abundances[idx])
                })
                .unwrap_or(0.0);

            if frame_abundance == 0.0 {
                continue;
            }

            // Collision energy from the PASEF meta
            let collision_energy = meta.collision_energy;
            // Resolve the fragment CE key, tolerant to ~0.1 eV quantization noise
            // (DDA pasef CE is full-precision; the stored key is f32-quantized).
            let Some(collision_energy_quantized) = crate::sim::handle::resolve_fragment_ce_key(
                fragment_ions, peptide_id, charge_state, collision_energy,
            ) else {
                // Fail loud if fragments exist for this precursor but none near the
                // applied CE: the prediction set does not cover this instrument's CE
                // (P5b). A precursor with NO predicted fragments at all is a
                // legitimate skip. (For Bruker after the ±0.1 eV probe this never
                // fires — verified 0 misses.)
                if crate::sim::handle::fragment_prefix_exists(fragment_ions, peptide_id, charge_state) {
                    panic!(
                        "DDA fragment lookup miss: peptide {} charge {} applied CE {:.4} eV has \
                         predicted fragments, but none within 0.1 eV — the prediction set does \
                         not cover this instrument's collision energy",
                        peptide_id, charge_state, collision_energy,
                    );
                }
                continue;
            };
            let (_, fragment_series_vec) = fragment_ions
                .get(&(peptide_id, charge_state, collision_energy_quantized))
                .expect("resolve_fragment_ce_key returned a present key");

            // Process scans within the PASEF selection window
            for (scan, scan_abundance) in all_scan_occurrence.iter().zip(all_scan_abundance.iter()) {
                // Check if scan is within the PASEF selection window
                let scan_i32 = *scan as i32;
                if scan_i32 < meta.scan_start || scan_i32 > meta.scan_end {
                    continue;
                }

                // Get transmitted isotope indices based on config mode
                let transmitted_indices = match self.isotope_transmission_config.mode {
                    IsotopeTransmissionMode::None => HashSet::new(),
                    IsotopeTransmissionMode::PrecursorScaling | IsotopeTransmissionMode::PerFragment => {
                        self.transmission_settings.get_transmission_set(
                            frame_id as i32,
                            scan_i32,
                            &spectrum.mz,
                            Some(self.isotope_transmission_config.min_probability),
                        )
                    },
                };

                // Calculate abundance factor
                let fraction_events = frame_abundance * scan_abundance * ion_abundance * total_events;

                // Cache scan mobility
                let scan_mobility = *self.precursor_frame_builder.scan_to_mobility.get(scan).unwrap() as f64;

                // Calculate transmission factor for PrecursorScaling mode
                let transmission_factor = if self.isotope_transmission_config.mode == IsotopeTransmissionMode::PrecursorScaling {
                    if let Some(comp_data) = self.fragment_ions_with_complementary.as_ref() {
                        if let Some(frag_data) = comp_data.get(&(peptide_id, charge_state, collision_energy_quantized)) {
                            calculate_precursor_transmission_factor(
                                &frag_data.precursor_isotope_distribution,
                                &transmitted_indices,
                            )
                        } else {
                            1.0
                        }
                    } else {
                        1.0
                    }
                } else {
                    1.0
                };

                // Complementary per-(peptide,charge,CE) data for PerFragment;
                // looked up once (same for all series) and passed to the shared
                // kernel — same lookup the transmission_factor block above uses.
                let frag_data = self
                    .fragment_ions_with_complementary
                    .as_ref()
                    .and_then(|c| c.get(&(peptide_id, charge_state, collision_energy_quantized)));

                for (series_idx, fragment_ion_series) in fragment_series_vec.iter().enumerate() {
                    let final_spectrum = crate::sim::dia::fragment_series_spectrum(
                        self.isotope_transmission_config.mode,
                        fragment_ion_series,
                        series_idx,
                        fraction_events,
                        transmission_factor,
                        frag_data,
                        &transmitted_indices,
                        self.isotope_transmission_config.max_isotopes,
                    );

                    let mz_spectrum = if mz_noise_fragment {
                        if uniform {
                            final_spectrum.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                        } else {
                            final_spectrum.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                        }
                    } else {
                        final_spectrum
                    };

                    let spectrum_len = mz_spectrum.mz.len();
                    tims_spectra.push(TimsSpectrum::new(
                        frame_id as i32,
                        *scan as i32,
                        rt,
                        scan_mobility,
                        ms_type.clone(),
                        IndexedMzSpectrum::from_mz_spectrum(
                            vec![0; spectrum_len],
                            mz_spectrum,
                        ).filter_ranged(100.0, 1700.0, 1.0, 1e9),
                    ));
                }

                // Add unfragmented precursor ions (survival) if configured
                if self.isotope_transmission_config.has_precursor_survival() {
                    let mut rng = rand::thread_rng();
                    let survival_fraction = rng.gen_range(
                        self.isotope_transmission_config.precursor_survival_min
                        ..=self.isotope_transmission_config.precursor_survival_max
                    );

                    if survival_fraction > 0.0 {
                        // Transmit the precursor spectrum through the quadrupole
                        let precursor_transmitted = self.transmission_settings.transmit_spectrum(
                            frame_id as i32,
                            *scan as i32,
                            spectrum.clone(),
                            Some(self.isotope_transmission_config.min_probability),
                        );

                        if !precursor_transmitted.mz.is_empty() {
                            // Scale by survival fraction and event count
                            let precursor_scaled = precursor_transmitted * (fraction_events as f64 * survival_fraction);

                            let precursor_mz_spectrum = if mz_noise_fragment {
                                if uniform {
                                    precursor_scaled.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                                } else {
                                    precursor_scaled.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                                }
                            } else {
                                precursor_scaled
                            };

                            let precursor_len = precursor_mz_spectrum.mz.len();
                            tims_spectra.push(TimsSpectrum::new(
                                frame_id as i32,
                                *scan as i32,
                                rt,
                                scan_mobility,
                                ms_type.clone(),
                                IndexedMzSpectrum::from_mz_spectrum(
                                    vec![0; precursor_len],
                                    precursor_mz_spectrum,
                                ).filter_ranged(100.0, 1700.0, 1.0, 1e9),
                            ));
                        }
                    }
                }
            }
        }

        if tims_spectra.is_empty() {
            return TimsFrame::new(frame_id as i32, ms_type, rt, vec![], vec![], vec![], vec![], vec![]);
        }

        let tims_frame = TimsFrame::from_tims_spectra(tims_spectra);
        tims_frame.filter_ranged(
            mz_min_val,
            mz_max_val,
            0,
            1000,
            0.0,
            10.0,
            intensity_min_val,
            1e9,
            0,
            i32::MAX,
        )
    }

    pub fn build_fragment_frame_annotated(
        &self,
        frame_id: u32,
        fragment_ions: &BTreeMap<
            (u32, i8, i32),
            (PeptideProductIonSeriesCollection, Vec<MzSpectrumAnnotated>),
        >,
        mz_noise_fragment: bool,
        uniform: bool,
        fragment_ppm: f64,
        mz_min: Option<f64>,
        mz_max: Option<f64>,
        intensity_min: Option<f64>,
        right_drag: Option<bool>,
    ) -> TimsFrameAnnotated {
        // One RNG per frame, keyed by the master seed and the frame id, so the m/z jitter
        // does not depend on which thread builds the frame. See mscore::simulation::noise_rng.
        let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
        // Cache frame-level lookups
        let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
            MsType::Unknown
        } else {
            MsType::FragmentDda
        };

        let rt = *self.precursor_frame_builder.frame_to_rt.get(&frame_id).unwrap() as f64;
        let right_drag_val = right_drag.unwrap_or(false);
        let mz_min_val = mz_min.unwrap_or(100.0);
        let mz_max_val = mz_max.unwrap_or(1700.0);
        let intensity_min_val = intensity_min.unwrap_or(1.0);

        // Get PASEF meta for this frame - these define which precursors were SELECTED
        let Some(pasef_meta) = self.transmission_settings.pasef_meta.get(&(frame_id as i32)) else {
            return TimsFrameAnnotated::new(frame_id as i32, rt, ms_type, vec![], vec![], vec![], vec![], vec![], vec![]);
        };

        // Preallocate with estimated capacity
        let estimated_capacity = pasef_meta.len() * 4;
        let mut tims_spectra: Vec<TimsSpectrumAnnotated> = Vec::with_capacity(estimated_capacity);

        // Iterate over PASEF meta entries - each represents a SELECTED precursor
        for meta in pasef_meta.iter() {
            // Get the selected ion_id from the PASEF precursor field
            let ion_id = meta.precursor as u32;

            // Look up peptide_id and charge from the ion_id
            let Some(&(peptide_id, charge_state)) = self.precursor_frame_builder.ion_id_to_peptide_charge.get(&ion_id) else {
                continue;
            };

            // Get peptide info
            let Some(peptide) = self.precursor_frame_builder.peptides.get(&peptide_id) else {
                continue;
            };

            // Check if this peptide is active in this frame
            if frame_id < peptide.frame_start || frame_id > peptide.frame_end {
                continue;
            }

            // Get ion info from peptide_to_ions
            let Some((ion_abundances, scan_occurrences, scan_abundances, charges, _)) = self
                .precursor_frame_builder
                .peptide_to_ions
                .get(&peptide_id)
            else {
                continue;
            };

            // Find the index of this specific charge state
            let Some(ion_index) = charges.iter().position(|&c| c == charge_state) else {
                continue;
            };

            let ion_abundance = ion_abundances[ion_index];
            let all_scan_occurrence = &scan_occurrences[ion_index];
            let all_scan_abundance = &scan_abundances[ion_index];
            let total_events = *self.precursor_frame_builder.peptide_to_events.get(&peptide_id).unwrap();

            // Get frame abundance for this peptide
            let frame_abundance = self.precursor_frame_builder.frame_to_abundances
                .get(&frame_id)
                .and_then(|(pep_ids, abundances)| {
                    pep_ids.iter().position(|&p| p == peptide_id)
                        .map(|idx| abundances[idx])
                })
                .unwrap_or(0.0);

            if frame_abundance == 0.0 {
                continue;
            }

            // Collision energy from the PASEF meta
            let collision_energy = meta.collision_energy;
            // Resolve the fragment CE key, tolerant to ~0.1 eV quantization noise
            // (DDA pasef CE is full-precision; the stored key is f32-quantized).
            let Some(collision_energy_quantized) = crate::sim::handle::resolve_fragment_ce_key(
                fragment_ions, peptide_id, charge_state, collision_energy,
            ) else {
                // Fail loud if fragments exist for this precursor but none near the
                // applied CE: the prediction set does not cover this instrument's CE
                // (P5b). A precursor with NO predicted fragments at all is a
                // legitimate skip. (For Bruker after the ±0.1 eV probe this never
                // fires — verified 0 misses.)
                if crate::sim::handle::fragment_prefix_exists(fragment_ions, peptide_id, charge_state) {
                    panic!(
                        "DDA fragment lookup miss: peptide {} charge {} applied CE {:.4} eV has \
                         predicted fragments, but none within 0.1 eV — the prediction set does \
                         not cover this instrument's collision energy",
                        peptide_id, charge_state, collision_energy,
                    );
                }
                continue;
            };
            let (_, fragment_series_vec) = fragment_ions
                .get(&(peptide_id, charge_state, collision_energy_quantized))
                .expect("resolve_fragment_ce_key returned a present key");

            // Create ion for annotation and precursor spectrum calculation
            let ion = PeptideIon::new(
                peptide.sequence.sequence.clone(),
                charge_state as i32,
                ion_abundance as f64,
                Some(peptide_id as i32),
            );
            // Calculate isotopic spectrum for precursor survival
            let precursor_spectrum_annotated = ion.calculate_isotopic_spectrum_annotated(1e-3, 1e-8, 200, 1e-4);

            // Process scans within the PASEF selection window
            for (scan, scan_abundance) in all_scan_occurrence.iter().zip(all_scan_abundance.iter()) {
                // Check if scan is within the PASEF selection window
                let scan_i32 = *scan as i32;
                if scan_i32 < meta.scan_start || scan_i32 > meta.scan_end {
                    continue;
                }

                // Calculate abundance factor
                let fraction_events = frame_abundance * scan_abundance * ion_abundance * total_events;

                // Cache scan mobility
                let scan_mobility = *self.precursor_frame_builder.scan_to_mobility.get(scan).unwrap() as f64;

                for fragment_ion_series in fragment_series_vec.iter() {
                    let scaled_spec = fragment_ion_series.clone() * fraction_events as f64;

                    let mz_spectrum = if mz_noise_fragment {
                        if uniform {
                            scaled_spec.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                        } else {
                            scaled_spec.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                        }
                    } else {
                        scaled_spec
                    };

                    let spectrum_len = mz_spectrum.mz.len();
                    tims_spectra.push(TimsSpectrumAnnotated::new(
                        frame_id as i32,
                        *scan,
                        rt,
                        scan_mobility,
                        ms_type.clone(),
                        vec![0; spectrum_len],
                        mz_spectrum,
                    ));
                }

                // Add unfragmented precursor ions (survival) if configured
                if self.isotope_transmission_config.has_precursor_survival() {
                    let mut rng = rand::thread_rng();
                    let survival_fraction = rng.gen_range(
                        self.isotope_transmission_config.precursor_survival_min
                        ..=self.isotope_transmission_config.precursor_survival_max
                    );

                    if survival_fraction > 0.0 {
                        // Create a non-annotated spectrum for transmission
                        let precursor_mz_spectrum = MzSpectrum::new(
                            precursor_spectrum_annotated.mz.clone(),
                            precursor_spectrum_annotated.intensity.clone(),
                        );

                        // Transmit through the quadrupole
                        let precursor_transmitted = self.transmission_settings.transmit_spectrum(
                            frame_id as i32,
                            scan_i32,
                            precursor_mz_spectrum,
                            Some(self.isotope_transmission_config.min_probability),
                        );

                        if !precursor_transmitted.mz.is_empty() {
                            // Scale by survival fraction and event count
                            let precursor_scaled = precursor_transmitted * (fraction_events as f64 * survival_fraction);

                            let precursor_final = if mz_noise_fragment {
                                if uniform {
                                    precursor_scaled.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
                                } else {
                                    precursor_scaled.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
                                }
                            } else {
                                precursor_scaled
                            };

                            // Convert to annotated spectrum (with precursor annotations)
                            let annotations: Vec<PeakAnnotation> = precursor_final.mz.iter()
                                .map(|_| PeakAnnotation { contributions: vec![] })
                                .collect();
                            let precursor_annotated = MzSpectrumAnnotated::new(
                                precursor_final.mz.to_vec(),
                                precursor_final.intensity.to_vec(),
                                annotations,
                            );

                            let precursor_len = precursor_annotated.mz.len();
                            tims_spectra.push(TimsSpectrumAnnotated::new(
                                frame_id as i32,
                                *scan,
                                rt,
                                scan_mobility,
                                ms_type.clone(),
                                vec![0; precursor_len],
                                precursor_annotated,
                            ));
                        }
                    }
                }
            }
        }

        if tims_spectra.is_empty() {
            return TimsFrameAnnotated::new(frame_id as i32, rt, ms_type, vec![], vec![], vec![], vec![], vec![], vec![]);
        }

        TimsFrameAnnotated::from_tims_spectra_annotated(tims_spectra).filter_ranged(
            mz_min_val, mz_max_val, 0.0, 10.0, 0, 1000, intensity_min_val, 1e9,
        )
    }

    pub fn get_collision_energy(&self, frame_id: i32, scan_id: i32) -> f64 {
        self.transmission_settings.get_collision_energy(frame_id, scan_id).unwrap_or(0.0)
    }

    pub fn get_collision_energies(&self, frame_ids: Vec<i32>, scan_ids: Vec<i32>) -> Vec<f64> {
        let mut collision_energies: Vec<f64> = Vec::new();
        for frame_id in frame_ids {
            for scan_id in &scan_ids {
                collision_energies.push(self.get_collision_energy(frame_id, *scan_id));
            }
        }
        collision_energies
    }
}