lindhard/tally/

electron.rs

//! The electron tally: energy deposition, emission yields and spectra, wired
//! to the electron transport hooks.
//!
//! [`FullElectronTally`] implements [`ElectronTally`]. Build one prototype
//! with [`FullElectronTally::new`] and pass clones to
//! [`Transport::run`](crate::electron::transport::Transport::run). The
//! merged tally turns into a plain-data, serializable [`ElectronReport`] with
//! [`FullElectronTally::report`].
//!
//! Units: positions and lengths in metres, energies in eV, angles in radians.
//! Totals are summed over all histories. Divide by [`ElectronReport::histories`]
//! for per-primary values (the yields are already per primary).
//!
//! # What is tallied
//!
//! - **Energy deposition.** An inelastic collision deposits the energy it
//!   leaves in the solid at the collision point: its whole loss `W` without a
//!   secondary model, and with one only the part no secondary carries away
//!   (see [the energy balance](#the-energy-balance)). An electron whose
//!   energy falls below the stopping threshold of its layer deposits what it
//!   has left at the point where it stopped. Both go into a per-layer total
//!   and into optional 3D grids: a
//!   Cartesian `x`-`y`-`z` grid ([`CartesianGrid`]) and a cylindrical `r`-`x`
//!   grid about the beam axis ([`CylindricalGrid`]). Energy deposited outside
//!   a grid is kept in that grid's `outside_ev`, so every grid sums to the
//!   total deposited energy. Placing the loss at the collision point and the
//!   sub-threshold remainder at the stopping point is this crate's bookkeeping
//!   convention for a point-collision transport loop, not a physics model.
//!   The same deposits optionally feed a radial profile in one depth slab,
//!   on log radial bins and with per-history errors
//!   ([`ElectronTallyConfig::psf`], [`super::psf`]), which is the input of
//!   the point-spread-function fits.
//! - **Ionisation and lattice channels.** Every deposit is also booked to
//!   one of two channels, per layer and per cell of both grids, see
//!   [the channel split](#the-ionisation-and-lattice-channels).
//! - **Emission.** Every electron the `escaped` hook reports is counted at its
//!   face, with its energy and its polar angle from the outward surface
//!   normal. It is classed as *slow* (`E < split`, the secondary-electron
//!   class) or *fast* (`E >= split`, the backscattered class), see
//!   [the split convention](#the-sebse-split). The front-face counts per
//!   primary are the backscatter yield `η` (fast) and the secondary yield `δ`
//!   (slow). The energy spectrum and the polar-angle spectrum (one per class)
//!   are histograms per face.
//! - **Generation-volume moments.** Energy-weighted mean and standard
//!   deviation of the deposition position along `x`, `y`, `z` and of the
//!   radial distance from the beam axis (the line `y = z = 0` through the
//!   entry point), plus the energy-weighted RMS radius. See
//!   [`GenerationVolume`].
//! - **Stopping points.** Unweighted moments of the depth and radial distance
//!   of electrons that fell below the stopping threshold (an electron range
//!   distribution), and the same moments for the primaries alone
//!   ([`PrimaryStoppingPoints`]: the penetration depth of stopped primaries).
//! - **Energy balance.** The incident energy plus the Fermi-sea source equals
//!   the deposited, escaped, trapped and barrier terms, see
//!   [the energy balance](#the-energy-balance) and [`ElectronEnergyBudget`].
//! - **Table coverage.** A numerical diagnostic, not a physical result: per
//!   layer and channel, how many rate evaluations of the transport read the
//!   cross-section table inside its energy grid and how many used the
//!   constant continuation below or above it
//!   ([`super::table_coverage`], [`ElectronReport::table_coverage`]). These
//!   are rate evaluations, not collisions.
//! - **Event-cap counts.** Another numerical diagnostic: how many secondary
//!   tracks the collision cap cut off, and how many primary histories had
//!   any electron (the primary or a secondary) cut off by it
//!   ([`EventCapCounts`], [`ElectronReport::event_caps`]). The primary's own
//!   cap stays in [`FateCounts::event_capped`]; zero there does not mean the
//!   secondary cascades completed.
//!
//! # The SE/BSE split
//!
//! An electron leaving the target with energy below 50 eV is by convention a
//! secondary electron and one at higher energy a backscattered electron:
//! H. Chen, Y. Zou, S. Mao, M. S. S. Khan, K. Tőkési and Z. J. Ding,
//! "Influence of energy loss function to the Monte Carlo simulated electron
//! backscattering coefficient", Sci. Rep. 12, 18201 (2022),
//! doi:10.1038/s41598-022-20466-3 (open access), Introduction ("These
//! electrons, with their energies higher than 50 eV, are defined as
//! backscattered electrons; while the lower energy electrons are the excited
//! secondary electrons", citing H. Niedrig, Scanning 1, 17 (1978)) and section
//! "Monte Carlo simulation" ("an escaped electron is counted either as a true
//! secondary electron (< 50 eV) or a backscattered electron (> 50 eV)"). The
//! class depends only on the energy, not on whether the electron is a primary
//! or a secondary.
//!
//! That paper leaves `E = 50 eV` itself unassigned. This tally puts it in the
//! fast (backscattered) class, as Nebula's detector flags do: `DETECTOR_LT50`
//! ("if energy < 50 eV") and `DETECTOR_GE50` ("if energy >= 50 eV") in
//! `source/core/cpu_material_manager.h`, Nebula commit `a50a8e8`
//! (<https://github.com/Nebula-simulator/Nebula>, BSD-3-Clause). Only the
//! convention is taken from there; no code is ported.
//!
//! The split is configurable ([`ElectronTallyConfig::se_bse_split_ev`],
//! default [`SE_BSE_SPLIT_EV`]) and is recorded, with its source, in
//! [`ElectronTallyMetadata`].
//!
//! # Event kinds and how later channels slot in
//!
//! The tally reads only the hooks of [`ElectronTally`] and does not depend on
//! how many electrons a history contains:
//!
//! - Emission is counted per `escaped` call, so secondaries that escape are
//!   counted with the primaries, as the convention requires.
//! - `stopped` is classified by energy against the stopping threshold of the
//!   electron's layer ([`Transport::stopping_thresholds_ev`]: the cutoff, or
//!   `U + cutoff` under
//!   [`CutoffReference::VacuumLevel`](crate::electron::transport::CutoffReference::VacuumLevel)):
//!   below it the electron stopped and its energy is a deposit, at or above
//!   it the electron is trapped (the loop calls `stopped` for
//!   [`Fate::Trapped`] too, and a trapped electron always has at least the
//!   threshold energy, because the loop checks the threshold after every
//!   energy change, a primary starts above it and a secondary below it is
//!   not created).
//! - `absorbed` (the back face under
//!   [`EscapeRule::FrontOnly`](crate::electron::transport::EscapeRule::FrontOnly))
//!   counts as trapped energy.
//! - [`Fate::EventCap`] has no hook of its own, so the energy of the last
//!   state seen for the capped electron is counted as trapped: at
//!   `end_secondary` for a secondary, at `end_history` for the primary. The
//!   last state is kept per electron; the primary's is set aside when the
//!   first secondary begins, so a secondary cannot overwrite it. Each capped
//!   secondary is also counted once in [`EventCapCounts::secondary_tracks`],
//!   and a history with any capped electron once in
//!   [`EventCapCounts::affected_histories`]; these counts do not touch the
//!   energy accounting.
//!
//! # The energy balance
//!
//! Inside a layer, kinetic energies are measured from the band bottom, and
//! with [`BoundaryModel::StepBarrier`](crate::electron::transport::BoundaryModel::StepBarrier)
//! an electron gains or loses the step `ΔU = U' - U` in inner potential at
//! every face it gets through ([`crate::electron::boundary`]; `U = E_F + Φ`
//! in a metal, Verduin's thesis Eq. 3.136, cited there). The primary's
//! incident energy is its energy where it starts, its vacuum energy with the
//! step barrier. Summed over all histories:
//!
//! ```text
//! incident + fermi_sea = deposited + escaped + trapped + barrier
//! ```
//!
//! - `barrier` is the sum of `-ΔU` over every face transmission (the
//!   `barrier` hook): kinetic energy taken by the potential steps. It is
//!   `-U` for a primary that enters (its vacuum energy becomes `E + U`
//!   inside), `+U` for any electron that leaves through a surface, and
//!   `-(U' - U)` at an interface. A primary reflected at entry changes
//!   nothing and is counted as escaped with its incident energy.
//! - With a secondary model, an inelastic event takes `W` from the electron
//!   and liberates, if anything, an electron of kinetic energy
//!   `E_SE = E_F + W - B` (Verduin Eq. 3.86, see
//!   [`crate::electron::secondary`]), which the tally follows like any other
//!   electron until it stops, escapes or is trapped. Per event
//!   ([`SecondaryEvent`]) `W = secondary + binding + deposited`, with
//!   `binding = B - E_F`. Only `deposited`, plus `binding` where it is
//!   positive (an initial state below the band bottom), stays in the solid
//!   at the event and is deposited there. Where `binding` is negative (a
//!   conduction electron, as in a metal, which already had `E_F - B` above
//!   the band bottom), the liberated electron brings `-binding` of its own:
//!   that is the `fermi_sea` source term. So `W + fermi_sea = secondary +
//!   deposit` event by event, and the secondary's energy is counted once,
//!   where it ends.
//!
//! - With the opt-in insulator channels of a layer
//!   ([`crate::electron::phonon`]), an LO-phonon emission gives `ħω` to the
//!   lattice (deposited where it happens), an absorption takes `ħω` from it
//!   (the `phonon_absorbed` source, added to the left side of the balance),
//!   and an electron trapped as a polaron deposits its whole remaining energy
//!   where it is trapped (the `polaron` deposit). The balance is then
//!   `incident + fermi_sea + phonon_absorbed = deposited + escaped + trapped +
//!   barrier`, with `deposited` including `phonon_emitted` and `polaron`.
//!   Polaron traps are not counted in the stopping points, which are
//!   electrons that fell below the stopping threshold.
//!
//! Without a secondary model and without the step barrier, `fermi_sea` and
//! `barrier` are zero and the balance is `incident = deposited + escaped +
//! trapped`.
//!
//! ## Deposits and the energy imparted
//!
//! Deposits (`deposited`, the per-layer totals and every grid cell) are
//! measured from the band bottom. With a secondary model in a layer with a
//! Fermi energy, a liberated conduction electron's kinetic energy includes
//! the `-binding` it already had (the `fermi_sea` source), and wherever that
//! electron stops, its whole remaining energy above the band bottom is
//! deposited. So deposits include Fermi-sea energy the beam never supplied,
//! and `deposited` can exceed the energy the beam imparted to the target.
//! Moving the sources of the balance above to the right side gives the net
//! energy imparted:
//!
//! ```text
//! incident - escaped = deposited + trapped + barrier - fermi_sea - phonon_absorbed
//! ```
//!
//! Without a secondary model and without the insulator channels the two
//! source terms are zero.
//!
//! # The ionisation and lattice channels
//!
//! The kind of a deposit is known at the call site and carried with it, not
//! reconstructed from totals. Two channels, with `ionisation + lattice =
//! deposited` per layer, per grid cell and in the budget (to floating-point
//! rounding; each is its own accumulator):
//!
//! - **Ionisation** (electronic excitation, the source of electron-hole
//!   pairs): the part of an inelastic event left in the solid (`W` without a
//!   secondary model; `deposited + max(binding, 0)` with one, from the
//!   `secondary` hook), and the **residual** energy of electrons that fell
//!   below the stopping threshold ([`ElectronEnergyBudget::residual_ev`]).
//! - **Lattice**: LO-phonon emission (`phonon_emitted`) and the whole
//!   remaining energy of an electron trapped as a polaron (`polaron`).
//!
//! The residual is the cutoff-dependent term. It is placed in the ionisation
//! channel by convention: a sub-cutoff electron is below the range where this
//! engine follows collisions, and its energy goes into further electronic
//! excitation as it thermalises. This is a bookkeeping convention, not a
//! computed split; the lattice channel is the lower bound of heating, and the
//! ionisation channel includes the residual and, like all deposits, Fermi-sea
//! energy (see [the energy imparted](#deposits-and-the-energy-imparted)).
//! Elastic nuclear recoil energy is not tallied by the electron engine and is
//! in neither channel.
//!
//! # Determinism
//!
//! Counts are integers, histograms have integer bins, and the grids, totals
//! and weighted moments are `f64` accumulators merged field by field in the
//! order [`crate::rng::run_particles`] calls `merge` (chunk order). For a
//! fixed `chunk_size` the report is bit-identical at any thread count
//! (`lindhard/tests/electron_tally.rs`). Nothing is stored in a hash map.

