Skip to main content

lindhard/tally/
ion.rs

1//! The full ion tally: ranges, damage and escapes, wired to the BCA engine.
2//!
3//! [`IonTally`] implements [`BcaTally`]. Build one prototype with
4//! [`IonTally::new`] and hand clones to [`crate::ion::bca::Bca::run`]; the
5//! merged tally turns into a plain-data [`IonReport`] with
6//! [`IonTally::report`].
7//!
8//! Units: positions and lengths in metres, energies in eV, angles in radians.
9//! Totals are summed over all histories; divide by
10//! [`IonReport::histories`] for per-ion values (the escape coefficients are
11//! already per ion).
12//!
13//! # Determinism
14//!
15//! Counts are integers, moments use the mergeable accumulator of
16//! [`super::moments`], and energy sums are plain `f64` sums; all are merged
17//! field by field in the order [`crate::rng::run_particles`] calls `merge`
18//! (chunk order). For a fixed `chunk_size` the result is bit-identical at any
19//! thread count (`lindhard/tests/tally.rs`). Nothing is stored in a hash map.
20
21use serde::{Deserialize, Serialize};
22
23use super::hist::{Binning, Histogram, WeightedHistogram};
24use super::moments::{MomentSummary, Moments};
25use super::pearson::{DualPearson, DualPearsonFit, PearsonError, PearsonIv};
26use crate::geometry::Stack;
27use crate::ion::bca::{BcaTally, EnergyBudget, Face, LatticeDeposit, Particle};
28use crate::ion::damage::{damage_energy_ev, kinchin_pease, nrt_displacements};
29use crate::material::MaterialError;
30
31/// Binnings for an [`IonTally`].
32#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
33pub struct IonTallyConfig {
34    /// Depth `x`, m: range histograms (overall and per layer) and the
35    /// defect depth histograms.
36    pub depth: Binning,
37    /// Lateral position `y` and `z`, m, of stopped primaries.
38    pub lateral: Binning,
39    /// Radial distance `(y^2 + z^2)^(1/2)`, m, of stopped primaries.
40    pub radial: Binning,
41    /// Kinetic energy outside the target of escaping particles, eV.
42    pub escape_energy: Binning,
43    /// Polar angle of escaping particles from the outward surface normal,
44    /// rad (`0` is along the normal, `π/2` grazing).
45    pub escape_polar: Binning,
46}
47
48/// Displacement-model estimates (NRT and Kinchin-Pease) from the PKA damage
49/// energies. **Not** a count of simulated defects: see [`CascadeDefects`]
50/// and [`crate::ion::damage`] for why the two are kept apart.
51///
52/// In a layer with more than one element, the Lindhard partition uses the
53/// PKA's own `Z1`/`A1` and the layer's atom-fraction mean `Z2`/`A2`. NRT is
54/// defined for monatomic targets, so this is an approximation and this
55/// crate's own convention: see
56/// [the compound-target section](crate::ion::damage#compound-and-multi-element-targets-an-approximation)
57/// of [`crate::ion::damage`].
58#[derive(Debug, Clone, Copy, PartialEq, Default, Serialize, Deserialize)]
59pub struct NrtDamage {
60    /// Primary knock-on atoms: target atoms displaced by the beam particle.
61    pub pka_count: u64,
62    /// Sum of the PKA energies `T` (energy transferred, before the lattice
63    /// binding is subtracted), eV.
64    pub pka_energy_ev: f64,
65    /// Sum of the PKA damage energies `T_dam` (Lindhard partition), eV.
66    pub damage_energy_ev: f64,
67    /// Sum of the NRT displacements `N_NRT(T_dam)` over PKAs.
68    pub nrt_displacements: f64,
69    /// Sum of the Kinchin-Pease displacements `N_KP(T_dam)` over PKAs.
70    pub kinchin_pease_displacements: f64,
71}
72
73impl NrtDamage {
74    fn merge(&mut self, o: &NrtDamage) {
75        self.pka_count += o.pka_count;
76        self.pka_energy_ev += o.pka_energy_ev;
77        self.damage_energy_ev += o.damage_energy_ev;
78        self.nrt_displacements += o.nrt_displacements;
79        self.kinchin_pease_displacements += o.kinchin_pease_displacements;
80    }
81}
82
83/// Defects counted event by event in the simulated cascades.
84///
85/// * `displacements`: target atoms given more than `E_d` (recoils created).
86/// * `replacements`: a moving atom came to rest, immediately after the
87///   collision in which it displaced an atom of its own element, at that
88///   atom's site, so it fills the site. In this engine a particle stops when
89///   its energy falls below its cutoff, so this count depends on the cutoffs
90///   (a recoil cutoff near `E_d` gives the most replacements). When a recoil
91///   is not followed (below its cutoff, or with
92///   [`BcaConfig::follow_recoils`](crate::ion::bca::BcaConfig::follow_recoils)
93///   `= false`), the engine reports `stopped(&r)` straight after
94///   `recoil(&r)`, which clears the pending site, so the displacer cannot
95///   score a replacement in that collision. With `follow_recoils = false` the
96///   replacement count is therefore always 0.
97/// * `vacancies = displacements - replacements`.
98/// * `interstitials`: recoils that came to rest in the target and did not
99///   fill a site. Implanted beam particles are not counted here (see the
100///   range tally). Recoils that leave the target create a vacancy and no
101///   interstitial, so the two counts differ.
102#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
103pub struct CascadeDefects {
104    /// Atoms displaced.
105    pub displacements: u64,
106    /// Replacement events.
107    pub replacements: u64,
108    /// Vacancies left.
109    pub vacancies: u64,
110    /// Interstitials left.
111    pub interstitials: u64,
112}
113
114#[derive(Debug, Clone, Copy, PartialEq, Default)]
115struct DefectCounts {
116    displacements: u64,
117    replacements: u64,
118    interstitials: u64,
119}
120
121impl DefectCounts {
122    fn merge(&mut self, o: &DefectCounts) {
123        self.displacements += o.displacements;
124        self.replacements += o.replacements;
125        self.interstitials += o.interstitials;
126    }
127
128    fn report(&self) -> CascadeDefects {
129        CascadeDefects {
130            displacements: self.displacements,
131            replacements: self.replacements,
132            vacancies: self.displacements - self.replacements,
133            interstitials: self.interstitials,
134        }
135    }
136}
137
138/// Per-layer constants for the damage models (not merged).
139#[derive(Debug, Clone, PartialEq)]
140struct LayerConsts {
141    /// `(Z, E_d, E_b)` of each element, eV.
142    elems: Vec<(u8, f64, f64)>,
143    /// Atom-fraction mean atomic number.
144    mean_z: f64,
145    /// Mean atomic weight, u.
146    mean_a: f64,
147}
148
149#[derive(Debug, Clone, PartialEq)]
150struct LayerAcc {
151    stopped: Moments,
152    depth_hist: Histogram,
153    nrt: NrtDamage,
154    defects: DefectCounts,
155}
156
157#[derive(Debug, Clone, PartialEq)]
158struct FaceAcc {
159    count: u64,
160    energy_ev: f64,
161    energy_hist: Histogram,
162    polar_hist: Histogram,
163}
164
165impl FaceAcc {
166    fn new(c: &IonTallyConfig) -> Self {
167        Self {
168            count: 0,
169            energy_ev: 0.0,
170            energy_hist: Histogram::new(c.escape_energy),
171            polar_hist: Histogram::new(c.escape_polar),
172        }
173    }
174
175    fn merge(&mut self, o: &FaceAcc) {
176        self.count += o.count;
177        self.energy_ev += o.energy_ev;
178        self.energy_hist.merge(&o.energy_hist);
179        self.polar_hist.merge(&o.polar_hist);
180    }
181
182    fn report(&self, histories: u64) -> FaceEscape {
183        FaceEscape {
184            count: self.count,
185            per_ion: per(self.count as f64, histories),
186            energy_ev: self.energy_ev,
187            mean_energy_ev: if self.count > 0 {
188                self.energy_ev / self.count as f64
189            } else {
190                0.0
191            },
192            energy_histogram: self.energy_hist.clone(),
193            polar_histogram: self.polar_hist.clone(),
194        }
195    }
196}
197
198fn per(x: f64, histories: u64) -> f64 {
199    if histories == 0 {
200        0.0
201    } else {
202        x / histories as f64
203    }
204}
205
206/// The site of the most recent displacement, for replacement detection.
207#[derive(Debug, Clone, Copy, PartialEq)]
208struct Site {
209    pos: [f64; 3],
210    z: u8,
211    species: usize,
212    displacer_generation: u32,
213}
214
215/// The replacement rule of [`CascadeDefects`], shared with the per-element
216/// voxel tally ([`crate::tally::VoxelTally`]) so both count the same events.
217///
218/// Feed it the engine's hooks: [`Self::recoil`] on every `recoil`,
219/// [`Self::clear`] on `begin_history`, `escaped`, `end_history` and on a
220/// `lattice` deposit of kind [`LatticeDeposit::Subthreshold`], and
221/// [`Self::stopped`] on every `stopped`.
222#[derive(Debug, Clone, Copy, PartialEq, Default)]
223pub(crate) struct ReplacementTracker {
224    site: Option<Site>,
225}
226
227impl ReplacementTracker {
228    /// Forget the pending displacement site.
229    pub(crate) fn clear(&mut self) {
230        self.site = None;
231    }
232
233    /// A target atom was displaced: remember its site.
234    pub(crate) fn recoil(&mut self, r: &Particle) {
235        self.site = Some(Site {
236            pos: r.pos,
237            z: r.z,
238            species: r.species,
239            displacer_generation: r.generation.saturating_sub(1),
240        });
241    }
242
243    /// A particle stopped. If it fills the pending site, returns the species
244    /// of the displaced atom; either way the pending site is cleared.
245    pub(crate) fn stopped(&mut self, p: &Particle) -> Option<usize> {
246        let hit = match self.site {
247            Some(s) if s.pos == p.pos && s.z == p.z && s.displacer_generation == p.generation => {
248                Some(s.species)
249            }
250            _ => None,
251        };
252        self.site = None;
253        hit
254    }
255}
256
257/// Range, damage and escape tally for the BCA engine. See the module docs.
258#[derive(Debug, Clone, PartialEq)]
259pub struct IonTally {
260    config: IonTallyConfig,
261    species_z: Vec<u8>,
262    consts: Vec<LayerConsts>,
263    histories: u64,
264    budget: EnergyBudget,
265    depth: Moments,
266    lateral_y: Moments,
267    lateral_z: Moments,
268    radial: Moments,
269    depth_hist: Histogram,
270    lateral_y_hist: Histogram,
271    lateral_z_hist: Histogram,
272    radial_hist: Histogram,
273    layers: Vec<LayerAcc>,
274    displacement_hist: Histogram,
275    replacement_hist: Histogram,
276    interstitial_hist: Histogram,
277    /// NRT displacements by PKA depth.
278    nrt_hist: WeightedHistogram,
279    /// Kinchin-Pease displacements by PKA depth.
280    kp_hist: WeightedHistogram,
281    /// `[species] -> (front, back)`.
282    escapes: Vec<(FaceAcc, FaceAcc)>,
283    /// Per-history scratch: not part of the result, reset every history.
284    last_site: ReplacementTracker,
285}
286
287impl IonTally {
288    /// Empty tally for `stack`, with the particle species of the engine
289    /// (`species_z`, as returned by [`crate::ion::bca::Bca::species_z`]:
290    /// index 0 is the beam).
291    ///
292    /// Every element of every layer must have `E_d` and `E_b` set (the engine
293    /// requires the same).
294    pub fn new(
295        stack: &Stack,
296        species_z: &[u8],
297        config: IonTallyConfig,
298    ) -> Result<Self, MaterialError> {
299        let mut consts = Vec::with_capacity(stack.layers().len());
300        for l in stack.layers() {
301            let m = l.material();
302            let mut elems = Vec::new();
303            let mut mean_z = 0.0;
304            for (z, x) in m.atom_fractions() {
305                elems.push((
306                    z,
307                    m.displacement_energy_ev(z)?,
308                    m.lattice_binding_energy_ev(z)?,
309                ));
310                mean_z += x * f64::from(z);
311            }
312            consts.push(LayerConsts {
313                elems,
314                mean_z,
315                mean_a: m.mean_atomic_weight(),
316            });
317        }
318        let layers = (0..consts.len())
319            .map(|_| LayerAcc {
320                stopped: Moments::new(),
321                depth_hist: Histogram::new(config.depth),
322                nrt: NrtDamage::default(),
323                defects: DefectCounts::default(),
324            })
325            .collect();
326        Ok(Self {
327            config,
328            species_z: species_z.to_vec(),
329            consts,
330            histories: 0,
331            budget: EnergyBudget::default(),
332            depth: Moments::new(),
333            lateral_y: Moments::new(),
334            lateral_z: Moments::new(),
335            radial: Moments::new(),
336            depth_hist: Histogram::new(config.depth),
337            lateral_y_hist: Histogram::new(config.lateral),
338            lateral_z_hist: Histogram::new(config.lateral),
339            radial_hist: Histogram::new(config.radial),
340            layers,
341            displacement_hist: Histogram::new(config.depth),
342            replacement_hist: Histogram::new(config.depth),
343            interstitial_hist: Histogram::new(config.depth),
344            nrt_hist: WeightedHistogram::new(config.depth),
345            kp_hist: WeightedHistogram::new(config.depth),
346            escapes: species_z
347                .iter()
348                .map(|_| (FaceAcc::new(&config), FaceAcc::new(&config)))
349                .collect(),
350            last_site: ReplacementTracker::default(),
351        })
352    }
353
354    /// Histories tallied so far.
355    pub fn histories(&self) -> u64 {
356        self.histories
357    }
358
359    /// `E_d` and `E_b` of element `z` in layer `layer`.
360    fn energies(&self, layer: usize, z: u8) -> (f64, f64) {
361        self.consts[layer]
362            .elems
363            .iter()
364            .find(|e| e.0 == z)
365            .map(|e| (e.1, e.2))
366            .expect("recoils are elements of their layer")
367    }
368
369    /// The plain-data result. With `fit_dual_pearson` the depth histogram is
370    /// also fitted with [`DualPearson::fit`] (slower: a few thousand
371    /// evaluations over the histogram).
372    pub fn report(&self, fit_dual_pearson: bool) -> IonReport {
373        let h = self.histories;
374        let depth = self.depth.summary();
375        let (pearson_iv, pearson_iv_error) = match depth {
376            None => (None, None),
377            Some(s) => match PearsonIv::from_moments(s.mean, s.std_dev, s.skewness, s.kurtosis) {
378                Ok(p) => (Some(p), None),
379                Err(e) => (None, Some(e)),
380            },
381        };
382        let (dual_pearson, dual_pearson_error) = if fit_dual_pearson {
383            match DualPearson::fit(&self.depth_hist) {
384                Ok(f) => (Some(f), None),
385                Err(e) => (None, Some(e)),
386            }
387        } else {
388            (None, None)
389        };
390        let range = RangeReport {
391            stopped: self.depth.n,
392            depth,
393            pearson_iv,
394            pearson_iv_error,
395            dual_pearson,
396            dual_pearson_error,
397            depth_histogram: self.depth_hist.clone(),
398            lateral_y: self.lateral_y.summary(),
399            lateral_z: self.lateral_z.summary(),
400            radial: self.radial.summary(),
401            lateral_y_histogram: self.lateral_y_hist.clone(),
402            lateral_z_histogram: self.lateral_z_hist.clone(),
403            radial_histogram: self.radial_hist.clone(),
404            layers: self
405                .layers
406                .iter()
407                .enumerate()
408                .map(|(i, l)| LayerRange {
409                    layer: i,
410                    stopped: l.stopped.n,
411                    depth: l.stopped.summary(),
412                    depth_histogram: l.depth_hist.clone(),
413                })
414                .collect(),
415        };
416        let mut nrt = NrtDamage::default();
417        let mut defects = DefectCounts::default();
418        for l in &self.layers {
419            nrt.merge(&l.nrt);
420            defects.merge(&l.defects);
421        }
422        let mut vacancy_hist = self.displacement_hist.clone();
423        for (v, r) in vacancy_hist
424            .counts
425            .iter_mut()
426            .zip(&self.replacement_hist.counts)
427        {
428            *v -= r;
429        }
430        vacancy_hist.underflow -= self.replacement_hist.underflow;
431        vacancy_hist.overflow -= self.replacement_hist.overflow;
432        let damage = DamageReport {
433            nrt,
434            cascade: defects.report(),
435            layers: self
436                .layers
437                .iter()
438                .enumerate()
439                .map(|(i, l)| LayerDamage {
440                    layer: i,
441                    nrt: l.nrt,
442                    cascade: l.defects.report(),
443                })
444                .collect(),
445            vacancy_histogram: vacancy_hist,
446            interstitial_histogram: self.interstitial_hist.clone(),
447            replacement_histogram: self.replacement_hist.clone(),
448            nrt_histogram: self.nrt_hist.clone(),
449            kinchin_pease_histogram: self.kp_hist.clone(),
450        };
451        let species: Vec<SpeciesEscape> = self
452            .escapes
453            .iter()
454            .enumerate()
455            .map(|(i, (f, b))| SpeciesEscape {
456                species: i,
457                z: self.species_z[i],
458                beam: i == 0,
459                front: f.report(h),
460                back: b.report(h),
461            })
462            .collect();
463        let sputtered: u64 = self.escapes.iter().skip(1).map(|e| e.0.count).sum();
464        let escapes = EscapeReport {
465            backscatter_coefficient: per(self.escapes[0].0.count as f64, h),
466            transmission_coefficient: per(self.escapes[0].1.count as f64, h),
467            energy_reflection_coefficient: if self.budget.incident > 0.0 {
468                self.budget.backscattered / self.budget.incident
469            } else {
470                0.0
471            },
472            sputter_yield: per(sputtered as f64, h),
473            species,
474        };
475        IonReport {
476            histories: h,
477            species_z: self.species_z.clone(),
478            budget: self.budget,
479            range,
480            damage,
481            escapes,
482        }
483    }
484}
485
486impl BcaTally for IonTally {
487    fn begin_history(&mut self, _index: u64) {
488        self.last_site.clear();
489    }
490
491    fn lattice(&mut self, _at: [f64; 3], _layer: usize, kind: LatticeDeposit, _energy_ev: f64) {
492        if kind == LatticeDeposit::Subthreshold {
493            self.last_site.clear();
494        }
495    }
496
497    fn recoil(&mut self, r: &Particle) {
498        let layer = r.layer;
499        let (e_d, e_b) = self.energies(layer, r.z);
500        let acc = &mut self.layers[layer];
501        acc.defects.displacements += 1;
502        self.displacement_hist.fill(r.pos[0]);
503        if r.generation == 1 {
504            // A PKA: the transfer T is the recoil energy plus the binding
505            // energy the engine left in the lattice.
506            let c = &self.consts[layer];
507            let t = r.energy_ev + e_b;
508            let t_dam = damage_energy_ev(t, f64::from(r.z), r.mass_amu, c.mean_z, c.mean_a);
509            acc.nrt.pka_count += 1;
510            acc.nrt.pka_energy_ev += t;
511            acc.nrt.damage_energy_ev += t_dam;
512            let n_nrt = nrt_displacements(t_dam, e_d);
513            let n_kp = kinchin_pease(t_dam, e_d);
514            acc.nrt.nrt_displacements += n_nrt;
515            acc.nrt.kinchin_pease_displacements += n_kp;
516            self.nrt_hist.fill(r.pos[0], n_nrt);
517            self.kp_hist.fill(r.pos[0], n_kp);
518        }
519        self.last_site.recoil(r);
520    }
521
522    fn stopped(&mut self, p: &Particle) {
523        let replaced = self.last_site.stopped(p).is_some();
524        let x = p.pos[0];
525        if replaced {
526            self.layers[p.layer].defects.replacements += 1;
527            self.replacement_hist.fill(x);
528        } else if !p.is_primary() {
529            self.layers[p.layer].defects.interstitials += 1;
530            self.interstitial_hist.fill(x);
531        }
532        if p.is_primary() {
533            let (y, z) = (p.pos[1], p.pos[2]);
534            let r = y.hypot(z);
535            self.depth.push(x);
536            self.lateral_y.push(y);
537            self.lateral_z.push(z);
538            self.radial.push(r);
539            self.depth_hist.fill(x);
540            self.lateral_y_hist.fill(y);
541            self.lateral_z_hist.fill(z);
542            self.radial_hist.fill(r);
543            let l = &mut self.layers[p.layer];
544            l.stopped.push(x);
545            l.depth_hist.fill(x);
546        }
547    }
548
549    fn escaped(&mut self, p: &Particle, face: Face) {
550        self.last_site.clear();
551        let polar = p.dir[0].abs().min(1.0).acos();
552        let pair = &mut self.escapes[p.species];
553        let acc = match face {
554            Face::Front => &mut pair.0,
555            Face::Back => &mut pair.1,
556            // Lateral faces exist only on voxel grids, which this stack
557            // tally does not describe.
558            Face::Side => return,
559        };
560        acc.count += 1;
561        acc.energy_ev += p.energy_ev;
562        acc.energy_hist.fill(p.energy_ev);
563        acc.polar_hist.fill(polar);
564    }
565
566    fn end_history(&mut self, _index: u64, budget: &EnergyBudget) {
567        self.histories += 1;
568        self.budget.add(budget);
569        self.last_site.clear();
570    }
571
572    /// Field-wise merge.
573    ///
574    /// # Panics
575    /// If `other` was built with a different stack, species list or
576    /// configuration.
577    fn merge(&mut self, o: Self) {
578        assert!(
579            self.config == o.config && self.species_z == o.species_z && self.consts == o.consts,
580            "merging ion tallies built for different runs"
581        );
582        self.histories += o.histories;
583        self.budget.add(&o.budget);
584        self.depth.merge(&o.depth);
585        self.lateral_y.merge(&o.lateral_y);
586        self.lateral_z.merge(&o.lateral_z);
587        self.radial.merge(&o.radial);
588        self.depth_hist.merge(&o.depth_hist);
589        self.lateral_y_hist.merge(&o.lateral_y_hist);
590        self.lateral_z_hist.merge(&o.lateral_z_hist);
591        self.radial_hist.merge(&o.radial_hist);
592        for (a, b) in self.layers.iter_mut().zip(&o.layers) {
593            a.stopped.merge(&b.stopped);
594            a.depth_hist.merge(&b.depth_hist);
595            a.nrt.merge(&b.nrt);
596            a.defects.merge(&b.defects);
597        }
598        self.displacement_hist.merge(&o.displacement_hist);
599        self.replacement_hist.merge(&o.replacement_hist);
600        self.interstitial_hist.merge(&o.interstitial_hist);
601        self.nrt_hist.merge(&o.nrt_hist);
602        self.kp_hist.merge(&o.kp_hist);
603        for (a, b) in self.escapes.iter_mut().zip(&o.escapes) {
604            a.0.merge(&b.0);
605            a.1.merge(&b.1);
606        }
607    }
608}
609
610/// Everything an [`IonTally`] measured, as plain data.
611#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
612pub struct IonReport {
613    /// Primary histories.
614    pub histories: u64,
615    /// Atomic number by species index (0 is the beam).
616    pub species_z: Vec<u8>,
617    /// Summed energy budget of all histories, eV.
618    pub budget: EnergyBudget,
619    /// Range distributions of the stopped beam particles.
620    pub range: RangeReport,
621    /// Damage: model estimates and cascade counts, side by side.
622    pub damage: DamageReport,
623    /// Backscattering, transmission and sputtering.
624    pub escapes: EscapeReport,
625}
626
627/// Range distributions of beam particles that came to rest in the target.
628#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
629pub struct RangeReport {
630    /// Beam particles at rest in the target.
631    pub stopped: u64,
632    /// Depth moments (`Rp`, `ΔRp`, `γ`, `β`), m. `None` with fewer than two
633    /// stopped particles.
634    pub depth: Option<MomentSummary>,
635    /// Pearson IV with the depth moments, if they are in the type IV region.
636    pub pearson_iv: Option<PearsonIv>,
637    /// Why `pearson_iv` is `None`, when there were moments to try.
638    pub pearson_iv_error: Option<PearsonError>,
639    /// Dual-Pearson fit to the depth histogram, if requested and possible.
640    pub dual_pearson: Option<DualPearsonFit>,
641    /// Why the requested dual-Pearson fit failed.
642    pub dual_pearson_error: Option<PearsonError>,
643    /// Depth histogram, m.
644    pub depth_histogram: Histogram,
645    /// Moments of the lateral coordinate `y`, m.
646    pub lateral_y: Option<MomentSummary>,
647    /// Moments of the lateral coordinate `z`, m.
648    pub lateral_z: Option<MomentSummary>,
649    /// Moments of the radial distance from the beam axis through the entry
650    /// point, m.
651    pub radial: Option<MomentSummary>,
652    /// Histogram of `y`, m.
653    pub lateral_y_histogram: Histogram,
654    /// Histogram of `z`, m.
655    pub lateral_z_histogram: Histogram,
656    /// Histogram of the radial distance, m.
657    pub radial_histogram: Histogram,
658    /// The same depth statistics per layer (by where the particle stopped).
659    pub layers: Vec<LayerRange>,
660}
661
662/// Depth statistics of the beam particles that stopped in one layer.
663#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
664pub struct LayerRange {
665    /// Layer index (0 is the front layer).
666    pub layer: usize,
667    /// Beam particles at rest in this layer.
668    pub stopped: u64,
669    /// Their depth moments, m.
670    pub depth: Option<MomentSummary>,
671    /// Their depth histogram, m (overall binning).
672    pub depth_histogram: Histogram,
673}
674
675/// Damage, as two separate kinds of result. See [`crate::ion::damage`].
676#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
677pub struct DamageReport {
678    /// NRT / Kinchin-Pease model estimates, all layers.
679    pub nrt: NrtDamage,
680    /// Simulated cascade defect counts, all layers.
681    pub cascade: CascadeDefects,
682    /// Both, per layer (by the layer of the event).
683    pub layers: Vec<LayerDamage>,
684    /// Depth histogram of vacancies, m.
685    pub vacancy_histogram: Histogram,
686    /// Depth histogram of interstitials, m.
687    pub interstitial_histogram: Histogram,
688    /// Depth histogram of replacements, m.
689    pub replacement_histogram: Histogram,
690    /// NRT displacement estimate by depth of the PKA, m. A model estimate
691    /// (damage-energy based, no cascade spread), not an event count; its sum
692    /// over bins, underflow and overflow is `nrt.nrt_displacements`.
693    pub nrt_histogram: WeightedHistogram,
694    /// Kinchin-Pease displacement estimate by PKA depth, m; sums to
695    /// `nrt.kinchin_pease_displacements`.
696    pub kinchin_pease_histogram: WeightedHistogram,
697}
698
699/// Damage in one layer.
700#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
701pub struct LayerDamage {
702    /// Layer index.
703    pub layer: usize,
704    /// Model estimates for PKAs created in this layer.
705    pub nrt: NrtDamage,
706    /// Cascade counts for events in this layer.
707    pub cascade: CascadeDefects,
708}
709
710/// Particles leaving the target.
711#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
712pub struct EscapeReport {
713    /// Beam particles leaving the front face, per incident ion.
714    pub backscatter_coefficient: f64,
715    /// Beam particles leaving the back face, per incident ion.
716    pub transmission_coefficient: f64,
717    /// Energy carried out of the front face by beam particles, as a fraction
718    /// of the incident energy.
719    pub energy_reflection_coefficient: f64,
720    /// Target atoms leaving the front face, all species, per incident ion.
721    pub sputter_yield: f64,
722    /// Per species (index 0 is the beam).
723    pub species: Vec<SpeciesEscape>,
724}
725
726/// Escapes of one species through each face.
727#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
728pub struct SpeciesEscape {
729    /// Species index.
730    pub species: usize,
731    /// Atomic number.
732    pub z: u8,
733    /// True for the beam species (front escapes are backscattering), false
734    /// for a target element (front escapes are sputtering).
735    pub beam: bool,
736    /// Through the front face (`x = 0`).
737    pub front: FaceEscape,
738    /// Through the back face of a finite stack.
739    pub back: FaceEscape,
740}
741
742/// Count and spectra of particles leaving through one face.
743#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
744pub struct FaceEscape {
745    /// Particles.
746    pub count: u64,
747    /// Particles per incident ion (backscatter or transmission coefficient,
748    /// or partial sputter yield).
749    pub per_ion: f64,
750    /// Sum of their energies outside the target, eV.
751    pub energy_ev: f64,
752    /// Mean energy outside the target, eV (0 if none escaped).
753    pub mean_energy_ev: f64,
754    /// Energy spectrum, eV.
755    pub energy_histogram: Histogram,
756    /// Polar angle from the outward normal, rad.
757    pub polar_histogram: Histogram,
758}