use std::collections::BTreeMap;

use serde::{Deserialize, Serialize};

use super::dose::{DoseAcc, DoseConfig, DoseError, DoseReport};
use super::exposure::{ExposureAccumulator, ExposureConfig, ExposureError, ExposureReport};
use super::hist::{Binning, Histogram};
use super::moments::{MomentSummary, Moments};
use super::psf::{PsfConfig, PsfError, RadialAccumulator, RadialProfile};
use super::table_coverage::{LayerTableCoverage, TableCoverageTally};
use super::voxel::{VoxelCell, VoxelRecord, VoxelReport, VoxelTallyGrid};
use crate::electron::secondary::{SecondaryEvent, SecondaryModel};
use crate::electron::transport::{
    Boundary, ElectronState, ElectronTally, Face, Fate, GridCoverage, PhononEvent, TableChannel,
    Transport,
};
use crate::geometry::Geometry;

/// Default SE/BSE energy split, eV: an escaping electron below this energy is
/// in the secondary (slow) class, at or above it in the backscattered (fast)
/// class. Chen et al., Sci. Rep. 12, 18201 (2022),
/// doi:10.1038/s41598-022-20466-3, Introduction and section "Monte Carlo
/// simulation"; the assignment of the boundary itself follows Nebula's
/// `DETECTOR_GE50` flag. See [the module docs](self#the-sebse-split).
pub const SE_BSE_SPLIT_EV: f64 = 50.0;

/// The source of the default split, as recorded in [`ElectronTallyMetadata`].
pub const SE_BSE_SPLIT_SOURCE: &str = "Chen, Zou, Mao, Khan, Tokesi and Ding, Sci. Rep. 12, 18201 \
     (2022), doi:10.1038/s41598-022-20466-3 (Introduction; section Monte Carlo simulation), citing \
     Niedrig, Scanning 1, 17 (1978); E = split is counted as backscattered, as Nebula's \
     DETECTOR_GE50 flag (source/core/cpu_material_manager.h, commit a50a8e8) does";

/// The classification rule, as recorded in [`ElectronTallyMetadata`].
pub const SE_BSE_SPLIT_RULE: &str =
    "slow (secondary) if E < split; fast (backscattered) if E >= split; E as reported by the escaped hook";

/// Errors from building a [`FullElectronTally`].
#[derive(Debug, Clone, PartialEq, thiserror::Error)]
pub enum ElectronTallyError {
    /// The SE/BSE split is not finite and non-negative.
    #[error("invalid SE/BSE split {0} eV: must be finite and non-negative")]
    Split(f64),
    /// The radial binning of a cylindrical grid starts below `r = 0`.
    #[error("the cylindrical grid's radial bins must start at r >= 0, got lo = {0} m")]
    NegativeRadius(f64),
    /// The escape polar-angle binning is not within `[0, π]`.
    #[error("the escape polar-angle bins must lie in [0, pi], got [{lo}, {hi})")]
    PolarRange {
        /// Lower edge.
        lo: f64,
        /// Upper edge.
        hi: f64,
    },
    /// A grid has more cells than fit in memory addressing.
    #[error("the {0} grid has too many cells")]
    GridTooLarge(&'static str),
    /// The radial-profile (PSF) configuration is invalid.
    #[error("invalid PSF configuration: {0}")]
    Psf(#[from] PsfError),
    /// The dose / depth-dose / transmission configuration is invalid.
    #[error(transparent)]
    Dose(#[from] DoseError),
    /// The exposure settings are invalid.
    #[error(transparent)]
    Exposure(#[from] ExposureError),
}

/// A Cartesian deposition grid over `x` (depth), `y` and `z`, m.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct CartesianGrid {
    /// Depth bins, m.
    pub x: Binning,
    /// Lateral `y` bins, m.
    pub y: Binning,
    /// Lateral `z` bins, m.
    pub z: Binning,
}

impl CartesianGrid {
    /// Number of voxels, or `None` on overflow.
    pub fn voxels(&self) -> Option<usize> {
        self.x
            .bins
            .checked_mul(self.y.bins)?
            .checked_mul(self.z.bins)
    }

    /// Flat index of voxel `(ix, iy, iz)`: `(ix * ny + iy) * nz + iz`.
    pub fn index(&self, ix: usize, iy: usize, iz: usize) -> usize {
        (ix * self.y.bins + iy) * self.z.bins + iz
    }

    /// Volume of one voxel, m³.
    pub fn voxel_volume_m3(&self) -> f64 {
        self.x.width() * self.y.width() * self.z.width()
    }

    fn locate(&self, pos: [f64; 3]) -> Option<usize> {
        let ix = self.x.locate(pos[0]).ok()?;
        let iy = self.y.locate(pos[1]).ok()?;
        let iz = self.z.locate(pos[2]).ok()?;
        Some(self.index(ix, iy, iz))
    }
}

/// A cylindrical deposition grid about the beam axis (the line `y = z = 0`
/// through the entry point, along `x`): radial distance `r = (y² + z²)^(1/2)`
/// and depth `x`, m.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct CylindricalGrid {
    /// Radial bins, m; `r.lo` must be at least 0.
    pub r: Binning,
    /// Depth bins, m.
    pub depth: Binning,
}

impl CylindricalGrid {
    /// Number of cells, or `None` on overflow.
    pub fn cells(&self) -> Option<usize> {
        self.r.bins.checked_mul(self.depth.bins)
    }

    /// Flat index of cell `(ir, ix)`: `ir * n_depth + ix`.
    pub fn index(&self, ir: usize, ix: usize) -> usize {
        ir * self.depth.bins + ix
    }

    /// Volume of a cell in radial bin `ir` (an annulus times a depth bin), m³.
    pub fn cell_volume_m3(&self, ir: usize) -> f64 {
        let (a, b) = (self.r.edge(ir), self.r.edge(ir + 1));
        std::f64::consts::PI * (b * b - a * a) * self.depth.width()
    }

    fn locate(&self, pos: [f64; 3]) -> Option<usize> {
        let ir = self.r.locate(pos[1].hypot(pos[2])).ok()?;
        let ix = self.depth.locate(pos[0]).ok()?;
        Some(self.index(ir, ix))
    }
}

/// Settings of a [`FullElectronTally`].
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct ElectronTallyConfig {
    /// Optional Cartesian deposition grid.
    pub cartesian: Option<CartesianGrid>,
    /// Optional cylindrical `r`-`x` deposition grid.
    pub cylindrical: Option<CylindricalGrid>,
    /// Energy of escaping electrons, eV.
    pub escape_energy: Binning,
    /// Polar angle of escaping electrons from the outward surface normal, rad
    /// (`0` along the normal, `π/2` grazing). Must lie within `[0, π]`.
    pub escape_polar: Binning,
    /// SE/BSE split, eV: escaping electrons below it are slow (secondary),
    /// at or above it fast (backscattered). See
    /// [the module docs](self#the-sebse-split).
    pub se_bse_split_ev: f64,
    /// Optional radial profile in a depth slab, log radial bins, with
    /// per-history errors: the input of the PSF fits of [`super::psf`].
    #[serde(default)]
    pub psf: Option<PsfConfig>,
    /// Optional dose per layer, depth-dose profile and interface
    /// transmission ([`super::dose`]).
    #[serde(default)]
    pub dose: Option<DoseConfig>,
    /// Optional exposure statistics: the per-electron deposit in target
    /// volumes and the compound-Poisson distribution of a finite exposure
    /// ([`super::exposure`]).
    #[serde(default)]
    pub exposure: Option<ExposureConfig>,
}

impl ElectronTallyConfig {
    /// Spectra with the given binnings, no deposition grids, no PSF profile,
    /// and the [`SE_BSE_SPLIT_EV`] split.
    pub fn new(escape_energy: Binning, escape_polar: Binning) -> Self {
        Self {
            cartesian: None,
            cylindrical: None,
            escape_energy,
            escape_polar,
            se_bse_split_ev: SE_BSE_SPLIT_EV,
            psf: None,
            dose: None,
            exposure: None,
        }
    }
}

/// Energy-weighted mean and centred second moment of one coordinate,
/// mergeable. The pairwise update is that of T. F. Chan, G. H. Golub and
/// R. J. LeVeque, Stanford report STAN-CS-79-773 (1979) (the update also used
/// by [`super::moments`]), with the counts replaced by weights, which is the
/// weighted form of D. H. D. West, "Updating mean and variance estimates: an
/// improved method", Commun. ACM 22, 532 (1979), doi:10.1145/359146.359153.
/// With `W = W_a + W_b` and `δ = mean_b - mean_a`:
/// `mean = mean_a + δ W_b / W`, `M2 = M2_a + M2_b + δ² W_a W_b / W`.
#[derive(Debug, Clone, Copy, PartialEq, Default)]
struct WeightedMoments {
    weight: f64,
    mean: f64,
    m2: f64,
}

impl WeightedMoments {
    fn push(&mut self, x: f64, w: f64) {
        self.merge_parts(w, x, 0.0);
    }

    fn merge(&mut self, o: &WeightedMoments) {
        self.merge_parts(o.weight, o.mean, o.m2);
    }

    fn merge_parts(&mut self, w_b: f64, mean_b: f64, m2_b: f64) {
        if w_b <= 0.0 {
            return;
        }
        if self.weight <= 0.0 {
            *self = Self {
                weight: w_b,
                mean: mean_b,
                m2: m2_b,
            };
            return;
        }
        let w = self.weight + w_b;
        let d = mean_b - self.mean;
        self.m2 += m2_b + d * d * self.weight * w_b / w;
        self.mean += d * w_b / w;
        self.weight = w;
    }

    fn variance(&self) -> f64 {
        if self.weight > 0.0 {
            (self.m2 / self.weight).max(0.0)
        } else {
            0.0
        }
    }
}

/// Energy sums, eV.
#[derive(Debug, Clone, Copy, PartialEq, Default)]
struct Budget {
    incident: f64,
    inelastic: f64,
    residual: f64,
    escaped_front: f64,
    escaped_back: f64,
    escaped_side: f64,
    no_interaction: f64,
    absorbed: f64,
    event_cap: f64,
    fermi_sea: f64,
    barrier: f64,
    phonon_emitted: f64,
    phonon_absorbed: f64,
    polaron: f64,
}

impl Budget {
    fn merge(&mut self, o: &Budget) {
        self.incident += o.incident;
        self.inelastic += o.inelastic;
        self.residual += o.residual;
        self.escaped_front += o.escaped_front;
        self.escaped_back += o.escaped_back;
        self.escaped_side += o.escaped_side;
        self.no_interaction += o.no_interaction;
        self.absorbed += o.absorbed;
        self.event_cap += o.event_cap;
        self.fermi_sea += o.fermi_sea;
        self.barrier += o.barrier;
        self.phonon_emitted += o.phonon_emitted;
        self.phonon_absorbed += o.phonon_absorbed;
        self.polaron += o.polaron;
    }

    fn report(&self) -> ElectronEnergyBudget {
        let deposited = self.inelastic + self.residual + self.phonon_emitted + self.polaron;
        let escaped = self.escaped_front + self.escaped_back + self.escaped_side;
        let trapped = self.no_interaction + self.absorbed + self.event_cap;
        let source = self.incident + self.fermi_sea + self.phonon_absorbed;
        let sink = deposited + escaped + trapped + self.barrier;
        let relative_imbalance = if source > 0.0 {
            (sink - source).abs() / source
        } else {
            0.0
        };
        ElectronEnergyBudget {
            incident_ev: self.incident,
            deposited_ev: deposited,
            escaped_ev: escaped,
            trapped_ev: trapped,
            inelastic_ev: self.inelastic,
            residual_ev: self.residual,
            escaped_front_ev: self.escaped_front,
            escaped_back_ev: self.escaped_back,
            escaped_side_ev: self.escaped_side,
            no_interaction_ev: self.no_interaction,
            absorbed_ev: self.absorbed,
            event_cap_ev: self.event_cap,
            fermi_sea_ev: self.fermi_sea,
            barrier_ev: self.barrier,
            phonon_emitted_ev: self.phonon_emitted,
            phonon_absorbed_ev: self.phonon_absorbed,
            polaron_ev: self.polaron,
            ionisation_ev: self.inelastic + self.residual,
            lattice_ev: self.phonon_emitted + self.polaron,
            relative_imbalance,
        }
    }
}

/// The channel a deposit is booked to (see
/// [the module docs](self#the-ionisation-and-lattice-channels)).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum Channel {
    Ionisation,
    Lattice,
}

/// One accumulator per channel; the plain total is kept alongside so the
/// existing total maps are unchanged.
#[derive(Debug, Clone, PartialEq)]
struct GridAcc {
    energy_ev: Vec<f64>,
    ionisation_ev: Vec<f64>,
    lattice_ev: Vec<f64>,
    outside_ev: f64,
    outside_ionisation_ev: f64,
    outside_lattice_ev: f64,
}

impl GridAcc {
    fn new(cells: usize) -> Self {
        Self {
            energy_ev: vec![0.0; cells],
            ionisation_ev: vec![0.0; cells],
            lattice_ev: vec![0.0; cells],
            outside_ev: 0.0,
            outside_ionisation_ev: 0.0,
            outside_lattice_ev: 0.0,
        }
    }

    fn add(&mut self, cell: Option<usize>, e: f64, channel: Channel) {
        match cell {
            Some(i) => {
                self.energy_ev[i] += e;
                match channel {
                    Channel::Ionisation => self.ionisation_ev[i] += e,
                    Channel::Lattice => self.lattice_ev[i] += e,
                }
            }
            None => {
                self.outside_ev += e;
                match channel {
                    Channel::Ionisation => self.outside_ionisation_ev += e,
                    Channel::Lattice => self.outside_lattice_ev += e,
                }
            }
        }
    }

    fn merge(&mut self, o: &GridAcc) {
        for (a, b) in self.energy_ev.iter_mut().zip(&o.energy_ev) {
            *a += b;
        }
        for (a, b) in self.ionisation_ev.iter_mut().zip(&o.ionisation_ev) {
            *a += b;
        }
        for (a, b) in self.lattice_ev.iter_mut().zip(&o.lattice_ev) {
            *a += b;
        }
        self.outside_ev += o.outside_ev;
        self.outside_ionisation_ev += o.outside_ionisation_ev;
        self.outside_lattice_ev += o.outside_lattice_ev;
    }
}

#[derive(Debug, Clone, PartialEq)]
struct ClassAcc {
    count: u64,
    energy_ev: f64,
    polar_hist: Histogram,
}

impl ClassAcc {
    fn new(c: &ElectronTallyConfig) -> Self {
        Self {
            count: 0,
            energy_ev: 0.0,
            polar_hist: Histogram::new(c.escape_polar),
        }
    }

    fn merge(&mut self, o: &ClassAcc) {
        self.count += o.count;
        self.energy_ev += o.energy_ev;
        self.polar_hist.merge(&o.polar_hist);
    }

    fn report(&self, histories: u64) -> EmissionClass {
        EmissionClass {
            count: self.count,
            per_primary: per(self.count as f64, histories),
            energy_ev: self.energy_ev,
            mean_energy_ev: per(self.energy_ev, self.count),
            polar_histogram: self.polar_hist.clone(),
        }
    }
}

#[derive(Debug, Clone, PartialEq)]
struct FaceAcc {
    slow: ClassAcc,
    fast: ClassAcc,
    energy_hist: Histogram,
}

impl FaceAcc {
    fn new(c: &ElectronTallyConfig) -> Self {
        Self {
            slow: ClassAcc::new(c),
            fast: ClassAcc::new(c),
            energy_hist: Histogram::new(c.escape_energy),
        }
    }

    fn merge(&mut self, o: &FaceAcc) {
        self.slow.merge(&o.slow);
        self.fast.merge(&o.fast);
        self.energy_hist.merge(&o.energy_hist);
    }

    fn report(&self, histories: u64) -> FaceEmission {
        let count = self.slow.count + self.fast.count;
        FaceEmission {
            count,
            per_primary: per(count as f64, histories),
            energy_ev: self.slow.energy_ev + self.fast.energy_ev,
            energy_histogram: self.energy_hist.clone(),
            slow: self.slow.report(histories),
            fast: self.fast.report(histories),
        }
    }
}

/// Which [`VoxelCell`] field a deposit goes to.
#[derive(Debug, Clone, Copy)]
enum DepositKind {
    /// Left at an inelastic event: `electronic_ev`.
    Electronic,
    /// Given to the lattice (LO-phonon emission, polaron trapping):
    /// `lattice_ev`.
    Lattice,
    /// The energy of an electron that fell below the stopping threshold:
    /// `rest_ev`.
    Rest,
}

impl DepositKind {
    /// The ionisation/lattice channel this deposit is booked to. The
    /// sub-cutoff residual goes to ionisation by convention.
    fn channel(self) -> Channel {
        match self {
            DepositKind::Electronic | DepositKind::Rest => Channel::Ionisation,
            DepositKind::Lattice => Channel::Lattice,
        }
    }
}

/// Sparse per-voxel energy deposition on a [`VoxelTallyGrid`], keyed on the
/// flat voxel index (`BTreeMap`, so merging and export follow one fixed
/// order).
#[derive(Debug, Clone, PartialEq)]
struct VoxelAcc {
    grid: VoxelTallyGrid,
    cells: BTreeMap<u64, VoxelCell>,
    outside: VoxelCell,
}

impl VoxelAcc {
    fn new(grid: VoxelTallyGrid) -> Self {
        Self {
            grid,
            cells: BTreeMap::new(),
            outside: VoxelCell::default(),
        }
    }

    fn add(&mut self, pos: [f64; 3], e: f64, kind: DepositKind) {
        let cell = match self.grid.locate(pos) {
            Some(ijk) => self.cells.entry(self.grid.flat_index(ijk)).or_default(),
            None => &mut self.outside,
        };
        match kind {
            DepositKind::Electronic => cell.electronic_ev += e,
            DepositKind::Lattice => cell.lattice_ev += e,
            DepositKind::Rest => cell.rest_ev += e,
        }
    }

    fn merge(&mut self, o: &VoxelAcc) {
        assert!(
            self.grid == o.grid,
            "merging electron tallies built on different voxel grids"
        );
        self.outside.add(&o.outside);
        for (&flat, cell) in &o.cells {
            self.cells.entry(flat).or_default().add(cell);
        }
    }
}

/// Electron counts at one place, split by the generation of the electron:
/// charge in units of the electron's, per [`ChargeReport`]'s convention.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
pub struct ChargeCounts {
    /// Primaries that terminated here (every terminal fate but escape).
    pub stopped_primary: u64,
    /// Secondaries that terminated here, including a liberated electron
    /// dropped at its creation point because it was below the stopping
    /// threshold.
    pub stopped_secondary: u64,
    /// Holes (liberated electrons) created here by collisions of a primary.
    pub created_by_primary: u64,
    /// Holes created here by collisions of a secondary.
    pub created_by_secondary: u64,
}

impl ChargeCounts {
    /// Electrons that terminated here.
    pub fn stopped(&self) -> u64 {
        self.stopped_primary + self.stopped_secondary
    }

    /// Electrons liberated here (holes left behind).
    pub fn created(&self) -> u64 {
        self.created_by_primary + self.created_by_secondary
    }

    /// Net electron number gained here: stopped minus created.
    pub fn net(&self) -> i64 {
        self.stopped() as i64 - self.created() as i64
    }

    fn merge(&mut self, o: &ChargeCounts) {
        self.stopped_primary += o.stopped_primary;
        self.stopped_secondary += o.stopped_secondary;
        self.created_by_primary += o.created_by_primary;
        self.created_by_secondary += o.created_by_secondary;
    }
}

/// Charge counts per cell of one grid, with the per-history net charge
/// squared and summed for the standard error. Everything is integer, so the
/// merge is exact and independent of the order of the parts.
#[derive(Debug, Clone, PartialEq)]
struct ChargeGrid {
    counts: Vec<ChargeCounts>,
    /// Sum over histories of the squared net charge the history left in the
    /// cell.
    net_sum_sq: Vec<u64>,
    outside: ChargeCounts,
    outside_net_sum_sq: u64,
    /// Per-history scratch: `(cell, net change)`; `usize::MAX` is outside.
    touched: Vec<(usize, i64)>,
}

impl ChargeGrid {
    fn new(cells: usize) -> Self {
        Self {
            counts: vec![ChargeCounts::default(); cells],
            net_sum_sq: vec![0; cells],
            outside: ChargeCounts::default(),
            outside_net_sum_sq: 0,
            touched: Vec::new(),
        }
    }

    /// An electron terminated in `cell`.
    fn stop(&mut self, cell: Option<usize>, secondary: bool) {
        let c = self.slot(cell, 1);
        if secondary {
            c.stopped_secondary += 1;
        } else {
            c.stopped_primary += 1;
        }
    }

    /// A hole was created in `cell` by an electron of the given generation.
    fn create(&mut self, cell: Option<usize>, by_secondary: bool) {
        let c = self.slot(cell, -1);
        if by_secondary {
            c.created_by_secondary += 1;
        } else {
            c.created_by_primary += 1;
        }
    }

    fn slot(&mut self, cell: Option<usize>, delta: i64) -> &mut ChargeCounts {
        self.touched.push((cell.unwrap_or(usize::MAX), delta));
        match cell {
            Some(i) => &mut self.counts[i],
            None => &mut self.outside,
        }
    }

    /// Fold the history's per-cell net charge into the squared sums.
    fn end_history(&mut self) {
        self.touched.sort_unstable_by_key(|t| t.0);
        let mut i = 0;
        while i < self.touched.len() {
            let cell = self.touched[i].0;
            let mut net = 0i64;
            while i < self.touched.len() && self.touched[i].0 == cell {
                net += self.touched[i].1;
                i += 1;
            }
            let sq = (net * net) as u64;
            if cell == usize::MAX {
                self.outside_net_sum_sq += sq;
            } else {
                self.net_sum_sq[cell] += sq;
            }
        }
        self.touched.clear();
    }

    fn merge(&mut self, o: &ChargeGrid) {
        for (a, b) in self.counts.iter_mut().zip(&o.counts) {
            a.merge(b);
        }
        for (a, b) in self.net_sum_sq.iter_mut().zip(&o.net_sum_sq) {
            *a += b;
        }
        self.outside.merge(&o.outside);
        self.outside_net_sum_sq += o.outside_net_sum_sq;
    }

    fn report(&self) -> ChargeCells {
        ChargeCells {
            counts: self.counts.clone(),
            net_sum_sq: self.net_sum_sq.clone(),
            outside: self.outside,
            outside_net_sum_sq: self.outside_net_sum_sq,
        }
    }
}

/// Charge bookkeeping of a tally: per layer, per grid cell, and escapes.
#[derive(Debug, Clone, PartialEq)]
struct ChargeAcc {
    layers: ChargeGrid,
    cartesian: Option<ChargeGrid>,
    cylindrical: Option<ChargeGrid>,
    /// Escapes through each face: `[primary, secondary]`.
    escaped_front: [u64; 2],
    escaped_back: [u64; 2],
    escaped_side: [u64; 2],
}

fn per(x: f64, n: u64) -> f64 {
    if n == 0 {
        0.0
    } else {
        x / n as f64
    }
}

/// Deposition, emission and energy-balance tally for the electron transport
/// loop. See the module docs.
#[derive(Debug, Clone, PartialEq)]
pub struct FullElectronTally {
    config: ElectronTallyConfig,
    cutoff_ev: f64,
    /// Stopping threshold per layer, eV, from the transport.
    thresholds: Vec<f64>,
    /// Whether the transport runs a secondary model (then the `secondary`
    /// hook, not `inelastic`, carries the deposit of an inelastic event).
    secondaries: bool,
    histories: u64,
    fates: FateCounts,
    event_caps: EventCapCounts,
    budget: Budget,
    layer_deposit_ev: Vec<f64>,
    layer_ionisation_ev: Vec<f64>,
    layer_lattice_ev: Vec<f64>,
    cartesian: Option<GridAcc>,
    cylindrical: Option<GridAcc>,
    /// Radial profile in a depth slab with per-history errors.
    psf: Option<RadialAccumulator>,
    /// Dose, depth-dose and transmission accumulators.
    dose: Option<DoseAcc>,
    /// Per-electron deposit in the exposure volumes.
    exposure: Option<ExposureAccumulator>,
    /// Energy-weighted deposition position: `x`, `y`, `z`, `r`.
    generation: [WeightedMoments; 4],
    stop_depth: Moments,
    stop_radial: Moments,
    /// The subset of `stop_depth` / `stop_radial` that are primaries.
    primary_stop_depth: Moments,
    primary_stop_radial: Moments,
    front: FaceAcc,
    back: FaceAcc,
    side: FaceAcc,
    /// Energy deposited per voxel of the optional voxel grid.
    voxel: Option<VoxelAcc>,
    /// Where electrons stopped and holes were created.
    charge: ChargeAcc,
    /// Rate evaluations against each table's energy grid.
    table_coverage: TableCoverageTally,
    /// Per-history scratch, not part of the result and cleared at every
    /// history boundary: the last state seen of the current electron, for
    /// [`Fate::EventCap`] (`None` once it ended through a hook).
    last: Option<ElectronState>,
    /// The primary's last state, set aside when the first secondary begins.
    primary_last: Option<ElectronState>,
    /// Whether the history has reached its secondaries.
    in_secondary: bool,
    /// Whether a secondary of the current history was cut off by the event
    /// cap (so the history is counted once in
    /// [`EventCapCounts::affected_histories`]).
    history_capped: bool,
}

impl FullElectronTally {
    /// Empty tally for runs of `transport` (its per-layer stopping thresholds
    /// classify `stopped` events, its secondary model decides which hook
    /// carries an inelastic deposit, its stack sets the number of layers).
    pub fn new<G: Geometry>(
        transport: &Transport<G>,
        config: ElectronTallyConfig,
    ) -> Result<Self, ElectronTallyError> {
        if !(config.se_bse_split_ev.is_finite() && config.se_bse_split_ev >= 0.0) {
            return Err(ElectronTallyError::Split(config.se_bse_split_ev));
        }
        let p = config.escape_polar;
        if p.lo < 0.0 || p.hi > std::f64::consts::PI {
            return Err(ElectronTallyError::PolarRange { lo: p.lo, hi: p.hi });
        }
        let cartesian = match config.cartesian {
            None => None,
            Some(g) => Some(GridAcc::new(
                g.voxels()
                    .ok_or(ElectronTallyError::GridTooLarge("Cartesian"))?,
            )),
        };
        let cylindrical = match config.cylindrical {
            None => None,
            Some(g) => {
                if g.r.lo < 0.0 {
                    return Err(ElectronTallyError::NegativeRadius(g.r.lo));
                }
                Some(GridAcc::new(
                    g.cells()
                        .ok_or(ElectronTallyError::GridTooLarge("cylindrical"))?,
                ))
            }
        };
        let psf = config.psf.map(RadialAccumulator::new).transpose()?;
        let dose = match config.dose {
            None => None,
            Some(d) => {
                let stack = transport.geometry().as_stack().ok_or_else(|| {
                    DoseError("the dose tally needs a layered stack target".into())
                })?;
                Some(DoseAcc::new(stack, d, config.escape_energy)?)
            }
        };
        let exposure = config
            .exposure
            .clone()
            .map(ExposureAccumulator::new)
            .transpose()?;
        let (front, back, side) = (
            FaceAcc::new(&config),
            FaceAcc::new(&config),
            FaceAcc::new(&config),
        );
        let charge = ChargeAcc {
            layers: ChargeGrid::new(transport.n_materials()),
            cartesian: cartesian
                .as_ref()
                .map(|a| ChargeGrid::new(a.energy_ev.len())),
            cylindrical: cylindrical
                .as_ref()
                .map(|a| ChargeGrid::new(a.energy_ev.len())),
            escaped_front: [0; 2],
            escaped_back: [0; 2],
            escaped_side: [0; 2],
        };
        Ok(Self {
            config,
            cutoff_ev: transport.config().cutoff_ev,
            thresholds: transport.stopping_thresholds_ev().to_vec(),
            secondaries: transport.config().secondaries != SecondaryModel::Off,
            histories: 0,
            fates: FateCounts::default(),
            event_caps: EventCapCounts::default(),
            budget: Budget::default(),
            layer_deposit_ev: vec![0.0; transport.n_materials()],
            layer_ionisation_ev: vec![0.0; transport.n_materials()],
            layer_lattice_ev: vec![0.0; transport.n_materials()],
            cartesian,
            cylindrical,
            psf,
            dose,
            exposure,
            generation: [WeightedMoments::default(); 4],
            stop_depth: Moments::new(),
            stop_radial: Moments::new(),
            primary_stop_depth: Moments::new(),
            primary_stop_radial: Moments::new(),
            front,
            back,
            side,
            voxel: None,
            charge,
            table_coverage: TableCoverageTally::new(transport),
            last: None,
            primary_last: None,
            in_secondary: false,
            history_capped: false,
        })
    }

    /// The configuration.
    pub fn config(&self) -> &ElectronTallyConfig {
        &self.config
    }

    /// Primary histories tallied so far.
    pub fn histories(&self) -> u64 {
        self.histories
    }

    /// Also tally the deposited energy per voxel of `grid`, for
    /// [`FullElectronTally::voxel_report`]. The grid is independent of the
    /// target: use [`VoxelTallyGrid::from_voxel_grid`] to tally on the voxels
    /// of a [`crate::geometry::VoxelGrid`] target.
    pub fn with_voxel_grid(mut self, grid: VoxelTallyGrid) -> Self {
        self.voxel = Some(VoxelAcc::new(grid));
        self
    }

    /// The per-voxel deposition as a [`VoxelReport`] (CSV and VTK export of
    /// [`crate::tally::voxel`]), or `None` without
    /// [`FullElectronTally::with_voxel_grid`].
    ///
    /// Each voxel's energy is split by how it was deposited: `electronic_ev`
    /// is the energy left at inelastic events (the
    /// [`ElectronEnergyBudget::inelastic_ev`] channel), `lattice_ev` the
    /// LO-phonon emission and polaron-trapping energy, and `rest_ev` the
    /// energy of electrons that fell below the stopping threshold
    /// ([`ElectronEnergyBudget::residual_ev`]); their sum over all voxels
    /// and `outside` is [`ElectronEnergyBudget::deposited_ev`].
    /// `implanted` and `displacements` are zero: the electron tally records
    /// neither, and the per-species defect fields (`species_z`, `defects`,
    /// `outside_defects`) are empty: the electron tally has no species axis.
    /// Deposits are point events at the collision (as for
    /// [`CartesianGrid`]), so a flight crossing voxel faces is not split.
    pub fn voxel_report(&self) -> Option<VoxelReport> {
        self.voxel.as_ref().map(|v| VoxelReport {
            grid: v.grid,
            histories: self.histories,
            voxels: v
                .cells
                .iter()
                .map(|(&flat, &cell)| VoxelRecord {
                    index: v.grid.unflatten(flat),
                    cell,
                    defects: Vec::new(),
                })
                .collect(),
            outside: v.outside,
            species_z: Vec::new(),
            outside_defects: Vec::new(),
        })
    }

    /// Every energy deposit goes through here. `kind` picks the voxel field
    /// and, through [`DepositKind::channel`], the ionisation or lattice
    /// channel of the layer and grid maps.
    fn deposit(&mut self, at: &ElectronState, e: f64, kind: DepositKind) {
        if e <= 0.0 {
            return;
        }
        let pos = at.pos;
        if let Some(v) = &mut self.voxel {
            v.add(pos, e, kind);
        }
        self.layer_deposit_ev[at.layer] += e;
        let channel = kind.channel();
        match channel {
            Channel::Ionisation => self.layer_ionisation_ev[at.layer] += e,
            Channel::Lattice => self.layer_lattice_ev[at.layer] += e,
        }
        if let (Some(acc), Some(g)) = (&mut self.cartesian, &self.config.cartesian) {
            acc.add(g.locate(pos), e, channel);
        }
        if let (Some(acc), Some(g)) = (&mut self.cylindrical, &self.config.cylindrical) {
            acc.add(g.locate(pos), e, channel);
        }
        if let Some(p) = &mut self.psf {
            p.deposit(pos, e);
        }
        if let Some(d) = &mut self.dose {
            d.deposit(at.layer, pos[0], e);
        }
        if let Some(x) = &mut self.exposure {
            x.deposit(pos, e);
        }
        let r = pos[1].hypot(pos[2]);
        for (m, v) in self.generation.iter_mut().zip([pos[0], pos[1], pos[2], r]) {
            m.push(v, e);
        }
    }

    /// An electron of the current generation terminated at `at`.
    fn charge_stop(&mut self, at: &ElectronState) {
        self.charge_stop_gen(at, self.in_secondary);
    }

    fn charge_stop_gen(&mut self, at: &ElectronState, secondary: bool) {
        let c = &mut self.charge;
        c.layers.stop(Some(at.layer), secondary);
        if let (Some(a), Some(g)) = (&mut c.cartesian, &self.config.cartesian) {
            a.stop(g.locate(at.pos), secondary);
        }
        if let (Some(a), Some(g)) = (&mut c.cylindrical, &self.config.cylindrical) {
            a.stop(g.locate(at.pos), secondary);
        }
    }

    /// A hole was left at `at` by a collision of the current electron.
    fn charge_create(&mut self, at: &ElectronState) {
        let by_secondary = self.in_secondary;
        let c = &mut self.charge;
        c.layers.create(Some(at.layer), by_secondary);
        if let (Some(a), Some(g)) = (&mut c.cartesian, &self.config.cartesian) {
            a.create(g.locate(at.pos), by_secondary);
        }
        if let (Some(a), Some(g)) = (&mut c.cylindrical, &self.config.cylindrical) {
            a.create(g.locate(at.pos), by_secondary);
        }
    }

    fn charge_end_history(&mut self) {
        let c = &mut self.charge;
        c.layers.end_history();
        if let Some(a) = &mut c.cartesian {
            a.end_history();
        }
        if let Some(a) = &mut c.cylindrical {
            a.end_history();
        }
    }

    fn charge_report(&self) -> ChargeReport {
        let c = &self.charge;
        let mut total = ChargeCounts::default();
        for l in &c.layers.counts {
            total.merge(l);
        }
        let esc = |e: [u64; 2]| EscapedCharge {
            primary: e[0],
            secondary: e[1],
        };
        let escaped = c
            .escaped_front
            .iter()
            .chain(&c.escaped_back)
            .chain(&c.escaped_side)
            .sum::<u64>();
        ChargeReport {
            incident: self.histories,
            total,
            per_layer: c.layers.report(),
            escaped_front: esc(c.escaped_front),
            escaped_back: esc(c.escaped_back),
            escaped_side: esc(c.escaped_side),
            imbalance: self.histories as i64 + total.created() as i64
                - total.stopped() as i64
                - escaped as i64,
            cartesian: self
                .config
                .cartesian
                .zip(c.cartesian.as_ref())
                .map(|(grid, a)| CartesianCharge {
                    grid,
                    cells: a.report(),
                }),
            cylindrical: self
                .config
                .cylindrical
                .zip(c.cylindrical.as_ref())
                .map(|(grid, a)| CylindricalCharge {
                    grid,
                    cells: a.report(),
                }),
        }
    }

    /// The plain-data result.
    pub fn report(&self) -> ElectronReport {
        let h = self.histories;
        let front = self.front.report(h);
        let back = self.back.report(h);
        let yields = Yields {
            backscatter_eta: front.fast.per_primary,
            secondary_delta: front.slow.per_primary,
            total_sigma: per((front.fast.count + front.slow.count) as f64, h),
            transmitted_fast: back.fast.per_primary,
            transmitted_slow: back.slow.per_primary,
        };
        let [gx, gy, gz, gr] = &self.generation;
        let generation_volume = (gx.weight > 0.0).then(|| GenerationVolume {
            energy_ev: gx.weight,
            mean_m: [gx.mean, gy.mean, gz.mean],
            std_dev_m: [
                gx.variance().sqrt(),
                gy.variance().sqrt(),
                gz.variance().sqrt(),
            ],
            mean_radius_m: gr.mean,
            radius_std_dev_m: gr.variance().sqrt(),
            rms_radius_m: (gy.variance() + gy.mean * gy.mean + gz.variance() + gz.mean * gz.mean)
                .sqrt(),
        });
        let deposition =
            DepositionReport {
                per_layer_ev: self.layer_deposit_ev.clone(),
                per_layer_ionisation_ev: self.layer_ionisation_ev.clone(),
                per_layer_lattice_ev: self.layer_lattice_ev.clone(),
                cartesian: self
                    .config
                    .cartesian
                    .zip(self.cartesian.as_ref())
                    .map(|(grid, a)| CartesianDeposition {
                        grid,
                        energy_ev: a.energy_ev.clone(),
                        outside_ev: a.outside_ev,
                        ionisation_ev: a.ionisation_ev.clone(),
                        lattice_ev: a.lattice_ev.clone(),
                        outside_ionisation_ev: a.outside_ionisation_ev,
                        outside_lattice_ev: a.outside_lattice_ev,
                    }),
                cylindrical: self.config.cylindrical.zip(self.cylindrical.as_ref()).map(
                    |(grid, a)| CylindricalDeposition {
                        grid,
                        energy_ev: a.energy_ev.clone(),
                        outside_ev: a.outside_ev,
                        ionisation_ev: a.ionisation_ev.clone(),
                        lattice_ev: a.lattice_ev.clone(),
                        outside_ionisation_ev: a.outside_ionisation_ev,
                        outside_lattice_ev: a.outside_lattice_ev,
                    },
                ),
                psf: self.psf.as_ref().map(RadialAccumulator::profile),
                exposure: self.exposure.as_ref().map(ExposureAccumulator::report),
            };
        let dose = self.dose.as_ref().map(|d| d.report(self.budget.incident));
        ElectronReport {
            histories: h,
            metadata: ElectronTallyMetadata {
                se_bse_split_ev: self.config.se_bse_split_ev,
                se_bse_split_rule: SE_BSE_SPLIT_RULE.to_string(),
                se_bse_split_source: SE_BSE_SPLIT_SOURCE.to_string(),
                cutoff_ev: self.cutoff_ev,
                stopping_threshold_ev: self.thresholds.clone(),
                layers: self.layer_deposit_ev.len(),
                config: self.config.clone(),
            },
            fates: self.fates,
            event_caps: self.event_caps,
            budget: self.budget.report(),
            yields,
            front,
            back,
            side: (self.side.slow.count + self.side.fast.count > 0).then(|| self.side.report(h)),
            deposition,
            dose,
            generation_volume,
            charge: self.charge_report(),
            stopping_points: StoppingPoints {
                stopped: self.stop_depth.n,
                depth: self.stop_depth.summary(),
                radial: self.stop_radial.summary(),
                primaries: PrimaryStoppingPoints {
                    stopped: self.primary_stop_depth.n,
                    depth: self.primary_stop_depth.summary(),
                    radial: self.primary_stop_radial.summary(),
                },
            },
            table_coverage: self.table_coverage.layers().to_vec(),
        }
    }
}

impl ElectronTally for FullElectronTally {
    fn begin_history(&mut self, index: u64, start: &ElectronState) {
        if let Some(x) = &mut self.exposure {
            x.begin_history(index, [start.pos[1], start.pos[2]]);
        }
        self.budget.incident += start.energy_ev;
        self.last = Some(*start);
        self.primary_last = None;
        self.in_secondary = false;
        self.history_capped = false;
    }

    fn step(&mut self, _from: [f64; 3], end: &ElectronState, _length_m: f64) {
        self.last = Some(*end);
    }

    fn table_lookup(&mut self, at: &ElectronState, channel: TableChannel, coverage: GridCoverage) {
        self.table_coverage.record(at.layer, channel, coverage);
    }

    fn elastic(&mut self, after: &ElectronState, _theta: f64) {
        self.last = Some(*after);
    }

    fn inelastic(&mut self, after: &ElectronState, w_ev: f64) {
        self.last = Some(*after);
        // With a secondary model the `secondary` hook that follows carries
        // the part of `w_ev` left in the solid.
        if !self.secondaries {
            self.budget.inelastic += w_ev;
            self.deposit(after, w_ev, DepositKind::Electronic);
        }
    }

    fn secondary(
        &mut self,
        primary: &ElectronState,
        event: &SecondaryEvent,
        created: Option<&ElectronState>,
    ) {
        self.last = Some(*primary);
        if event.liberated {
            // A hole is left where the electron was lifted out. A liberated
            // electron that is not followed was below the stopping threshold
            // and is dropped here, so it stops where its hole is.
            self.charge_create(primary);
            if created.is_none() {
                self.charge_stop_gen(primary, true);
            }
        }
        // `W = secondary + binding + deposited` (module docs): the created
        // secondary is counted where it ends, a positive binding stays in
        // the solid, a negative one is energy the liberated electron brings.
        let local = event.deposited_ev + event.binding_ev.max(0.0);
        self.budget.inelastic += local;
        self.budget.fermi_sea += (-event.binding_ev).max(0.0);
        self.deposit(primary, local, DepositKind::Electronic);
    }

    fn phonon(
        &mut self,
        after: &ElectronState,
        kind: PhononEvent,
        hbar_omega_ev: f64,
        _theta: f64,
    ) {
        self.last = Some(*after);
        match kind {
            PhononEvent::Emission => {
                self.budget.phonon_emitted += hbar_omega_ev;
                self.deposit(after, hbar_omega_ev, DepositKind::Lattice);
            }
            PhononEvent::Absorption => self.budget.phonon_absorbed += hbar_omega_ev,
        }
    }

    fn polaron_trapped(&mut self, at: &ElectronState) {
        self.last = None;
        self.charge_stop(at);
        self.budget.polaron += at.energy_ev;
        self.deposit(at, at.energy_ev, DepositKind::Lattice);
    }

    fn interface(&mut self, at: &ElectronState, from_layer: usize, to_layer: usize) {
        self.last = Some(*at);
        if let Some(d) = &mut self.dose {
            d.crossing(from_layer, to_layer, at.energy_ev, self.in_secondary);
        }
    }

    fn barrier(&mut self, at: &ElectronState, _boundary: Boundary, delta_u_ev: f64) {
        self.last = Some(*at);
        self.budget.barrier -= delta_u_ev;
    }

    fn reflected(&mut self, at: &ElectronState, _boundary: Boundary) {
        self.last = Some(*at);
    }

    fn begin_secondary(&mut self, start: &ElectronState, _generation: u32) {
        if !self.in_secondary {
            // The primary has ended; keep its last state for `end_history`.
            self.primary_last = self.last.take();
            self.in_secondary = true;
        }
        self.last = Some(*start);
    }

    fn end_secondary(&mut self, fate: Fate) {
        if fate == Fate::EventCap {
            if let Some(s) = self.last {
                self.budget.event_cap += s.energy_ev;
            }
            if let Some(s) = self.last {
                self.charge_stop_gen(&s, true);
            }
            self.event_caps.secondary_tracks += 1;
            self.history_capped = true;
        }
        self.last = None;
    }

    fn stopped(&mut self, at: &ElectronState) {
        self.last = None;
        // Below the threshold or trapped without an interaction, the
        // electron ends inside the target either way.
        self.charge_stop(at);
        let e = at.energy_ev;
        if e < self.thresholds[at.layer] {
            self.budget.residual += e;
            self.deposit(at, e, DepositKind::Rest);
            let r = at.pos[1].hypot(at.pos[2]);
            self.stop_depth.push(at.pos[0]);
            self.stop_radial.push(r);
            if !self.in_secondary {
                // The primary is followed to its end before any secondary
                // begins, so this is the primary stopping.
                self.primary_stop_depth.push(at.pos[0]);
                self.primary_stop_radial.push(r);
            }
        } else {
            self.budget.no_interaction += e;
        }
    }

    fn escaped(&mut self, at: &ElectronState, face: Face) {
        self.last = None;
        let e = at.energy_ev;
        let polar = at.dir[0].abs().min(1.0).acos();
        let g = usize::from(self.in_secondary);
        let acc = match face {
            Face::Front => {
                self.budget.escaped_front += e;
                self.charge.escaped_front[g] += 1;
                &mut self.front
            }
            Face::Back => {
                self.budget.escaped_back += e;
                self.charge.escaped_back[g] += 1;
                &mut self.back
            }
            Face::Side => {
                self.budget.escaped_side += e;
                self.charge.escaped_side[g] += 1;
                &mut self.side
            }
        };
        acc.energy_hist.fill(e);
        let class = if e < self.config.se_bse_split_ev {
            &mut acc.slow
        } else {
            &mut acc.fast
        };
        class.count += 1;
        class.energy_ev += e;
        class.polar_hist.fill(polar);
    }

    fn absorbed(&mut self, at: &ElectronState) {
        self.last = None;
        self.charge_stop(at);
        self.budget.absorbed += at.energy_ev;
    }

    fn end_history(&mut self, _index: u64, fate: Fate) {
        self.histories += 1;
        if let Some(p) = &mut self.psf {
            p.end_history();
        }
        if let Some(d) = &mut self.dose {
            d.end_history();
        }
        if let Some(x) = &mut self.exposure {
            x.end_history();
        }
        match fate {
            Fate::Stopped => self.fates.stopped += 1,
            Fate::Escaped(Face::Front) => self.fates.escaped_front += 1,
            Fate::Escaped(Face::Back) => self.fates.escaped_back += 1,
            Fate::Escaped(Face::Side) => self.fates.escaped_side += 1,
            Fate::Absorbed => self.fates.absorbed += 1,
            Fate::Trapped => self.fates.trapped += 1,
            Fate::PolaronTrapped => self.fates.polaron_trapped += 1,
            Fate::EventCap => {
                self.fates.event_capped += 1;
                let primary = if self.in_secondary {
                    self.primary_last
                } else {
                    self.last
                };
                if let Some(s) = primary {
                    self.budget.event_cap += s.energy_ev;
                    self.charge_stop_gen(&s, false);
                }
            }
        }
        self.charge_end_history();
        if self.history_capped || fate == Fate::EventCap {
            // Once per history, however many of its electrons were capped.
            self.event_caps.affected_histories += 1;
        }
        self.last = None;
        self.primary_last = None;
        self.in_secondary = false;
        self.history_capped = false;
    }

    /// Field-wise merge.
    ///
    /// # Panics
    /// If `other` was built with a different configuration, cutoff,
    /// thresholds, secondary model or stack.
    fn merge(&mut self, o: Self) {
        assert!(
            self.config == o.config
                && self.cutoff_ev == o.cutoff_ev
                && self.thresholds == o.thresholds
                && self.secondaries == o.secondaries
                && self.layer_deposit_ev.len() == o.layer_deposit_ev.len(),
            "merging electron tallies built for different runs"
        );
        self.histories += o.histories;
        self.fates.merge(&o.fates);
        self.event_caps.merge(&o.event_caps);
        self.budget.merge(&o.budget);
        for (a, b) in self.layer_deposit_ev.iter_mut().zip(&o.layer_deposit_ev) {
            *a += b;
        }
        for (a, b) in self
            .layer_ionisation_ev
            .iter_mut()
            .zip(&o.layer_ionisation_ev)
        {
            *a += b;
        }
        for (a, b) in self.layer_lattice_ev.iter_mut().zip(&o.layer_lattice_ev) {
            *a += b;
        }
        if let (Some(a), Some(b)) = (&mut self.cartesian, &o.cartesian) {
            a.merge(b);
        }
        if let (Some(a), Some(b)) = (&mut self.cylindrical, &o.cylindrical) {
            a.merge(b);
        }
        if let (Some(a), Some(b)) = (&mut self.psf, &o.psf) {
            a.merge(b);
        }
        if let (Some(a), Some(b)) = (&mut self.dose, &o.dose) {
            a.merge(b);
        }
        if let (Some(a), Some(b)) = (&mut self.exposure, &o.exposure) {
            a.merge(b);
        }
        for (a, b) in self.generation.iter_mut().zip(&o.generation) {
            a.merge(b);
        }
        self.stop_depth.merge(&o.stop_depth);
        self.stop_radial.merge(&o.stop_radial);
        self.primary_stop_depth.merge(&o.primary_stop_depth);
        self.primary_stop_radial.merge(&o.primary_stop_radial);
        self.front.merge(&o.front);
        self.back.merge(&o.back);
        self.side.merge(&o.side);
        match (&mut self.voxel, &o.voxel) {
            (Some(a), Some(b)) => a.merge(b),
            (None, None) => {}
            _ => panic!("merging electron tallies built for different runs"),
        }
        self.charge.layers.merge(&o.charge.layers);
        if let (Some(a), Some(b)) = (&mut self.charge.cartesian, &o.charge.cartesian) {
            a.merge(b);
        }
        if let (Some(a), Some(b)) = (&mut self.charge.cylindrical, &o.charge.cylindrical) {
            a.merge(b);
        }
        for g in 0..2 {
            self.charge.escaped_front[g] += o.charge.escaped_front[g];
            self.charge.escaped_back[g] += o.charge.escaped_back[g];
            self.charge.escaped_side[g] += o.charge.escaped_side[g];
        }
        self.table_coverage.merge_from(&o.table_coverage);
    }
}

/// Everything a [`FullElectronTally`] measured, as plain data.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct ElectronReport {
    /// Primary histories.
    pub histories: u64,
    /// The tally settings and conventions.
    pub metadata: ElectronTallyMetadata,
    /// How the primary histories ended.
    pub fates: FateCounts,
    /// Electrons cut off by the collision cap: capped secondary tracks and
    /// histories with any capped electron. A numerical diagnostic; all zero
    /// in reports written before it existed.
    #[serde(default)]
    pub event_caps: EventCapCounts,
    /// Summed energy balance of all histories, eV.
    pub budget: ElectronEnergyBudget,
    /// Emission yields per primary.
    pub yields: Yields,
    /// Electrons leaving the front face (`x = 0`).
    pub front: FaceEmission,
    /// Electrons leaving the back face of a finite stack.
    pub back: FaceEmission,
    /// Electrons leaving through a lateral face of a voxel target; `None`
    /// when none did (always, for a stack).
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub side: Option<FaceEmission>,
    /// Where the energy was deposited.
    pub deposition: DepositionReport,
    /// Dose per layer, depth-dose and interface transmission, if configured
    /// ([`ElectronTallyConfig::dose`]).
    #[serde(default)]
    pub dose: Option<DoseReport>,
    /// Energy-weighted moments of the deposition position; `None` if nothing
    /// was deposited.
    pub generation_volume: Option<GenerationVolume>,
    /// Where electrons stopped and where holes were left, and the balance
    /// against the incident electrons. Empty in reports written before it
    /// existed.
    #[serde(default)]
    pub charge: ChargeReport,
    /// Where electrons fell below the stopping threshold.
    pub stopping_points: StoppingPoints,
    /// Per layer, the rate evaluations of the elastic and inelastic tables
    /// inside and beyond each table's energy grid, with the grid bounds
    /// (primaries and secondaries together). These count evaluations, not
    /// collisions; see [`super::table_coverage`]. Empty in reports written
    /// before it existed.
    #[serde(default)]
    pub table_coverage: Vec<LayerTableCoverage>,
}

/// Where the electrons of a run came to rest and where holes were left, as
/// integer counts, and the balance against the incident electrons.
///
/// Charge is counted in units of the electron's, per electron number:
///
/// - every electron that ends inside the target counts `+1` where it ends:
///   every terminal fate except escape (stopped below the layer's threshold,
///   trapped without an interaction, absorbed at the back face, polaron
///   trapped, cut off by the collision cap);
/// - every liberated electron (a secondary-model event with
///   `liberated = true`) counts `-1` at the event, the hole it leaves, whether
///   or not the electron is followed. A liberated electron below the stopping
///   threshold is not followed, it is dropped and counted `+1` where it was
///   created, so its `-1` and `+1` cancel;
/// - an electron leaving through a face carries its charge out and is counted
///   in `escaped_front` / `escaped_back`.
///
/// The balance is exact in integers: `incident + created = stopped +
/// escaped`, so `imbalance` is zero on every run and the net electron number
/// deposited in the target, `stopped - created`, equals `incident - escaped`.
/// Generation: `primary` is the incident electron, `secondary` every electron
/// created by the secondary model. A secondary below the cutoff is not
/// followed and is therefore counted where it is dropped (see
/// [`ElectronTallyMetadata::cutoff_ev`]).
#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
pub struct ChargeReport {
    /// Incident electrons (primary histories).
    pub incident: u64,
    /// Totals over the target.
    pub total: ChargeCounts,
    /// Per layer.
    pub per_layer: ChargeCells,
    /// Electrons that left through the front face.
    pub escaped_front: EscapedCharge,
    /// Electrons that left through the back face.
    pub escaped_back: EscapedCharge,
    /// Electrons that left through the lateral boundary of a non-slab
    /// geometry.
    pub escaped_side: EscapedCharge,
    /// `incident + created - stopped - escaped`; zero if every electron was
    /// counted exactly once.
    pub imbalance: i64,
    /// The Cartesian grid, if the run has one.
    pub cartesian: Option<CartesianCharge>,
    /// The cylindrical grid, if the run has one.
    pub cylindrical: Option<CylindricalCharge>,
}

/// Electrons that left through one face, by generation.
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)]
pub struct EscapedCharge {
    /// Primaries.
    pub primary: u64,
    /// Secondaries.
    pub secondary: u64,
}

/// Charge counts per cell (or per layer) with the data for the standard
/// error of the net charge per primary.
#[derive(Debug, Clone, Default, PartialEq, Eq, Serialize, Deserialize)]
pub struct ChargeCells {
    /// Counts per cell.
    pub counts: Vec<ChargeCounts>,
    /// Per cell, the sum over histories of the square of the net charge
    /// (stopped minus created) the history left there.
    pub net_sum_sq: Vec<u64>,
    /// Counts outside the grid.
    pub outside: ChargeCounts,
    /// Squared net charge sum outside the grid.
    pub outside_net_sum_sq: u64,
}

impl ChargeCells {
    /// Net charge per incident electron in cell `i` and its standard error
    /// (of the mean over `histories`).
    pub fn net_per_primary(&self, i: usize, histories: u64) -> (f64, f64) {
        net_mean_se(self.counts[i].net(), self.net_sum_sq[i], histories)
    }
}

/// Mean and standard error of a per-history integer quantity from its sum
/// and sum of squares.
pub fn net_mean_se(sum: i64, sum_sq: u64, n: u64) -> (f64, f64) {
    if n == 0 {
        return (0.0, 0.0);
    }
    let nf = n as f64;
    let mean = sum as f64 / nf;
    if n == 1 {
        return (mean, 0.0);
    }
    let var = ((sum_sq as f64 - nf * mean * mean) / (nf - 1.0)).max(0.0);
    (mean, (var / nf).sqrt())
}

/// Charge on the Cartesian grid, in [`CartesianGrid::index`] order.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct CartesianCharge {
    /// The grid.
    pub grid: CartesianGrid,
    /// Per voxel.
    pub cells: ChargeCells,
}

/// Charge on the cylindrical `(r, depth)` grid, in
/// [`CylindricalGrid::index`] order.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct CylindricalCharge {
    /// The grid.
    pub grid: CylindricalGrid,
    /// Per cell.
    pub cells: ChargeCells,
}

/// The settings and conventions a report was made with.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct ElectronTallyMetadata {
    /// SE/BSE split, eV.
    pub se_bse_split_ev: f64,
    /// How an escaping electron is classed against the split.
    pub se_bse_split_rule: String,
    /// Where the split convention comes from.
    pub se_bse_split_source: String,
    /// Transport cutoff, eV, as configured (measured from the band bottom or
    /// the vacuum level, see the run's
    /// [`RunMetadata`](crate::electron::transport::RunMetadata)).
    pub cutoff_ev: f64,
    /// Stopping threshold per layer, eV: a `stopped` electron below its
    /// layer's threshold is a deposit, at or above it a trapped electron.
    pub stopping_threshold_ev: Vec<f64>,
    /// Layers in the stack.
    pub layers: usize,
    /// The full tally configuration (grids and spectrum binnings).
    pub config: ElectronTallyConfig,
}

/// How the primary histories ended (one per history, from its
/// [`Fate`]).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
pub struct FateCounts {
    /// Fell below the stopping threshold in the target.
    pub stopped: u64,
    /// Left through the front face.
    pub escaped_front: u64,
    /// Left through the back face.
    pub escaped_back: u64,
    /// Left through a lateral face of a voxel target.
    #[serde(default)]
    pub escaped_side: u64,
    /// Absorbed at the back face.
    pub absorbed: u64,
    /// No interaction available and no face to reach.
    pub trapped: u64,
    /// Cut off by the collision cap.
    pub event_capped: u64,
    /// Trapped as a polaron (an opt-in insulator channel).
    #[serde(default)]
    pub polaron_trapped: u64,
}

impl FateCounts {
    fn merge(&mut self, o: &FateCounts) {
        self.stopped += o.stopped;
        self.escaped_front += o.escaped_front;
        self.escaped_back += o.escaped_back;
        self.escaped_side += o.escaped_side;
        self.absorbed += o.absorbed;
        self.trapped += o.trapped;
        self.event_capped += o.event_capped;
        self.polaron_trapped += o.polaron_trapped;
    }
}

/// Electrons cut off by the collision cap ([`Fate::EventCap`]), a numerical
/// truncation diagnostic rather than a physical fate. The primaries' own caps
/// are [`FateCounts::event_capped`] (one fate per primary); these counts add
/// the secondaries. Both are plain counts (tracks, histories), not energies;
/// the energy the capped electrons still carried is
/// [`ElectronEnergyBudget::event_cap_ev`].
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
pub struct EventCapCounts {
    /// Secondary tracks cut off by the cap, one per capped secondary.
    #[serde(default)]
    pub secondary_tracks: u64,
    /// Primary histories in which the primary or at least one secondary was
    /// cut off by the cap, each counted once.
    #[serde(default)]
    pub affected_histories: u64,
}

impl EventCapCounts {
    fn merge(&mut self, o: &EventCapCounts) {
        self.secondary_tracks += o.secondary_tracks;
        self.affected_histories += o.affected_histories;
    }
}

/// The energy balance, eV, summed over all histories and every electron of
/// them (see [the module docs](self#the-energy-balance)):
/// `incident + fermi_sea + phonon_absorbed = deposited + escaped + trapped +
/// barrier`, with
///
/// - `deposited = inelastic + residual + phonon_emitted + polaron`: energy
///   left in the solid at inelastic events, plus the energy electrons had
///   left when they fell below the stopping threshold, plus the `ħω` given to
///   the lattice at LO-phonon emissions, plus the whole remaining energy of
///   electrons trapped as polarons;
/// - `phonon_absorbed`: the `ħω` electrons took from the lattice at LO-phonon
///   absorptions (a source; zero without the insulator channels);
/// - `escaped = escaped_front + escaped_back + escaped_side`;
/// - `trapped = no_interaction + absorbed + event_cap`: energy still carried
///   by electrons that ended in the target without falling below the
///   threshold (no interaction available, absorbed at the back face, or cut
///   off by the collision cap);
/// - `fermi_sea`: the kinetic energy liberated secondaries already had (a
///   source; zero without a secondary model);
/// - `barrier`: the kinetic energy taken by the potential steps at the faces,
///   the sum of `-ΔU` (signed; zero with transparent faces).
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct ElectronEnergyBudget {
    /// Energy of the primaries.
    pub incident_ev: f64,
    /// Deposited in the target, measured from the band bottom. With a
    /// secondary model in a layer with a Fermi energy this includes the
    /// Fermi-sea energy of liberated conduction electrons (`fermi_sea_ev`),
    /// so it can exceed the energy the beam imparted; the net energy imparted
    /// is `incident_ev - escaped_ev = deposited_ev + trapped_ev + barrier_ev -
    /// fermi_sea_ev - phonon_absorbed_ev` (see
    /// [the module docs](self#deposits-and-the-energy-imparted)).
    pub deposited_ev: f64,
    /// Carried out of the target.
    pub escaped_ev: f64,
    /// Carried by electrons that ended in the target at or above the
    /// stopping threshold.
    pub trapped_ev: f64,
    /// Energy left in the solid at inelastic events (part of
    /// `deposited_ev`): the loss `W` without a secondary model; with one,
    /// [`SecondaryEvent::deposited_ev`] plus any positive
    /// [`SecondaryEvent::binding_ev`].
    pub inelastic_ev: f64,
    /// Remaining energy of electrons that fell below the stopping threshold
    /// (part of `deposited_ev`).
    pub residual_ev: f64,
    /// Out of the front face (part of `escaped_ev`).
    pub escaped_front_ev: f64,
    /// Out of the back face (part of `escaped_ev`).
    pub escaped_back_ev: f64,
    /// Out of a lateral face of a voxel target (part of `escaped_ev`); zero
    /// for a stack.
    #[serde(default)]
    pub escaped_side_ev: f64,
    /// Electrons with no interaction and no face to reach (part of
    /// `trapped_ev`).
    pub no_interaction_ev: f64,
    /// Electrons absorbed at the back face (part of `trapped_ev`).
    pub absorbed_ev: f64,
    /// Electrons cut off by the collision cap (part of `trapped_ev`).
    pub event_cap_ev: f64,
    /// Energy the liberated secondaries brought with them, the sum of
    /// `-binding_ev` over events where it is negative (a source).
    pub fermi_sea_ev: f64,
    /// Kinetic energy taken by the potential steps, the sum of `-ΔU` over
    /// face transmissions (signed).
    pub barrier_ev: f64,
    /// `ħω` given to the lattice by LO-phonon emission (part of
    /// `deposited_ev`).
    #[serde(default)]
    pub phonon_emitted_ev: f64,
    /// `ħω` taken from the lattice by LO-phonon absorption (a source).
    #[serde(default)]
    pub phonon_absorbed_ev: f64,
    /// Remaining energy of electrons trapped as polarons (part of
    /// `deposited_ev`).
    #[serde(default)]
    pub polaron_ev: f64,
    /// Ionisation channel of `deposited_ev`: `inelastic_ev + residual_ev`
    /// (the residual is placed here by convention, see
    /// [the channel split](self#the-ionisation-and-lattice-channels)).
    #[serde(default)]
    pub ionisation_ev: f64,
    /// Lattice channel of `deposited_ev`: `phonon_emitted_ev + polaron_ev`.
    #[serde(default)]
    pub lattice_ev: f64,
    /// `|deposited + escaped + trapped + barrier - incident - fermi_sea| /
    /// (incident + fermi_sea)` (0 with no incident energy).
    pub relative_imbalance: f64,
}

/// Emission yields per primary.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct Yields {
    /// Backscatter yield `η`: fast electrons out of the front face per
    /// primary.
    pub backscatter_eta: f64,
    /// Secondary yield `δ`: slow electrons out of the front face per primary.
    pub secondary_delta: f64,
    /// Total yield `σ = η + δ`.
    pub total_sigma: f64,
    /// Fast electrons out of the back face per primary.
    pub transmitted_fast: f64,
    /// Slow electrons out of the back face per primary.
    pub transmitted_slow: f64,
}

/// Electrons leaving through one face.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct FaceEmission {
    /// Electrons.
    pub count: u64,
    /// Electrons per primary.
    pub per_primary: f64,
    /// Sum of their energies, eV.
    pub energy_ev: f64,
    /// Energy spectrum, eV, both classes.
    pub energy_histogram: Histogram,
    /// Slow electrons (`E < split`, the secondary class).
    pub slow: EmissionClass,
    /// Fast electrons (`E >= split`, the backscattered class).
    pub fast: EmissionClass,
}

/// One energy class of the electrons leaving through a face.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct EmissionClass {
    /// Electrons.
    pub count: u64,
    /// Electrons per primary.
    pub per_primary: f64,
    /// Sum of their energies, eV.
    pub energy_ev: f64,
    /// Mean energy, eV (0 if none).
    pub mean_energy_ev: f64,
    /// Polar angle from the outward surface normal, rad.
    pub polar_histogram: Histogram,
}

/// Deposited energy by layer and on the configured grids, eV.
///
/// Deposits are measured from the band bottom. With a secondary model in a
/// layer with a Fermi energy they include the Fermi-sea energy of liberated
/// conduction electrons, so they can sum to more than the energy the beam
/// imparted; see [`ElectronEnergyBudget::deposited_ev`] and
/// [the module docs](self#deposits-and-the-energy-imparted).
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct DepositionReport {
    /// By the layer the deposit was made in.
    pub per_layer_ev: Vec<f64>,
    /// The ionisation part of `per_layer_ev`, eV (electronic excitation,
    /// including the residual of sub-cutoff electrons; see
    /// [the channel split](self#the-ionisation-and-lattice-channels)). Empty
    /// in reports written before it existed.
    #[serde(default)]
    pub per_layer_ionisation_ev: Vec<f64>,
    /// The lattice part of `per_layer_ev`, eV (LO-phonon emission and polaron
    /// traps); `ionisation + lattice = per_layer_ev` per layer.
    #[serde(default)]
    pub per_layer_lattice_ev: Vec<f64>,
    /// On the Cartesian grid, if configured.
    pub cartesian: Option<CartesianDeposition>,
    /// On the cylindrical grid, if configured.
    pub cylindrical: Option<CylindricalDeposition>,
    /// Radial profile in the PSF depth slab, with per-history errors, if
    /// configured ([`ElectronTallyConfig::psf`]); fit it with
    /// [`super::psf::fit_psf`].
    #[serde(default)]
    pub psf: Option<RadialProfile>,
    /// Exposure statistics, if configured ([`ElectronTallyConfig::exposure`]).
    #[serde(default)]
    pub exposure: Option<ExposureReport>,
}

/// Deposited energy on a [`CartesianGrid`].
///
/// Measured from the band bottom: with a secondary model in a layer with a
/// Fermi energy the cells include the Fermi-sea energy of liberated
/// conduction electrons, so they can sum to more than the energy the beam
/// imparted (see [`ElectronEnergyBudget::deposited_ev`]).
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct CartesianDeposition {
    /// The grid.
    pub grid: CartesianGrid,
    /// Energy per voxel, eV, at [`CartesianGrid::index`].
    pub energy_ev: Vec<f64>,
    /// Energy deposited outside the grid, eV.
    pub outside_ev: f64,
    /// Ionisation part of `energy_ev` per voxel, eV; `ionisation + lattice =
    /// energy` per voxel. Empty in reports written before it existed.
    #[serde(default)]
    pub ionisation_ev: Vec<f64>,
    /// Lattice part of `energy_ev` per voxel, eV.
    #[serde(default)]
    pub lattice_ev: Vec<f64>,
    /// Ionisation part of `outside_ev`, eV.
    #[serde(default)]
    pub outside_ionisation_ev: f64,
    /// Lattice part of `outside_ev`, eV.
    #[serde(default)]
    pub outside_lattice_ev: f64,
}

impl CartesianDeposition {
    /// Energy in voxel `(ix, iy, iz)`, eV.
    pub fn at(&self, ix: usize, iy: usize, iz: usize) -> f64 {
        self.energy_ev[self.grid.index(ix, iy, iz)]
    }

    /// Energy per volume in each voxel, eV/m³.
    pub fn density_ev_per_m3(&self) -> Vec<f64> {
        let v = self.grid.voxel_volume_m3();
        self.energy_ev.iter().map(|e| e / v).collect()
    }
}

/// Deposited energy on a [`CylindricalGrid`].
///
/// Measured from the band bottom: with a secondary model in a layer with a
/// Fermi energy the cells include the Fermi-sea energy of liberated
/// conduction electrons, so they can sum to more than the energy the beam
/// imparted (see [`ElectronEnergyBudget::deposited_ev`]).
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct CylindricalDeposition {
    /// The grid.
    pub grid: CylindricalGrid,
    /// Energy per cell, eV, at [`CylindricalGrid::index`].
    pub energy_ev: Vec<f64>,
    /// Energy deposited outside the grid, eV.
    pub outside_ev: f64,
    /// Ionisation part of `energy_ev` per cell, eV; `ionisation + lattice =
    /// energy` per cell. Empty in reports written before it existed.
    #[serde(default)]
    pub ionisation_ev: Vec<f64>,
    /// Lattice part of `energy_ev` per cell, eV.
    #[serde(default)]
    pub lattice_ev: Vec<f64>,
    /// Ionisation part of `outside_ev`, eV.
    #[serde(default)]
    pub outside_ionisation_ev: f64,
    /// Lattice part of `outside_ev`, eV.
    #[serde(default)]
    pub outside_lattice_ev: f64,
}

impl CylindricalDeposition {
    /// Energy in cell `(ir, ix)`, eV.
    pub fn at(&self, ir: usize, ix: usize) -> f64 {
        self.energy_ev[self.grid.index(ir, ix)]
    }

    /// Energy per volume in each cell, eV/m³.
    pub fn density_ev_per_m3(&self) -> Vec<f64> {
        self.energy_ev
            .iter()
            .enumerate()
            .map(|(i, e)| e / self.grid.cell_volume_m3(i / self.grid.depth.bins))
            .collect()
    }
}

/// Energy-weighted moments of the deposition position (weights: the energy
/// of each deposit). Lengths in metres.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct GenerationVolume {
    /// Total weight: the deposited energy, eV.
    pub energy_ev: f64,
    /// Weighted mean of `[x, y, z]`.
    pub mean_m: [f64; 3],
    /// Weighted standard deviation of `[x, y, z]`.
    pub std_dev_m: [f64; 3],
    /// Weighted mean distance from the beam axis.
    pub mean_radius_m: f64,
    /// Weighted standard deviation of the distance from the beam axis.
    pub radius_std_dev_m: f64,
    /// Weighted RMS distance from the beam axis, `<y² + z²>^(1/2)`.
    pub rms_radius_m: f64,
}

/// Where electrons fell below the stopping threshold (unweighted, one sample
/// each).
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct StoppingPoints {
    /// Electrons that fell below the stopping threshold.
    pub stopped: u64,
    /// Depth moments, m; `None` with fewer than two.
    pub depth: Option<MomentSummary>,
    /// Moments of the distance from the beam axis, m.
    pub radial: Option<MomentSummary>,
    /// The primaries among `stopped`: those that fell below the threshold
    /// before any secondary of their history began (the transport follows a
    /// history's primary to its end first), i.e. histories whose fate is
    /// [`Fate::Stopped`]. Absent from reports written before it existed.
    #[serde(default)]
    pub primaries: PrimaryStoppingPoints,
}

/// Where the primaries fell below the stopping threshold (unweighted, one
/// sample per stopped primary): the penetration depth of stopped primaries,
/// as compared with other codes in `docs/validation.md`.
#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
pub struct PrimaryStoppingPoints {
    /// Primaries that fell below the stopping threshold.
    pub stopped: u64,
    /// Depth moments, m; `None` with fewer than two.
    pub depth: Option<MomentSummary>,
    /// Moments of the distance from the beam axis, m.
    pub radial: Option<MomentSummary>,
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn weighted_moments_merge_matches_direct() {
        let xs = [(1.0, 2.0), (3.0, 1.0), (-2.0, 0.5), (7.0, 4.0), (0.5, 3.0)];
        let mut all = WeightedMoments::default();
        for &(x, w) in &xs {
            all.push(x, w);
        }
        let mut a = WeightedMoments::default();
        let mut b = WeightedMoments::default();
        for &(x, w) in &xs[..2] {
            a.push(x, w);
        }
        for &(x, w) in &xs[2..] {
            b.push(x, w);
        }
        a.merge(&b);
        let wsum: f64 = xs.iter().map(|p| p.1).sum();
        let mean = xs.iter().map(|p| p.0 * p.1).sum::<f64>() / wsum;
        let var = xs.iter().map(|p| p.1 * (p.0 - mean).powi(2)).sum::<f64>() / wsum;
        for m in [all, a] {
            assert!((m.weight - wsum).abs() < 1e-12);
            assert!((m.mean - mean).abs() < 1e-12);
            assert!((m.variance() - var).abs() < 1e-12);
        }
        // Zero weights are ignored.
        let before = all;
        all.push(100.0, 0.0);
        assert_eq!(all, before);
    }

    #[test]
    fn grid_indices() {
        let g = CartesianGrid {
            x: Binning::new(0.0, 1.0, 2).unwrap(),
            y: Binning::new(-1.0, 1.0, 3).unwrap(),
            z: Binning::new(-1.0, 1.0, 4).unwrap(),
        };
        assert_eq!(g.voxels(), Some(24));
        assert_eq!(g.locate([0.75, 0.9, -0.9]), Some(g.index(1, 2, 0)));
        assert_eq!(g.index(1, 2, 0), 20);
        assert_eq!(g.locate([1.5, 0.0, 0.0]), None);
        let c = CylindricalGrid {
            r: Binning::new(0.0, 2.0, 2).unwrap(),
            depth: Binning::new(0.0, 1.0, 5).unwrap(),
        };
        assert_eq!(c.locate([0.3, 0.6, 0.8]), Some(c.index(1, 1)));
        let total: f64 = (0..2).map(|i| c.cell_volume_m3(i)).sum::<f64>() * 5.0;
        assert!((total - std::f64::consts::PI * 4.0).abs() < 1e-12);
    }
}